Repository navigation
Use SLICOT's TB05AD in Statespace.freqresp #116
Description
StateSpace.freqresp current converts to a transfer-function matrix (TFM) to evaluate the frequency response (see [1]) .
Polynomials can give numerical problems (see below), so this method should probably use [2] (unfortunately not currently exposed in Slycot). This technique uses an orthogonal similarity transform to make the system quicker to evaluate.
One possible problem is that if one evaluates a decoupled (diagonal TFM) system in which one system has poles on the imaginary axis, I imagine the whole frequency response will be infinite. I'm not sure if this is really important.
[1]
python-control/control/statesp.py
Line 398 in cdd3e73
[2] http://slicot.org/objects/software/shared/doc/TB05AD.html
A script to demonstrate the problem:
import numpy as np
import control
wn = 1;
zeta = 0.01;
g = control.tf(wn**2, [1,2*zeta*wn,wn**2])
def trivial_evalfr(gss,w):
I = np.eye(*gss.A.shape)
return np.dot(gss.C,np.linalg.solve(1j*w*I-gss.A,gss.B)) + gss.D
def trivpow(g,n):
from functools import reduce
from operator import mul
assert g.inputs == g.outputs
I = np.eye(g.inputs)
return reduce(mul,[g for i in range(n)])
gss = control.ss(g)
print()
print('{:3s} {:>9s} {:>9s} {:>9s} {:>9s}'.format('exp','tfv/ref','tfv/ref','ssv/ref','ssv/ref'))
print('{:3s} {:>9s} {:>9s} {:>9s} {:>9s}'.format('','1','dB','1','dB'))
exps = range(1,21)
gtfs = [trivpow(g,i)
for i in exps]
gsss = [trivpow(gss,i)
for i in exps]
for exp,gtfi,gssi in zip(exps,gtfs,gsss):
ref = g.evalfr(1)**exp
tfv = gtfi.evalfr(1)
ssv = trivial_evalfr(gssi,1)
tferr = np.abs(tfv/ref)[0,0]
sserr = np.abs(ssv/ref)[0,0]
db = lambda x: 20*np.log10(x)
print('{:3d} {:9.2e} {:9.2e} {:9.2e} {:9.2e}'.format(exp,tferr,db(tferr),sserr,db(sserr)))result is
exp tfv/ref tfv/ref ssv/ref ssv/ref
1 dB 1 dB
1 1.00e+00 0.00e+00 1.00e+00 -6.75e-15
2 1.00e+00 9.57e-13 1.00e+00 -1.54e-14
3 1.00e+00 -8.69e-12 1.00e+00 -2.41e-14
4 1.00e+00 -4.36e-08 1.00e+00 -3.18e-14
5 1.00e+00 1.10e-07 1.00e+00 -3.95e-14
6 1.00e+00 7.16e-05 1.00e+00 -4.82e-14
7 1.00e+00 1.65e-03 1.00e+00 -5.40e-14
8 9.01e-01 -9.08e-01 1.00e+00 -6.27e-14
9 2.88e-02 -3.08e+01 1.00e+00 -7.04e-14
10 2.51e-04 -7.20e+01 1.00e+00 -8.00e-14
11 2.27e-05 -9.29e+01 1.00e+00 -8.68e-14
12 7.62e-08 -1.42e+02 1.00e+00 -9.35e-14
13 6.38e-10 -1.84e+02 1.00e+00 -1.03e-13
14 1.38e-11 -2.17e+02 1.00e+00 -1.09e-13
15 4.16e-13 -2.48e+02 1.00e+00 -1.20e-13
16 1.76e-15 -2.95e+02 1.00e+00 -1.25e-13
17 1.13e-17 -3.39e+02 1.00e+00 -1.33e-13
18 2.74e-20 -3.91e+02 1.00e+00 -1.41e-13
19 6.49e-22 -4.24e+02 1.00e+00 -1.49e-13
20 1.72e-24 -4.75e+02 1.00e+00 -1.56e-13
Activity
rabraker commentedon Mar 29, 2017
Hi all,
I would like to work on this.
I agree that evaluation of the frequency response by converting to a transfer function is problematic.
I have in the last few days written a wrapper module that exposes two sets of functionality in slycots TB05AD. This needs to be tested further but the current version has been pushed to my github repo. It works for all the examples I've thrown at it though.
I've also been working modyfing the statespace class so that sys.freqresp uses TB05AD. Because this is pretty basic functionality, I think getting the frequency response should still work if somebody doesn't have slycot installed. Thus, I wrote a fallback method that emulates the functionality TB05AD in python. It works pretty much the same way using scipy.linalg.hessenberg().
Since both TB05AD and linalg.hessengberg() and linalg.solv() will all need to be called the same number of times, it's unclear to me if there is anything to be gained by using the slycot function...
I've never contributed to an open source project before. Do I just make a pull request when I think it's ready?
murrayrm commentedon Mar 29, 2017
That would be great if you want to take a cut at addressing this. We definitely want to make sure that have a fallback solution that doesn't require slycot, which will make the functionality much more accessible.
When you are ready, you can just issue a pull request and someone will take a look at the code and make sure it fits the programming style of the project, has appropriate unit tests, etc. We haven't put together a style guide for the project yet, but you can get a sense of common problems by looking at some of the recently merged pull requests.
ilayn commentedon Mar 29, 2017
I think the most convenient solution for this is given by Misra, Patel SIMAX 1988 Algo. 3.1 via staircase forms. If I'm not mistaken SLICOT is using a similar approach but Misra et al's solution is faster. I've tested the algorithm and it indeed provides a fast enough and accurate response.
roryyorke commentedon Mar 29, 2017
I'm not sure that numpy or scipy offer solvers that take advantage of the Hessenberg form (see, e.g. numpy.linalg and scipy.linalg). LAPACK itself apparently doesn't offer a Hessenberg-specific solver either, see e.g. this stackoverflow answer.
So, the non-Slycot fallback could be a straight solution of C*(sI-A)\B + D; anyone who cares enough about having it run faster can install Slycot.
ilayn commentedon Mar 29, 2017
The article I've linked to offers a solution precisely to that restriction, namely, you avoid solving a hessenberg form but instead, using the structure, only use LU decomposition updates which are fast to compute in any domain.
rabraker commentedon Mar 30, 2017
Roryyorke, thanks for pointing out that there isn't a Hessenberg specific solver in LAPACK or scipy. That hadn't actually occurred to me yet.
I've been working on this a bit this evening. For a random state space system with 50 states, 10 inputs and 5 outputs evaluated at 1000 frequency points, the slycot implementation takes about 0.1 seconds, while directly solving C(sI-A)\B takes about 0.6 seconds. Errors between the two are on the order 10^-11 for magnitude and 10^-14 for phases.
In the same scenario but with 500 states, slycot finishes in about 3 seconds while direct evaluation takes about 35 seconds. Errors between the two are on the order of 10^-9 for magnitude and 10^-12 for phase.
My plan is to finish writing unit tests for my current implementation. My personal feeling is that the method in the article linked by Ilayn is a good idea but it will take some time to understand and implement. I think that what I currently have will be acceptable to many users and it would be good to fix the robustness issue associated with the current implementation before optimizing for speed.
murrayrm commentedon Dec 30, 2017
@rabraker Are you still working on this?
rabraker commentedon Jan 2, 2018
murrayrm commentedon Jul 18, 2018
Fixed in PR #173.