Wrapper for SB10MD

[giannibi]

Dear all,

I am currently teaching a robust control course at the University of Siena, Italy, and I am a happy Slycot user. It would be of utmost importance to me to have a wrapper for SB10MD (which implements the D step in the D-K iteration process in mu-synthesis), which in turn uses AB13MD, which as far as I can tell is wrapped already. I tried to no avail to do it myself but apparently I am lacking the required expertise.
I would be very grateful for any advice on the matter, or better for the inclusion of SB10MD in the near future.

Thanks in advance and keep up the good work,
Gianni

[roryyorke]

I had a go at this here: https://github.com/roryyorke/Slycot/tree/rory/sb10md

It's incomplete, and I can't promise any more time on it. Anyone is free to take over and submit.

I tried it on Skogestad&Postlethwaite first edition section 8.12.4, but couldn't get the same result. SB01MD gave out the same $\mu$ values for the first iteration, but not the same $D$; and this $D$ gave a worse $\mu$, so I must have done something wrong. Script below.

Besides this wrong result, there are two big outstanding todos:

• doc string
• unit test (could perhaps be based on S&P 8.12.4. example I've tried -- but we'd need to decouple it from python-control. Shouldn't be too hard: output the ABCD from the interconnect, and hard-code them)

Two things I noticed while doing this:

• for f2py, I needed to "break" the dependency cycle between a and lda (etc.). I tried to declare a as dimension(*) but that didn't work; maybe dimension(*,*) ? (is that meaningful?)
• conda-forge (or baseline?) python-control seems to depend on slycot; I tried mamba remove --force slycot, but it still wanted to remove python-control. This is probably some misunderstanding on my part, but regardless I think we need some instruction on installing python-control with development Slycot.

# all references are to:
#  Skogestad & Postelthwaite, Multivariable Feedback Control, first edition (1996)
import sys

# adapt to your system
if '/home/rory/src/slycot' not in sys.path[:2]:
    sys.path.insert(0, '/home/rory/src/slycot')

import matplotlib.pyplot as plt
from itertools import product
from slycot import sb10ad, sb10md

import numpy as np
import control as ct

# eq (8.143)
ggain = ct.StateSpace([],[],[],np.array([[87.8, -86.4],
                                         [108.2, -109.6]]))
gdyn = ct.tf2ss(ct.tf([1],[75,1]))

g = ct.append(gdyn,gdyn) * ggain
# we'll call the inputs to g 'up', for "u-prime"; and 'yp' for output
g.input_labels = ['up[0]','up[1]']
g.output_labels = ['yp[0]','yp[1]']

# eq (1.132)
widyn = ct.tf2ss(ct.tf([1,0.2],[0.5,1]))
wi = ct.append(widyn, widyn)
wi.output_labels = ['yd[0]', 'yd[1]']

#wpdyn = ct.tf2ss(ct.tf([0.5,0.05],[1,1e-3])) # table 8.2
wpdyn = 0.5*ct.tf2ss(ct.tf([10,1],[10,1e-5])) # table 8.2
wp = ct.append(wpdyn, wpdyn)
wp.input_labels = ['y[0]', 'y[1]']
wp.output_labels = ['z[0]', 'z[1]']

# fig 8.7, eq. (8.29)

# summing junction into G
s1 = ct.summing_junction(inputs=['u','ud'], output='up', dimension=2)

# summing junction to y
s2 = ct.summing_junction(inputs=['yp','w'], output='y', dimension=2)

# negative feedback: v = -y
fbk = ct.StateSpace([],[],[],-np.eye(2))
fbk.input_labels=['y[0]','y[1]']
fbk.output_labels=['v[0]','v[1]']

# put it all together

P = ct.interconnect([g, wi, wp, s1, s2, fbk],
                    inputs=['ud[0]','ud[1]','w[0]','w[1]','u[0]','u[1]'],
                    outputs=['yd[0]','yd[1]','z[0]','z[1]','v[0]','v[1]'])


# n = np.size(P.A, 0)
# m = np.size(P.B, 1)
# np_ = np.size(P.C, 0)
ncon = 2
nmeas = 2

k, cl0, gamma0, rcond = ct.hinfsyn(P, 2, 2)

# gamma0 = 100
# out = sb10ad(n, m, np_, ncon, nmeas, gamma0, P.A, P.B, P.C, P.D)

print(f'{gamma0=:.3f}')

# D-step
f = 2
order = 4
nblock = np.array([1,1,2])
itype = np.array([2,2,2])
qutol = 2

omega = np.geomspace(1e-3, 1e3, 5*100 + 1)
ordout, aout, bout, cout, dout, totord, ad, bd, cd, dd, mju, dwork = sb10md(f, order, nblock, itype, qutol, cl0.A, cl0.B, cl0.C, cl0.D, omega)

Dsys = ct.minreal(ct.StateSpace(ad, bd, cd, dd))

resp = Dsys(1j * omega)

# find DPD^-1
def invss(g):
    # find inverse of state-space system
    # eq (4.27)
    # inverse of d
    dinv = np.linalg.inv(g.D)
    # inverse system matrices
    ainv = g.A - g.B @ dinv @ g.C
    binv = g.B @ dinv
    cinv = -dinv @ g.C

    return ct.StateSpace(ainv, binv, cinv, dinv)

DPD = Dsys * P * invss(Dsys)

respinv = invss(Dsys)(1j * omega)

fig, ax = plt.subplots()
ax.loglog(omega, abs(resp[0,0]), lw=5, label='d1')
ax.loglog(omega, abs(respinv[0,0]), label='d1inv')
ax.legend()
fig.tight_layout()

def hinf_conda3(sys, omega, ncon, nmeas):
    # section 9.3.1
    A = sys.A
    B2 = sys.B[:,-ncon:]
    C1 = sys.C[:-nmeas]
    D12 = sys.D[:-nmeas,-ncon:]
    return np.vstack([
        np.hstack([A - 1j*omega*np.eye(A.shape[0]), B2]),
        np.hstack([C1 , D12]),
        ])
    return A,B2,C1,D12
    
k, cl1, gamma1, rcond = ct.hinfsyn(DPD, 2, 2)
print(gamma1)

ordout, aout, bout, cout, dout, totord, ad, bd, cd, dd, mju2, dwork = sb10md(f, order, nblock, itype, qutol, cl1.A, cl1.B, cl1.C, cl1.D, omega)

# cf. Fig 8.17
fig, ax = plt.subplots()
ax.loglog(omega, mju, label='iter1')
ax.loglog(omega, mju2, label='iter2')
ax.legend()
fig.tight_layout()

plt.show()

This is $\mu$ for the first two iterations:

And $d_1$ and $d_1^{-1}$ for the first iteration:

[giannibi]

Dear Rory,

first off, thank you so very much for your effort and your prompt response.

I was able to exactly reproduce your results. However, in order to build the package I had to merge your diffs in the master branch of Slycot and build from there. I could not build from your fork because the compilation of SB04OD.f failed with the following errors:

---8<---
[ 66%] Building Fortran object slycot/CMakeFiles/_wrapper.dir/src/SB04OD.f.o
/home/giannibi/Slycot/slycot/src/SB04OD.f:630:29:

570 |      $               D, LDD, SDIM, DWORK, DWORK(M+1), DWORK(2*M+1), P, LDP, Q,
    |                     2

......
630 | $ LDE, SDIM, DWORK, DWORK(N+1), DWORK(2*N+1), U, LDU, V,
| 1
Error: Type mismatch between actual argument at (1) and actual argument at (2) (INTEGER(4)/REAL(8)).
/home/giannibi/Slycot/slycot/src/SB04OD.f:630:35:

570 |      $               D, LDD, SDIM, DWORK, DWORK(M+1), DWORK(2*M+1), P, LDP, Q,
    |                        2

......
630 | $ LDE, SDIM, DWORK, DWORK(N+1), DWORK(2*N+1), U, LDU, V,
| 1
Error: Type mismatch between actual argument at (1) and actual argument at (2) (REAL(8)/INTEGER(4)).
/home/giannibi/Slycot/slycot/src/SB04OD.f:630:42:

570 |      $               D, LDD, SDIM, DWORK, DWORK(M+1), DWORK(2*M+1), P, LDP, Q,
    |                             2

......
630 | $ LDE, SDIM, DWORK, DWORK(N+1), DWORK(2*N+1), U, LDU, V,
| 1
Error: Type mismatch between actual argument at (1) and actual argument at (2) (REAL(8)/INTEGER(4)).
/home/giannibi/Slycot/slycot/src/SB04OD.f:630:71:

570 |      $               D, LDD, SDIM, DWORK, DWORK(M+1), DWORK(2*M+1), P, LDP, Q,
    |                                                      2

......
630 | $ LDE, SDIM, DWORK, DWORK(N+1), DWORK(2*N+1), U, LDU, V,
| 1
Error: Type mismatch between actual argument at (1) and actual argument at (2) (INTEGER(4)/REAL(8)).
/home/giannibi/Slycot/slycot/src/SB04OD.f:631:55:

571 |      $               LDQ, DWORK(3*M+1), LDWORK-3*M, 0, INFO )
    |                          2

......
631 | $ LDV, DWORK(3N+1), LDWORK-3N, 0, INFO )
| 1
Error: Type mismatch between actual argument at (1) and actual argument at (2) (INTEGER(4)/REAL(8)).
---8<---

Merging the diffs into the master branch worked ok and I was able to manually install the built package into conda by first installing the conda-forge version of python-control and slycot and then manually installing the local version of slycot.

I am now trying to investigate the mismatch between the results in S&P and the output of your test script. I will let you know if there is any progress on this.

Thank you once more for your precious help.

Sincerely,
Gianni

[roryyorke]

However, in order to build the package I had to merge your diffs in the master branch of Slycot and build from there. I could not build from your fork because the compilation of SB04OD.f failed with the following errors:

That's strange - I just checked, and I branched off 9485d94 which is current master head. I think to affect the Fortran build I would need to have changed CMakeLists.txt, which I didn't.

You got it working, which is the main thing.

[giannibi]

Did a few more experiments. I'm attaching two more scripts: the first is a rework of the S&P example, and the second implements the example seen at https://it.mathworks.com/help/robust/ref/uss.musyn.html.

import matplotlib.pyplot as plt
from itertools import product
from slycot import sb10ad, sb10md

import numpy as np
import control as ct

# Combine a matrix of transfer functions in a MIMO one
def combine(systems):
    rrows=[]
    for srow in systems:
        s1 = srow[0]
        if not isinstance(s1,ct.StateSpace):
            s1=ct.ss(s1)            
            
        for s2 in srow[1:]:
            if not isinstance(s2, ct.StateSpace):
                s2 = ct.ss(s2)
           
            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=ct.StateSpace(A,B,C,D,s1.dt)
        rrows.append(s1)
    r1=rrows[0]
    for r2 in rrows[1:]:      
        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=ct.StateSpace(A,B,C,D,r1.dt)
    return r1

