Repository navigation
modal canonical form does not work with repeated eigenvalues #318
Description
Activity
Just a couple thoughts on this.
It looks like matlab's strategy to deal with repeated eigenvalues is to reduce the
Amatrix 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
MB03RDover in slycot. We can usescipy.linalg.schurto putAinto real schur form before the call toMB03RD.- 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
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.060711659854526on 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)).
The
modal_form()function appears to assume that all eigenvalues are distinct. The following example gives strange results:which yields
The problem is that the transformation
Tas computed inmodal_formis singular: