Repository navigation
Elegant way to generate a transfer function matrix #325
Description
Activity
You could try using
parallel,seriesandfeedbackto construct the system. Alternatively, take a look at theappendandconnectcommands? All of those should produce a system that doesn't involve creating extra states (unless you list one of the systems multiple times, which I guess is what is causing the problem above?).The following should work:
# Create a 3 input, 3 output system consisting of component systems sysvec = append(ss(W1), ss(W2), Gss) # Add interconnections between the components (NB: indexing starts at 1) # Input of G (3) is output of W1 (1), input of W2 (2) is G (3) # Keep all inputs and outputs P_ = connect(sysvec, [[3, 1], [2, 3]], [1, 2, 3], [1, 2, 3])Note that this generates a minimal representation in the state space because we never create a copy of the plant (G) or weighting filters. So
P_has 4 states, 3 inputs, 3 outputs.One annoying (and I think confusing) aspect of this construction is that the
control.connectfunction uses indexing starting at 1 (to be compatible with MATLAB but also to allow negative feedback), while Python indexing normally starts at 0. This is documented in theconnectfunction, but I still think it is confusing.I may have misunderstood the transfer function matrix you were trying to create (I sketched it out on a napkin that I have since tossed -:). The main point is that you can use
appendto stack the system dynamics and then useconnectto interconnect them.I'm unfortunately not able to produce a minimal example (Kai Müller, Entwurf robubster Regelungen Teubner Springer 1996, Chapter 12.3.1, p. 172).
from control import *
W1 = tf([2], [1])
W2 = tf([3], [1])
G = tf([1], [1, 1])
sysvec = append(ss(G), ss(W1), ss(W2))
P = connect(sysvec, [[1, 2] ], [2, 1], [2, 3, 1])
print(P)A = [[-1.]]
B = [[1. 0.]]
C = [[1.]
[0.]
[ 0.]]
D = [[0. 0.]
[0. 2.]
[0. 0.]]The right answer is: A=[-1], B = [0 1], C=[2 0 1]^T, D=[2 0; 0 3; 1 0]
But I'm not able to produce it.So there is a connection between G and W1. And I have to inputs: u and z, And I have three outputs: v1, v2 and y.
I put this example into MATLAB (same exact commands) and got the output
P = A = x1 x1 -1 B = u1 u2 x1 2 1 C = x1 y1 0 y2 0 y3 1 D = u1 u2 y1 2 0 y2 0 0 y3 0 0When I run those commands in Python 3.5 using the current master (which should be the same as v0.8.2), I get output
A = [[-1.]] B = [[2. 1.]] C = [[0.] [0.] [1.]] D = [[2. 0.] [0. 0.] [0. 0.]]This matches the output of MATLAB, but is different than what you showed as your output. What version of
python-controlare you running?There is still the issue that this is not the right system description for the block diagram. I think that part of the problem is that your interconnection matrix is not right. You have
[[1, 2]]for the interconnection, which sets input 1 (G) to output 2 (W1). From the block diagram, what you want is[[2, 1]], so that the input of W1 (= input #2) is the output from G (output #1).Running the command with the updated interconnection matrix gives
P = connect(sysvec, [[2, 1]], [2, 1], [2, 3, 1])gives
A = [[-1.]] B = [[0. 1.]] C = [[2.] [0.] [1.]] D = [[2. 0.] [0. 0.] [0. 0.]]This is almost what you expect, except that the
Dmatrix is not correct.The problem is because of a limitations in the
connectcommand (in MATLAB and Python):-
The outputs for the
connectcommand can only correspond to outputs from one of the subsystems. The output that you want for y in your diagram is the output of G plus the input z. But you can't actually do that withconnect(in MATLAB or Python). So one problem is that the input u should be going to output v2 through a gain of 3. (That entry is missing in both the MATLAB and Python versions, for reasons that don't make senses.) -
There is a similar constrain on the inputs. It is not possible to take an input from one system and map it to the input of another. So the only input to W2 is its natural input and we removed this from the transfer function.
-
Try this:
from control import tf,ss,connect,append UU = tf([1],[1]) #1 u ->u_o Z = tf([1],[1]) #2 z ->z_o G = tf([1], [1, 1]) #3 ug ->y W1 = tf([2], [1]) #4 uw1 ->v1 W2 = tf([3], [1]) #5 uw2 ->v2 ZY = tf([1],[1]) #6 zy_i->zy_o sysvec = append(ss(UU),ss(Z),ss(G),ss(W1),ss(W2),ss(ZY)) Q=[[6,2],[6,3], # zy_i=z_o+y [4,6], # uw1=zy_o [3,1], # ug = u_o [5,1], # uw2 =u_o ] P = connect(sysvec,Q,[2,1],[4,5,6]) # z,u -> v1,v2,z+y, print(P)
which gives your "right answer"
A = [[-1.]] B = [[0. 1.]] C = [[2.] [0.] [1.]] D = [[2. 0.] [0. 3.] [1. 0.]]I am pretty sure your D does not reflect your desired output of
(v1,v2,y)but it is(v1,v2,z+y)and your inputs are ordered(z,u)Thanks @bnavigator. It would be useful to think about whether there is a nice way to do these sorts of manipulations. The construct that @iljastas gave at the start is really the best way since it creates copies of the system/weighting filters that are then "minimized" away, but it is also tedious to put in the extra elements that you used to get things to work.
I'd lover to hear suggestions about what a cleaner implementation would look like.
Thanks @murrayrm and @bnavigator I will try it later or tomorrow and post also a more complex example with and hinfsyn-controller!
Maybe it would be useful to have a function
combine()or something similar to mimic the Matrix combination of systems like in MATLAB. Here is a hack how that could look like. Of course it still needs all the type checks etc. to have it included into the StateSpace class:import numpy as np from control import StateSpace,tf,minreal from control.statesp import _convertToStateSpace def combine(systems): """ systems: 2D array of systems to combine """ rrows=[] for srow in systems: s1 = srow[0] if not isinstance(s1,StateSpace): s1=_convertToStateSpace(s1) for s2 in srow[1:]: if not isinstance(s2, StateSpace): s2 = _convertToStateSpace(s2) if s1.dt != s2.dt: raise ValueError("Systems must have the same time step") n = s1.states + s2.states m = s1.inputs + s2.inputs p = s1.outputs if s2.outputs != p: raise ValueError('inconsistent systems') A = np.zeros((n, n)) B = np.zeros((n, m)) C = np.zeros((p, n)) D = np.zeros((p, m)) A[:s1.states, :s1.states] = s1.A A[s1.states:, s1.states:] = s2.A B[:s1.states, :s1.inputs] = s1.B B[s1.states:, s1.inputs:] = s2.B C[:, :s1.states] = s1.C C[:, s1.states:] = s2.C D[:, :s1.inputs] = s1.D D[:, s1.inputs:] = s2.D s1=StateSpace(A,B,C,D,s1.dt) rrows.append(s1) r1=rrows[0] for r2 in rrows[1:]: if r1.dt != r2.dt: raise ValueError("Systems must have the same time step") n = r1.states + r2.states m = r1.inputs if r2.inputs != m: raise ValueError('inconsistent systems') p = r1.outputs + r2.outputs A = np.zeros((n, n)) B = np.zeros((n, m)) C = np.zeros((p, n)) D = np.zeros((p, m)) A[:r1.states, :r1.states] = r1.A A[r1.states:, r1.states:] = r2.A B[:r1.states, :] = r1.B B[r1.states:, :] = r2.B C[:r1.outputs, :r1.states] = r1.C C[r1.outputs:, r1.states:] = r2.C D[:r1.outputs, :] = r1.D D[r1.outputs:, :] = r2.D r1=StateSpace(A,B,C,D,r1.dt) return r1 Gss = StateSpace([[-1.0, -10.0], [1, 0]], [[10.0], [0]], [0, 1], [0]) G = tf(Gss) W1 = tf([100], [0.1, 1]) W2 = tf([100], [0.1, 1]) W2G=W2*G P11=np.block([[0,0],[W2G,W2]]) P12=np.block([[W1],[W2G]]) P21=np.block([[G,1]]) P22=np.block([[G]]) P_=np.block([[P11,P12],[P21,P22]]) Pc=combine(P_) Pctf=minreal(tf(Pc))
output:
8 states have been removed from the model Input 1 to output 1: 0 - 1 Input 1 to output 2: 1e+04 ------------------------- s^3 + 11 s^2 + 20 s + 100 Input 1 to output 3: 10 ------------ s^2 + s + 10 Input 2 to output 1: 0 - 1 Input 2 to output 2: 1000 ------ s + 10 Input 2 to output 3: 1 - 1 Input 3 to output 1: 1000 ------ s + 10 Input 3 to output 2: 1.756e-12 s + 1e+04 ------------------------- s^3 + 11 s^2 + 20 s + 100 Input 3 to output 3: 10 ------------ s^2 + s + 10MATLAB:
>> G=ss([-1,-10;1,0],[10;0],[0,1],[0]); >> W1=tf(100,[0.1 1]); >> W2=W1; >> W2G=W2*G; >> P11 = [0 0; W2G W2]; >> P12 = [W1; W2G]; >> P21 = [G 1]; >> P22 = [G]; >> Pc = [P11 P12; P21 P22]; >> Pctf=tf(Pc) Pctf = From input 1 to output... 1: 0 10000 2: ------------------------- s^3 + 11 s^2 + 20 s + 100 10 3: ------------ s^2 + s + 10 From input 2 to output... 1: 0 1000 2: ------ s + 10 3: 1 From input 3 to output... 1000 1: ------ s + 10 10000 2: ------------------------- s^3 + 11 s^2 + 20 s + 100 10 3: ------------ s^2 + s + 10 Continuous-time transfer function.
A
combine()for TransferFunction would be even simpler, but that is just what has been done manually in the ininital question:In [113]: from control import * ...: Gss = ss([[-1.0, -10.0], [1, 0]], [[10.0], [0]], [0, 1], [0]) ...: G = tf(Gss) ...: W1 = tf([100], [0.1, 1]) ...: W2 = tf([100], [0.1, 1]) ...: W2G = W2*G ...: num = [ [ [0], [0], W1.num[0][0] ], [W2G.num[0][0], W2.num[0][0], G.num[0][0]], [G.num[0][0], [1], G.num[0][0] ] ] ...: den = [ [ [1], [1], W1.den[0][0] ], [W2G.den[0][0], W2.den[0][0], G.den[0][0]], [G.den[0][0], [1], G.den[0][0] ] ] ...: P = tf( num, den ) ...: print(minreal(P)) -2 states have been removed from the model Input 1 to output 1: 0 - 1 Input 1 to output 2: 1e+04 ------------------------- s^3 + 11 s^2 + 20 s + 100 Input 1 to output 3: 10 ------------ s^2 + s + 10 Input 2 to output 1: 0 - 1 Input 2 to output 2: 1000 ------ s + 10 Input 2 to output 3: 1 - 1 Input 3 to output 1: 1000 ------ s + 10 Input 3 to output 2: 10 ------------ s^2 + s + 10 Input 3 to output 3: 10 ------------ s^2 + s + 10I was not able to get the plant in a more complex example by the append, connect procedure.
Because I have a MIMO-System (G) with a disturance z2 after G: z2, G element of R^2 and W1 has to 2 inputs and two outputsBut thanks to @bnavigator: the combine-method is great. Thanks a lot! Maybe my example will help others.
import control from control import * import matplotlib.pyplot as plt import numpy as np from control.statesp import _convertToStateSpace
def combine(systems): rrows=[] for srow in systems: s1 = srow[0] if not isinstance(s1,StateSpace): s1=_convertToStateSpace(s1) for s2 in srow[1:]: if not isinstance(s2, StateSpace): s2 = _convertToStateSpace(s2) if s1.dt != s2.dt: raise ValueError("Systems must have the same time step") n = s1.states + s2.states m = s1.inputs + s2.inputs p = s1.outputs if s2.outputs != p: raise ValueError('inconsistent systems') A = np.zeros((n, n)) B = np.zeros((n, m)) C = np.zeros((p, n)) D = np.zeros((p, m)) A[:s1.states, :s1.states] = s1.A A[s1.states:, s1.states:] = s2.A B[:s1.states, :s1.inputs] = s1.B B[s1.states:, s1.inputs:] = s2.B C[:, :s1.states] = s1.C C[:, s1.states:] = s2.C D[:, :s1.inputs] = s1.D D[:, s1.inputs:] = s2.D s1=StateSpace(A,B,C,D,s1.dt) rrows.append(s1) r1=rrows[0] for r2 in rrows[1:]: if r1.dt != r2.dt: raise ValueError("Systems must have the same time step") n = r1.states + r2.states m = r1.inputs if r2.inputs != m: raise ValueError('inconsistent systems') p = r1.outputs + r2.outputs A = np.zeros((n, n)) B = np.zeros((n, m)) C = np.zeros((p, n)) D = np.zeros((p, m)) A[:r1.states, :r1.states] = r1.A A[r1.states:, r1.states:] = r2.A B[:r1.states, :] = r1.B B[r1.states:, :] = r2.B C[:r1.outputs, :r1.states] = r1.C C[r1.outputs:, r1.states:] = r2.C D[:r1.outputs, :] = r1.D D[r1.outputs:, :] = r2.D r1=StateSpace(A,B,C,D,r1.dt) return r1
def make_lowpass(dc, crossw): return tf([dc], [crossw, 1]) def make_highpass(hf_gain, m3db_frq): return tf([hf_gain, 0], [1, m3db_frq]) def make_weight(dc, crossw, hf): s = tf([1, 0], [1]) return (s * hf + crossw) / (s + crossw / dc)
v = 1.1 lo = 2.2 A = [ [0, v], [0, 0] ] B = [ [lo], [1] ] C = [ [1, 0], [0, 1] ] D = [ [0], [0]] Gss = ss(A, B, C, D) G = tf(Gss) G1 = G[0,0] G2 = G[1,0] Wy1 = make_weight(dc=1, crossw=2.1, hf=0) Wy2 = make_highpass(hf_gain=0.99, m3db_frq=2.1) Wu = make_highpass(hf_gain=0.3, m3db_frq=2.1) # P11 = [ [W1, 0, W1*G1], [0, W2, W2*G2], [0, 0, 0] ] # P12 = [ [W1*G1], [W2*G2], [Wu] ] # P21 = [ [-1, 0, -G1], [0, -1, -G2] ] # P22 = [ [-G1], [-G2] ] # P = [ [P11, P21], [P21, P22]] P11=np.block( [ [Wy1, 0, Wy1*G1], [0, Wy2, Wy2*G2], [0, 0,0] ] ) P12=np.block( [ [Wy1*G1], [Wy2*G2], [Wu] ]) P21=np.block( [ [-1, 0, -G1], [0, -1, -G2] ] ) P22=np.block( [ [-G1], [-G2] ] ) P_ = np.block( [ [P11, P12], [P21, P22]] ) Pc = combine(P_) print(Pc.A.shape) Pcss = minreal( Pc ) print(Pcss) K, CL, gam, rcond = hinfsyn(Pcss, 2, 1)
Please format your posts to include code within code tags. Mark the multiline code and click on the '<>' button. It will enclose the code in
```tags``` def yourcode(): pass ```results in
def yourcode(): passwith the python tag you even get some highlighting:
```python def yourcode(): pass ```def yourcode(): pass
See also https://guides.github.com/features/mastering-markdown/
Much better now :)
BTW have you looked at
control.robust.augw()? It does something similar using theconnect()function and intermediate "distribution" and "summer" systemsIeandIu.Of course I tried this with "control" and in MATLAB but I got the errror:
ValueError: The matrix [A-jomegaI, B1 ; C2, D21] does not have full row rank with respect to the tolerance eps.
for:
from control import * def make_weight(dc, crossw, hf): s = tf([1, 0], [1]) return (s * hf + crossw) / (s + crossw / dc) def merge_tfs(tf1, tf2): num = [ [tf1.num[0][0], [0] ], [ [0], tf2.num[0][0]] ] den = [ [tf1.den[0][0], [1] ], [ [1], tf2.den[0][0]] ] return tf( num, den ) v, lo = 1.1, 1.4 A = [ [0, v], [0, 0] ] B = [ [lo], [1] ] C = [ [1, 0], [0, 1] ] D = [ [0], [0] ] Gss = ss(A, B, C, D) w1y = make_weight(dc=100.0, crossw=0.1, hf=0) w1Y = make_weight(dc=0.9, crossw=0.5, hf=0) w1 = merge_tfs(w1y, w1Y) w2 = make_weight(dc=0.1, crossw=0.8, hf=1.0) print("self.Gss:\n", Gss) Paug = augw(g=Gss, w1=ss(w1) ) print(Paug) Kss, CL, gamma, rcond = hinfsyn(Paug, 1, 1) print("hinfsyn gamma:", gamma)
And further on my control system has only to control disturbances, so that this formulation is more suited to my problem
In matlab, use
lmimethod to avoid the rank conditions but still always penalize your control output even the slightest to avoid the ill-conditioned solutions.@iljastas Just have a look at the hinfsyn options. You can change the method. If nothing has changed in the last years, the default should be riccati based but you can change it to
lmi. See this page https://nl.mathworks.com/help/robust/ref/hinfsynoptions.htmlThis issue has been idle for ~6 months and I don't think there is anything actionable in the comments except some ideas for future functionality. Closing for now, but someone should create a new issue if there is something specific that we would like to request as an enhancement.
Hi everybody,
I want to know if there is a more elegant way to generate a transfer-function-matrix for a hinf-controller.
In matlab I just do:
P11 = [0 0; W2G W2]
P12 = [W1; W2G]
P21 = [G 1]
P22 = [G]
P_ = [P11 P12; P21 P22]
sys = ss(P_, 'minimal');
In python with control, I habe to do:
import control
from control import *
Gss = ss([[-1.0, -10.0], [1, 0]], [[10.0], [0]], [0, 1], [0])
G = tf(Gss)
W1 = tf([100], [0.1, 1])
W2 = tf([100], [0.1, 1])
W2G = W2*G
num = [ [ [0], [0], W1.num[0][0] ], [W2G.num[0][0], W2.num[0][0], G.num[0][0]], [G.num[0][0], [1], G.num[0][0] ] ]
den = [ [ [1], [1], W1.den[0][0] ], [W2G.den[0][0], W2.den[0][0], G.den[0][0]], [G.den[0][0], [1], G.den[0][0] ] ]
P = tf( num, den )
print(P)
The problem is that the return of the numerator and denominator is not directly the list.
And I get a high order system. Should I use just the balanced trunction?
Pss = ss(P)
Pbr = balred(Pss, 4)
Matlab is using minreal with exist also in "control", but I have to find manually a suitable tolerance...