# For computing the inverse of scalings
def invss(d):
    # inverse of D
    dinv = np.linalg.inv(d.D)
    # inverse system matrices
    ainv = d.A - d.B @ dinv @ d.C
    binv = d.B @ dinv
    cinv = -dinv @ d.C

    return ct.StateSpace(ainv, binv, cinv, dinv)


numP = [[[87.8], [-86.4]],
        [[108.2], [-109.6]]]
den = [75, 1]
denP = [[den, den],
        [den, den]]
P = ct.tf(numP, denP)
print("Plant transfer function: ", P)

Id = ct.tf(ct.ss([],[],[],np.eye(2)))
Zd = ct.tf(ct.ss([],[],[],np.zeros([2,2])))

# Uncertainty weight
widyn = ct.tf([1,0.2],[0.5,1])
wi = ct.tf(ct.append(widyn, widyn))
print("Uncertainty weight: ", wi)

# Performance weight
wpdyn = 0.5*ct.tf([10,1],[10,1e-5])
wp = ct.tf(ct.append(wpdyn, wpdyn))
print("Performance weight: ", wp)

# Make the LFT
G_ = np.block([ [Zd, Zd, wi],
                [wp*P, wp, wp*P],
                [-P, -Id, -P] ])
G = ct.minreal(combine(G_))
print("G" ,G)

# Controller I/O
nu = ny = f = 2
# Uncertainty structure and type
nblock = np.array([1,1,2])
itype = np.array([2,2,2])
# Scaling computation parameters
qutol = 2.0
order = 4
# Frequency range
omega = np.logspace(-3, 3, 61)


DPDI = G
gammaold = np.infty
epsi = 0.01
maxiter = 10

