modal canonical form does not work with repeated eigenvalues #318

Activity

1. rabraker commented on Jul 21, 2019

   

   rabraker

   on Jul 21, 2019

   Contributor

   More actions

   Just a couple thoughts on this.

   It looks like matlab's strategy to deal with repeated eigenvalues is to reduce the A matrix to block schur form. The MB03RD routine in slicot can do this, given a matrix in already in real schur form.

   So maybe the strategy to solve this would be to first wrap MB03RD over in slycot. We can use scipy.linalg.schur to put A into real schur form before the call to MB03RD.
2. repagh commented on Aug 2, 2019

   

   repagh

   on Aug 2, 2019 via email

   Member

   More actions

   I created a PR to add mb03rd to Slycot. Have fun with it.

   …

   On Sun, 21 Jul 2019 at 20:50, Arnold ***@***.***> wrote: Just a couple thoughts on this. It looks like matlab's strategy to deal with repeated eigenvalues is to reduce the A matrix to block schur form. The MB03RD <http://slicot.org/objects/software/shared/doc/MB03RD.html> routine in slicot can do this, given a matrix in already in real schur form. So maybe the strategy to solve this would be to first wrap MB03RD over in slycot. We can use scipy.linalg.schur to put A into real schur form before the call to MB03RD. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#318?email_source=notifications&email_token=AAVKIDNWPA6WU5VV44TCFULQASVYPA5CNFSM4H22AJ72YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD2OJJ4A#issuecomment-513578224>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAVKIDLB635XB2QZV3A3R3LQASVYPANCNFSM4H22AJ7Q> .

   -- René van Paassen | ______o____/_| Rene.vanPaassen@gmail.com <[___\_\_-----< t: +31 15 2628685 | o' mobile: +31 6 39846891
3. murrayrm

   added

   bug

   help wanted

   on Dec 26, 2020
4. murrayrm

   added this to the 0.9.x milestone on Dec 26, 2020
5. roryyorke commented on Dec 30, 2020

   

   roryyorke

   on Dec 30, 2020

   Contributor

   More actions

   Some references:

   • https://www.mathworks.com/help/control/ref/bdschur.html , in particular for the "condmax" parameter default
   • http://www.engr.mun.ca/~millan/Eng6825/canonicals.pdf - definition of canonical forms

   How important is the ordering of the eigenvalues in the result? If it is important, do we have eigenvalue sorting routines (I did a search for "sort" in the code, couldn't spot anything obvious) ?

   With this:

   def modal_form(xsys, pmax=None):
       """Convert a system into modal canonical form
   
       Parameters
       ----------
       xsys : StateSpace object
           System to be transformed, with state `x`
       pmax: float
           An upper bound on individual transformations.  Default eps ** -0.25.
   
       Returns
       -------
       zsys : StateSpace object
           System in modal canonical form, with state `z`
       T : matrix
           Coordinate transformation: z = T * x
   
       Notes
       -----
       pmax is passed to slycot.mb03rd.  The pmax value allows for a
       trade-off: larger values give smaller modal blocks, but
       worse-conditioned T.
       """
   
       if pmax is None:
           pmax = finfo(float64).eps ** -0.5
   
       aschur, tschur = schur(xsys.A)
       _, tmodal, _, _ = mb03rd(aschur.shape[0], aschur, tschur, pmax=pmax)
       return similarity_transform(xsys, tmodal, inverse=True), tmodal
   
   
   def similarity_transform(xsys, T, timescale=1, inverse=False):
       """Perform a similarity transformation, with option time rescaling.
   
       Transform a linear state space system to a new state space representation
       z = T x, where T is an invertible matrix.
   
       Parameters
       ----------
       xsys : StateSpace object
              System to transform
       T : 2D invertible array
           The matrix `T` defines the new set of coordinates z = T x.
       timescale : float
           If present, also rescale the time unit to tau = timescale * t
       inverse: boolean
           If True (default), transform so z = T x.  If False, transform
           so x = T z.
   
       Returns
       -------
       zsys : StateSpace object
           System in transformed coordinates, with state 'z'
   
       """
       # Create a new system, starting with a copy of the old one
       zsys = StateSpace(xsys)
   
       # Define a function to compute the right inverse (solve x M = y)
       def rsolve(M, y):
           return transpose(solve(transpose(M), transpose(y)))
   
       # Update the system matrices
       if not inverse:
           zsys.A = rsolve(T, dot(T, zsys.A)) / timescale
           zsys.B = dot(T, zsys.B) / timescale
           zsys.C = rsolve(T, zsys.C)
       else:
           zsys.A = solve(T, zsys.A).dot(T) / timescale
           zsys.B = solve(T, zsys.B) / timescale
           zsys.C = zsys.C.dot(T)
   
       return zsys

   I can find a modal decomposition of 1 / (s+1)**3 / (s**2 + s + 100) ** 2
   (triply-repeated eigenvalues at -1, and double repeated complex pole pair).

   from control import canonical_form, ss
   import numpy as np
   
   from scipy.linalg import schur
   from slycot import mb03rd
   from control import similarity_transform
   
   import control as ct
   
   s = ct.tf([1,0],[1])
   g_tf = 1 / (s+1)**3 * 1 / (s**2 + s + 100)**2
   
   g = ct.ss(g_tf)
   gmod, tmod = ct.modal_form(g)
   
   with np.printoptions(precision=3, suppress=True):
       print(gmod.A)
   
   print(f'{np.linalg.cond(tmod) = }')

   [[-0.5   13.936 49.989 -1.919  0.     0.    -0.   ]
    [-7.158 -0.5    3.282 -7.172  0.     0.     0.   ]
    [-0.    -0.    -0.5   -1.999  0.    -0.    -0.   ]
    [ 0.     0.    49.901 -0.5    0.    -0.     0.   ]
    [-0.    -0.     0.    -0.    -1.     0.758  0.388]
    [ 0.     0.    -0.     0.    -0.    -1.    -1.398]
    [ 0.     0.    -0.     0.    -0.     0.    -1.   ]]
   np.linalg.cond(tmod) = 13.060711659854526
   

   on pmax: cond(A B) <= cond(A) cond(B); the mb03rd docs state that the condition numbers on the order of pmax ** 2, so to get an upper bound of on cond(T) of 1/sqrt(eps), I've set pmax to 1/sqrt(sqrt(eps)).
6. murrayrm

   closed this as completedin #495on Jan 15, 2021