René's URL Explorer Experiment

{"@context":"https://schema.org","@type":"DiscussionForumPosting","headline":"Use SLICOT's TB05AD in Statespace.freqresp","articleBody":"StateSpace.freqresp current converts to a transfer-function matrix (TFM) to evaluate the frequency response (see [1]) .  \n\nPolynomials 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.\n\nOne 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.\n\n[1] https://github.com/python-control/python-control/blob/cdd3e73838cea0b26a7155d0b2084aa698e5eefd/control/statesp.py#L398\n[2] http://slicot.org/objects/software/shared/doc/TB05AD.html\n\nA script to demonstrate the problem:\n\n``` python\nimport numpy as np\nimport control\n\nwn = 1;\nzeta = 0.01;\ng = control.tf(wn**2, [1,2*zeta*wn,wn**2])\n\ndef trivial_evalfr(gss,w):\n    I = np.eye(*gss.A.shape)\n    return np.dot(gss.C,np.linalg.solve(1j*w*I-gss.A,gss.B)) + gss.D\n\ndef trivpow(g,n):\n    from functools import reduce\n    from operator import mul\n    assert g.inputs == g.outputs\n    I = np.eye(g.inputs)\n    return reduce(mul,[g for i in range(n)])\n\ngss = control.ss(g)\n\nprint()\nprint('{:3s} {:\u003e9s} {:\u003e9s} {:\u003e9s} {:\u003e9s}'.format('exp','tfv/ref','tfv/ref','ssv/ref','ssv/ref'))\nprint('{:3s} {:\u003e9s} {:\u003e9s} {:\u003e9s} {:\u003e9s}'.format('','1','dB','1','dB'))\n\nexps = range(1,21)\n\ngtfs = [trivpow(g,i)\n        for i in exps]\ngsss = [trivpow(gss,i)\n        for i in exps]\n\nfor exp,gtfi,gssi in zip(exps,gtfs,gsss):\n    ref = g.evalfr(1)**exp\n    tfv = gtfi.evalfr(1)\n    ssv = trivial_evalfr(gssi,1)\n\n    tferr = np.abs(tfv/ref)[0,0]\n    sserr = np.abs(ssv/ref)[0,0]\n\n    db = lambda x: 20*np.log10(x)\n    print('{:3d} {:9.2e} {:9.2e} {:9.2e} {:9.2e}'.format(exp,tferr,db(tferr),sserr,db(sserr)))\n```\n\nresult is\n\n```\nexp   tfv/ref   tfv/ref   ssv/ref   ssv/ref\n            1        dB         1        dB\n  1  1.00e+00  0.00e+00  1.00e+00 -6.75e-15\n  2  1.00e+00  9.57e-13  1.00e+00 -1.54e-14\n  3  1.00e+00 -8.69e-12  1.00e+00 -2.41e-14\n  4  1.00e+00 -4.36e-08  1.00e+00 -3.18e-14\n  5  1.00e+00  1.10e-07  1.00e+00 -3.95e-14\n  6  1.00e+00  7.16e-05  1.00e+00 -4.82e-14\n  7  1.00e+00  1.65e-03  1.00e+00 -5.40e-14\n  8  9.01e-01 -9.08e-01  1.00e+00 -6.27e-14\n  9  2.88e-02 -3.08e+01  1.00e+00 -7.04e-14\n 10  2.51e-04 -7.20e+01  1.00e+00 -8.00e-14\n 11  2.27e-05 -9.29e+01  1.00e+00 -8.68e-14\n 12  7.62e-08 -1.42e+02  1.00e+00 -9.35e-14\n 13  6.38e-10 -1.84e+02  1.00e+00 -1.03e-13\n 14  1.38e-11 -2.17e+02  1.00e+00 -1.09e-13\n 15  4.16e-13 -2.48e+02  1.00e+00 -1.20e-13\n 16  1.76e-15 -2.95e+02  1.00e+00 -1.25e-13\n 17  1.13e-17 -3.39e+02  1.00e+00 -1.33e-13\n 18  2.74e-20 -3.91e+02  1.00e+00 -1.41e-13\n 19  6.49e-22 -4.24e+02  1.00e+00 -1.49e-13\n 20  1.72e-24 -4.75e+02  1.00e+00 -1.56e-13\n```\n","author":{"url":"https://github.com/roryyorke","@type":"Person","name":"roryyorke"},"datePublished":"2016-10-09T19:40:27.000Z","interactionStatistic":{"@type":"InteractionCounter","interactionType":"https://schema.org/CommentAction","userInteractionCount":9},"url":"https://github.com/116/python-control/issues/116"}