for i in range(maxiter):
    # K-step
    k, cl0, gamma, rcond = ct.hinfsyn(DPDI, ny, nu)
    print(f'{gamma=:.3f}')

    if np.abs(gammaold-gamma) <= epsi:
        break

    # D-step
    ordout, aout, bout, cout, dout, totord, ad, bd, cd, dd, mju, dwork = sb10md(f, order, nblock, itype, qutol, cl0.A, cl0.B, cl0.C, cl0.D, omega)
    Dsys = ct.minreal(ct.StateSpace(ad, bd, cd, dd))
    #resp = Dsys(1j * omega)
    DPDI =  Dsys * G * invss(Dsys)
    #respinv = invss(Dsys)(1j * omega)

    gammaold = gamma

fig, ax = plt.subplots()
ax.loglog(omega, mju, label='Final mu')
ax.legend()
fig.tight_layout()

plt.show()

import matplotlib.pyplot as plt
from itertools import product
from slycot import sb10ad, sb10md

import numpy as np
import control as ct

# Combine a matrix of transfer functions in a MIMO one
def combine(systems):
    rrows=[]
    for srow in systems:
        s1 = srow[0]
        if not isinstance(s1,ct.StateSpace):
            s1=ct.ss(s1)            
            
        for s2 in srow[1:]:
            if not isinstance(s2, ct.StateSpace):
                s2 = ct.ss(s2)
           
            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=ct.StateSpace(A,B,C,D,s1.dt)
        rrows.append(s1)
    r1=rrows[0]
    for r2 in rrows[1:]:      
        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=ct.StateSpace(A,B,C,D,r1.dt)
    return r1

# For computing the inverse of scalings
def invss(d):
    # inverse of D
    dinv = np.linalg.inv(d.D)
    # inverse system matrices
    ainv = d.A - d.B @ dinv @ d.C
    binv = d.B @ dinv
    cinv = -dinv @ d.C

    return ct.StateSpace(ainv, binv, cinv, dinv)



P = ct.tf(1, [1,-1])
print("Plant transfer function: ", P)

Id = ct.tf(ct.ss([],[],[],np.eye(1)))
Zd = ct.tf(ct.ss([],[],[],np.zeros([1,1])))

# Uncertainty weight
wi = 0.25*ct.tf([1/2,1], [1/32,1])
print("Uncertainty weight: ", wi)

# Performance weight
wp = 0.25*ct.tf([100],[1,0.5])
print("Performance weight: ", wp)

# Make the LFT
G_ = np.block([ [Zd, Zd, wi],
                [wp*P, wp, wp*P],
                [-P, -Id, -P] ])
G = ct.minreal(combine(G_))
print("G" ,G)

# Controller I/O
nu = ny = f = 1
# Uncertainty structure and type
nblock = np.array([1,1])
itype = np.array([2,2])
# Scaling computation parameters
qutol = 2.0
order = 4
# Frequency range
omega = np.logspace(-3, 3, 61)


DPDI = G
gammaold = np.infty
epsi = 0.01
maxiter = 10

for i in range(maxiter):
    # K-step
    k, cl0, gamma, rcond = ct.hinfsyn(DPDI, ny, nu)
    print(f'{gamma=:.3f}')
    
    if np.abs(gammaold-gamma) <= epsi:
        break

    # D-step
    ordout, aout, bout, cout, dout, totord, ad, bd, cd, dd, mju, dwork = sb10md(f, order, nblock, itype, qutol, cl0.A, cl0.B, cl0.C, cl0.D, omega)
    Dsys = ct.minreal(ct.StateSpace(ad, bd, cd, dd))
    #print(ct.tf(Dsys))
    #resp = Dsys(1j * omega)
    DPDI =  Dsys * G * invss(Dsys)
    #respinv = invss(Dsys)(1j * omega)

    gammaold = gamma

fig, ax = plt.subplots()
ax.loglog(omega, mju, label='Final mu')
ax.legend()
fig.tight_layout()

plt.show()

There can be mistakes on my side but in both examples it turns out that the achieved gamma in the K step behaves almost erratically across iterations. Moreover, it seems that only the "first" blocks of the scaling D at each iteration evaluate to something different from one. More specifically, in the first example the delta block corresponding to plant uncertainty is of size [1,1] and the delta block pertaining to performance is of size [2], while in the second they are both of size [1]. In both cases, the entries of the scaling system D related to the performance block do not seem to be "touched" by sb10md (i.e., they seem to be always the identity).
I will continue investigating.