td04ad (convert tf to ss) returns incorrect poles

[murrayrm]

A bug was reported in python-control issue #935 in which tf2ss returns the incorrect poles.

It looks like the underlying problem is in td04ad. The following code illustrates the bug:

from numpy import array, roots, sort
from numpy.linalg import eig
from numpy.testing import assert_allclose
from slycot import td04ad
from scipy.signal import tf2ss

# Test case #1 (works)

num = array([[[1, 1, 1, 1, 1]]])
den = array([[1, 2, 1, 2, 1]])
denorder = np.array([den.size - 1])

ssout = td04ad('C', 1, 1, denorder, den, num, tol=0)

tf_poles = sort(roots(den[0]))
ss_poles = sort(eig(ssout[1])[0])
assert_allclose(tf_poles, ss_poles)

# Test case #2 (doesn't work)

num = array([[[1.05351324e-02, 2.83623596e-01, 2.44891698e+00, 1.22011448e+01,
               5.34846570e+01, 1.42092310e+02, 2.24627039e+02, 2.54305034e+02,
               2.18224432e+02, 1.24102950e+02, 4.99448809e+01, 1.26651596e+01,
               2.02239050e+00, 1.06681479e-01, 1.56212741e-03, 8.90622718e-06,
               1.77205330e-08]]])
den = array([[1.00000000e+00, 6.41551095e+02, 8.40213016e+03, 5.55121354e+04,
              1.93428590e+05, 1.31800818e+05, 9.20183802e+04, 2.90211403e+04,
              4.54824418e+03, 4.08636727e+02, 1.72367976e+01, 8.91452881e-02,
              1.04426479e-04, 4.32051569e-08, 5.88752872e-12, 3.24913608e-16,
              6.34764375e-21]])
denorder = np.array([den.size - 1])

# Try scipy
A, B, C, D = tf2ss(num[0, 0], den[0])
tf_poles = sort(roots(den[0]))
ss_poles = sort(eig(A)[0])
assert_allclose(tf_poles, ss_poles)

# Now try slycot
ssout = td04ad('C', 1, 1, denorder, den, num, tol=0)
ss_poles = sort(eig(ssout[1])[0])
assert_allclose(tf_poles, ss_poles)

For the second system, the scipy-generated poles match, but some of the the slycot-generated poles are not the same (and one of them is unstable).

It's possible this is a numerical issue (the coefficients span 21 orders of magnitude!), but odd the SciPy gets it right and slycot doesn't. (SciPy uses an observable canonical form-based realization.)

[bnavigator]

At first I also would have dismissed as numerical issues. But then when I was tinkering around with different parameters for td04ad, I found that there is a bug in the wrapper: n is not supposed to be an input parameter.

Slycot/slycot/transform.py

Line 624 in 9485d94

out = _wrapper.td04ad_c(m,p,index,dcoeff,ucoeff,n,tol,ldwork)

Slycot/slycot/src/transform.pyf

Line 352 in 9485d94

integer intent(in,out) :: nr != sum(index_bn)

>>> print(_wrapper.td04ad_c.__doc__)
nr,a,b,c,d,info = td04ad_c(m,p,index_bn,dcoeff,ucoeff,nr,[tol,ldwork,overwrite_dcoeff,overwrite_ucoeff])

Wrapper for ``td04ad_c``.

Parameters
----------
m : input int
p : input int
index_bn : input rank-1 array('i') with bounds (m)
dcoeff : input rank-2 array('d') with bounds (max(1, m),*)
ucoeff : input rank-3 array('d') with bounds (max(1, max(m, p)),max(1, max(m, p)),*)
nr : input int

Other Parameters
----------------
overwrite_dcoeff : input int, optional
    Default: 0
overwrite_ucoeff : input int, optional
    Default: 0
tol : input float, optional
    Default: 0
ldwork : input int, optional
    Default: max(1,nr+max(nr,max(3*m,3*p)))

Returns
-------
nr : int
a : rank-2 array('d') with bounds (max(1, nr),max(1, nr))
b : rank-2 array('d') with bounds (max(1, nr),max(m, p))
c : rank-2 array('d') with bounds (max(1, max(m, p)),max(1, nr))
d : rank-2 array('d') with bounds (max(1, max(m, p)),max(m, p))
info : int

https://github.com/python-control/SLICOT-Reference/blob/979f39d7863628407b0f9cae6804efc2833849ab/src/TD04AD.f#L71-L73

C     NR      (output) INTEGER
C             The order of the resulting minimal realization, i.e. the
C             order of the state dynamics matrix A.

I don't know if this is related. I will try to find time to fix the wrapper in the coming days and see if that changes results.

[roryyorke]

I think this is numerical:

• the zeros from td04ad are closer to the roots of the original numerator than the zeros of the observable canonical representation, so td04ad is not worse in all respects
• reconstructing a denominator polynomial from td04ad 's roots gives small absolute errors, though the relative errors get large towards the low-degree coefficients; but then the coefficient values are getting very small too

