td04ad (convert tf to ss) returns incorrect poles #222

Activity

1. bnavigator commented on Oct 25, 2023

   

   bnavigator

   on Oct 25, 2023

   Collaborator

   More actions

   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.
2. roryyorke commented on Oct 28, 2023

   

   roryyorke

   on Oct 28, 2023

   Collaborator

   More actions

   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])}')