Use SLICOT's TB05AD in Statespace.freqresp #116

roryyorke · 2016-10-09T19:40:27Z

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

# when evaluating at many frequencies, much faster to convert to

[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

The text was updated successfully, but these errors were encountered:

rabraker · 2017-03-29T04:27:19Z

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 · 2017-03-29T05:59:16Z

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 · 2017-03-29T07:43:12Z

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 · 2017-03-29T19:31:13Z

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 · 2017-03-29T19:51:08Z

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 · 2017-03-30T05:28:03Z

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 · 2017-12-30T16:23:00Z

@rabraker Are you still working on this?

rabraker · 2018-01-02T03:07:15Z

Thanks for the reminder, I had essentially forgotten about this. I think I have the function written. Should be able to submit a PR before winter break is over.

…

On 12/30/17, Richard Murray wrote: @rabraker Are you still working on this? -- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: #116 (comment)

murrayrm · 2018-07-18T04:06:12Z

Fixed in PR #173.

rabraker mentioned this issue Apr 18, 2017

TB05AD wrapper python-control/Slycot#10

Merged

murrayrm added the enhancement label Dec 30, 2017

This was referenced Jan 2, 2018

Use slycot's tb05ad for faster and more accurate frequency response e… #172

Closed

Use slycot's tb05ad for faster and more accurate frequency response #173

Merged

murrayrm closed this as completed Jul 18, 2018

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Use SLICOT's TB05AD in Statespace.freqresp #116

Use SLICOT's TB05AD in Statespace.freqresp #116

roryyorke commented Oct 9, 2016

rabraker commented Mar 29, 2017

murrayrm commented Mar 29, 2017

ilayn commented Mar 29, 2017 •

edited

Loading

roryyorke commented Mar 29, 2017

ilayn commented Mar 29, 2017

rabraker commented Mar 30, 2017

murrayrm commented Dec 30, 2017

rabraker commented Jan 2, 2018 via email

murrayrm commented Jul 18, 2018

Use SLICOT's TB05AD in Statespace.freqresp #116

Use SLICOT's TB05AD in Statespace.freqresp #116

Comments

roryyorke commented Oct 9, 2016

rabraker commented Mar 29, 2017

murrayrm commented Mar 29, 2017

ilayn commented Mar 29, 2017 • edited Loading

roryyorke commented Mar 29, 2017

ilayn commented Mar 29, 2017

rabraker commented Mar 30, 2017

murrayrm commented Dec 30, 2017

rabraker commented Jan 2, 2018 via email

murrayrm commented Jul 18, 2018

ilayn commented Mar 29, 2017 •

edited

Loading