GeezerVison is showing github.com — modal canonical form does not work with repeated eigenvalues · Issue #318 · python-control/python-control
modal canonical form does not work with repeated eigenvalues #318DescriptionThe modal_form() function appears to assume that all eigenvalues are distinct. The following example gives strange results:
import control as ct
sys_repeat = ct.StateSpace([[-1, 1], [0, -1]], [[0], [1]], [[1, 0]], [[0]])
sys_modal, T = ct.modal_form(sys_repeat)
print(sys_modal)
which yields
A = [[-1.00000000e+00 0.00000000e+00]
[ 2.22044605e-16 -1.00000000e+00]]
B = [[4.50359963e+15]
[4.50359963e+15]]
C = [[-1. 1.]]
D = [[0.]]
The problem is that the transformation T as computed in modal_form is singular:
[[-1.00000000e+00 1.00000000e+00]
[ 2.22044605e-16 0.00000000e+00]]
Activityrabraker commented on Jul 21, 2019 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. repagh commented on Aug 2, 2019 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____/_| [email protected]
<[___\_\_-----< t: +31 15 2628685
| o' mobile: +31 6 39846891 roryyorke commented on Dec 30, 2020 Some references:
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 zsysI 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)).
MetadataMetadataDevelopmentNo branches or pull requestsIssue actions
|