You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
# Implementation inspired by scipy https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.cdf2rdf.html
22
+
# with the licence https://github.com/scipy/scipy/blob/v1.6.3/LICENSE.txt
23
+
D,V = E
24
+
n =length(D)
25
+
26
+
# get indices for each first pair of complex eigenvalues
27
+
complex_inds =complex_indices(D)
28
+
29
+
# eigvals are sorted so conjugate pairs are next to each other
30
+
j = complex_inds[1:2:end]
31
+
k = complex_inds[2:2:end]
32
+
33
+
# put real parts on diagonal
34
+
Db =zeros(n, n)
35
+
Db[diagind(Db)] .=real(D)
36
+
37
+
# compute eigenvectors for real block diagonal eigenvalues
38
+
U =zeros(eltype(D), n, n)
39
+
U[diagind(U)] .=1.0
40
+
41
+
# transform complex eigvals to real blockdiag form
42
+
for (k,j) inzip(k,j)
43
+
Db[j, k] =imag(D[j]) # put imaginary parts in blocks
44
+
Db[k, j] =imag(D[k])
45
+
46
+
U[j, j] =0.5im
47
+
U[j, k] =0.5
48
+
U[k, j] =-0.5im
49
+
U[k, k] =0.5
50
+
end
51
+
Vb =real(V*U)
52
+
53
+
return Db, Vb
54
+
end
55
+
56
+
"""
57
+
sysm, T = modal_form(sys; C1 = false)
58
+
59
+
Bring `sys` to modal form.
60
+
61
+
The modal form is characterized by being tridiagonal with the real values of eigenvalues of `A` on the main diagonal and the complex parts on the first sub and super diagonals. `T` is the similarity transform applied to the system such that
62
+
```julia
63
+
sysm ≈ similarity_transform(sys, T)
64
+
```
65
+
66
+
If `C1`, then an additional convention for SISO systems is used, that the `C`-matrix coefficient of real eigenvalues is 1. If `C1 = false`, the `B` and `C` coefficients are chosen in a balanced fashion.
67
+
68
+
See also [`hess_form`](@ref) and [`schur_form`](@ref)
69
+
"""
70
+
functionmodal_form(sys; C1 =false)
71
+
Ab,T =blockdiagonalize(sys.A)
72
+
# Calling similarity_transform looks like a detour, but this implementation allows modal_form to work with any AbstractStateSpace which implements a custom method for similarity transform
73
+
sysm =similarity_transform(sys, T)
74
+
sysm.A .= Ab # sysm.A should already be Ab after similarity_transform, but Ab has less numerical noise
75
+
if ControlSystems.issiso(sysm)
76
+
# This enforces a convention: the C matrix entry for the first component in each mode is positive. This allows SISO systems on modal form to be interpolated in a meaningful way by interpolating their coefficients.
77
+
# Ref: "New Metrics Between Rational Spectra and their Connection to Optimal Transport" , Bagge Carlson, Chitre
78
+
ci =complex_indices(sysm.A)
79
+
flips =ones(sysm.nx)
80
+
for i in ci
81
+
if sysm.C[1, i] <0
82
+
flips[i] =-1
83
+
flips[i .+1] =-1
84
+
end
85
+
end
86
+
ri =real_indices(sysm.A)
87
+
for i in ri
88
+
c = sysm.C[1, i]
89
+
if C1
90
+
if c !=0
91
+
flips[i] /= c
92
+
end
93
+
else
94
+
b = sysm.B[i, 1]
95
+
flips[i] *=sqrt(abs(b))/(sqrt(abs(c)) +eps(b))
96
+
end
97
+
end
98
+
T2 =diagm(flips)
99
+
sysm =similarity_transform(sysm, T2)
100
+
T = T*T2
101
+
sysm.A .= Ab # Ab unchanged by diagonal T
102
+
end
103
+
sysm, T
104
+
end
105
+
106
+
"""
107
+
sysm, T = schur_form(sys)
108
+
109
+
Bring `sys` to Schur form.
110
+
111
+
The Schur form is characterized by `A` being Schur with the real values of eigenvalues of `A` on the main diagonal. `T` is the similarity transform applied to the system such that
112
+
```julia
113
+
sysm ≈ similarity_transform(sys, T)
114
+
```
115
+
116
+
See also [`modal_form`](@ref) and [`hess_form`](@ref)
117
+
"""
118
+
functionschur_form(sys)
119
+
SF =schur(sys.A)
120
+
A = SF.T
121
+
B = SF.Z'*sys.B
122
+
C = sys.C*SF.Z
123
+
ss(A,B,C,sys.D, sys.timeevol), SF.Z, SF
124
+
end
125
+
126
+
"""
127
+
sysm, T = hess_form(sys)
128
+
129
+
Bring `sys` to Hessenberg form form.
130
+
131
+
The Hessenberg form is characterized by `A` having upper Hessenberg structure. `T` is the similarity transform applied to the system such that
132
+
```julia
133
+
sysm ≈ similarity_transform(sys, T)
134
+
```
135
+
136
+
See also [`modal_form`](@ref) and [`schur_form`](@ref)
Copy file name to clipboardExpand all lines: src/glover_mcfarlane.jl
+2-2Lines changed: 2 additions & 2 deletions
Original file line number
Diff line number
Diff line change
@@ -154,8 +154,8 @@ end
154
154
155
155
For discrete systems, the `info` tuple contains also feedback gains `F, L` and observer gain `Hkf` such that the controller on observer form is given by
156
156
```math
157
-
x^+ = Ax + Bu + H_{kf}*(Cx - y)\\\\
158
-
u = Fx + L*(Cx - y)
157
+
x^+ = Ax + Bu + H_{kf}(Cx - y)\\\\
158
+
u = Fx + L(Cx - y)
159
159
```
160
160
Note, this controller is *not* strictly proper, i.e., it has a non-zero D matrix.
161
161
The controller can be transformed to observer form for the scaled plant (`info.Gs`)
iszero(P.D11) &&iszero(P.D22) ||error("D11 and D22 must be zero for the standard H2 formulation to be used, try calling with γ=1000 to use the H∞ formulation.")
13
14
X2, Y2, F2, L2 =_solvematrixequations2(P)
14
15
Â2 = P.A + P.B2*F2 + L2*P.C2
15
-
K =ss(Â2, -L2, F2, 0)
16
+
# Af2 = P.A + P.B2*F2
17
+
# C1f2 = P.C1 + P.D12*F2
18
+
# Gc = ss(Af2, I(size(Af2, 1)), C1f2, 0)
19
+
K =ss(Â2, -L2, F2, 0, P.timeevol)
16
20
return K, lft(ss(P), K)
17
21
end
18
22
23
+
iscontinuous(P) ||throw(ArgumentError("h2syn with specified γ is only supported for continuous systems."))
0 commit comments