roryyorke
commented
12 months ago 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
aandlda(etc.). I tried to declareaasdimension(*)but that didn't work; maybedimension(*,*)? (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 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