I think td04ad is balancing to get a minimum realization for MIMO TF->SS conversion, where you'd expect to drop states. I don't know why minreal of the Scipy conversion is different from the output of td04ad, but I'm also not sure it's all that important.

This is a high degree polynomial, so numerical issues are not too surprising. It would be good to give a warning when the user attempts a fragile TF->SS conversion, but we'd need to know how to quantify the sensitivity.

Script output:

0 states have been removed from the model
0 states have been removed from the model

max scipy pole diff         0
max slycot pole diff        0.000578
max scipy minreal pole diff 0.000182

max scipy pole real val     -6.25e-05
max slycot pole diff        0.000415
max scipy minreal pole diff 9.03e-05

max scipy zero diff         9.59e-05
max slycot zero diff        8.6e-09
max scipy minreal zero diff 2.76e-05

poly & reconstructed coeff difference
[1.00e+00 6.42e+02 8.40e+03 5.55e+04 1.93e+05 1.32e+05 9.20e+04 2.90e+04 4.55e+03 4.09e+02 1.72e+01 8.91e-02 1.04e-04 4.32e-08 5.89e-12 3.25e-16 6.35e-21]
[0.00e+00 1.14e-13 1.82e-11 1.24e-10 1.25e-09 9.60e-10 5.82e-11 2.07e-10 3.37e-10 1.56e-11 2.05e-11 3.43e-12 7.31e-13 8.35e-14 7.93e-16 1.99e-17 8.55e-18]

Script:

from numpy import array, roots, sort, poly
from numpy.linalg import eigvals
from numpy.testing import assert_allclose
from slycot import td04ad
from scipy.signal import tf2ss
from control import zeros, ss, poles, minreal

num = array([[[1.05351324e-02, 2.83623596e-01, 2.44891698e+00, 1.22011448e+01,
               5.34846570e+01, 1.42092310e+02, 2.24627039e+02, 2.54305034e+02,
               2.18224432e+02, 1.24102950e+02, 4.99448809e+01, 1.26651596e+01,
               2.02239050e+00, 1.06681479e-01, 1.56212741e-03, 8.90622718e-06,
               1.77205330e-08]]])
den = array([[1.00000000e+00, 6.41551095e+02, 8.40213016e+03, 5.55121354e+04,
              1.93428590e+05, 1.31800818e+05, 9.20183802e+04, 2.90211403e+04,
              4.54824418e+03, 4.08636727e+02, 1.72367976e+01, 8.91452881e-02,
              1.04426479e-04, 4.32051569e-08, 5.88752872e-12, 3.24913608e-16,
              6.34764375e-21]])
denorder = array([den.size - 1])

# ref poles & zeros
ref_poles = sort(roots(den[0]))
ref_zeros = sort(roots(num[0,0]))

A0, B0, C0, D0 = tf2ss(num[0, 0], den[0])
scipy_poles = sort(eigvals(A0))
scipy_zeros = sort(zeros(ss(A0,B0,C0,D0)))
scipy_minreal_poles = sort(poles(minreal(ss(A0,B0,C0,D0))))
scipy_minreal_zeros = sort(zeros(minreal(ss(A0,B0,C0,D0))))

# on my system, sorting doesn't *quite* work; the real parts aren't equal
# maybe fixable with a sort with a tolerance for equality
if max(abs(ref_zeros-scipy_minreal_zeros)) > 0.1:
    scipy_minreal_zeros[[6, 7, -2,-1]] = scipy_minreal_zeros[[7, 6, -1,-2]]
    if max(abs(ref_zeros-scipy_minreal_zeros)) > 0.1:
        raise ValueError

nr,A1,B1,C1,D1 = td04ad('C', 1, 1, denorder, den, num, tol=0)
slycot_poles = sort(eigvals(A1))
slycot_zeros = sort(zeros(ss(A1,B1,C1,D1)))

print()
print(f'max scipy pole diff         {max(abs(ref_poles-scipy_poles)):.3g}')
print(f'max slycot pole diff        {max(abs(ref_poles-slycot_poles)):.3g}')
print(f'max scipy minreal pole diff {max(abs(ref_poles-scipy_minreal_poles)):.3g}')
print()

print(f'max scipy pole real val     {max(scipy_poles.real):.3g}')
print(f'max slycot pole diff        {max(slycot_poles.real):.3g}')
print(f'max scipy minreal pole diff {max(scipy_minreal_poles.real):.3g}')
print()

print(f'max scipy zero diff         {max(abs(ref_zeros-scipy_zeros)):.3g}')
print(f'max slycot zero diff        {max(abs(ref_zeros-slycot_zeros)):.3g}')
print(f'max scipy minreal zero diff {max(abs(ref_zeros-scipy_minreal_zeros)):.3g}')
print()

# reconstruct denominator from scipy poles
recon_slycot_den = poly(slycot_poles)

print('poly & reconstructed coeff difference')
print(den[0])
print(f'{abs(recon_slycot_den - den[0])}')