Elegant way to generate a transfer function matrix

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...

You could try using parallel, series and feedback to construct the system. Alternatively, take a look at the append and connect commands? 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.connect function 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 the connect function, but I still think it is confusing.

@murrayrm

Thanks for your help @murrayrm !
But I don't get it. Input of G is the controller output or better say the disturbed controller output.

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 append to stack the system dynamics and then use connect to 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

When 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-control are 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 D matrix is not correct.

The problem is because of a limitations in the connect command (in MATLAB and Python):

The outputs for the connect command 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 with connect (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)

@bnavigator

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.

@murrayrm

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 + 10

MATLAB:

>> 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 + 10

@bnavigator

I 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 outputs

But 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():
    pass

with 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 the connect() function and intermediate "distribution" and "summer" systems Ieand Iu.

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 lmi method to avoid the rank conditions but still always penalize your control output even the slightest to avoid the ill-conditioned solutions.

@ilyan

@ilyan do you have an example or a literature source for solving this as a lmi formulation?

@iljastas

@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.html

This 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.

Navigation Menu

Repository navigation

Elegant way to generate a transfer function matrix #325

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Elegant way to generate a transfer function matrix #325

Description

Activity

murrayrm commented on Jul 11, 2019

murrayrm commented on Jul 13, 2019

iljastas commented on Jul 15, 2019

murrayrm commented on Jul 15, 2019

iljastas commented on Jul 16, 2019

murrayrm commented on Jul 16, 2019

bnavigator commented on Jul 16, 2019

murrayrm commented on Jul 16, 2019

iljastas commented on Jul 16, 2019

bnavigator commented on Jul 16, 2019

iljastas commented on Jul 17, 2019

bnavigator commented on Jul 18, 2019

bnavigator commented on Jul 18, 2019

iljastas commented on Jul 19, 2019

ilayn commented on Jul 19, 2019

iljastas commented on Jul 20, 2019

ilayn commented on Jul 20, 2019

murrayrm commented on Dec 29, 2019

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions