Skip to content

Commit e4eefce

Browse files
baggepinnenbaggepinnen
committed
add some canonical forms
1 parent 9909007 commit e4eefce

4 files changed

Lines changed: 265 additions & 0 deletions

File tree

‎src/RobustAndOptimalControl.jl‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,9 @@ include("descriptor.jl")
3030
export ExtendedStateSpace, system_mapping, performance_mapping, noise_mapping, ssdata_e, partition, ss
3131
include("ExtendedStateSpace.jl")
3232

33+
export modal_form, schur_form, hess_form
34+
include("canonical.jl")
35+
3336
export δ, δr, δc, δss, nominal, UncertainSS, uss, blocksort
3437
include("uncertainty_interface.jl")
3538

‎src/canonical.jl‎

Lines changed: 146 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,146 @@
1+
blockdiagonalize(A::AbstractMatrix) = cdf2rdf(eigen(A, sortby=eigsortby))
2+
3+
eigsortby(λ::Real) = λ
4+
eigsortby(λ::Complex) = (abs(imag(λ)),real(λ))
5+
6+
function complex_indices(A::Matrix) # assumes A on block diagonal form
7+
findall(diag(A, -1) .!= 0)
8+
end
9+
10+
function real_indices(A::Matrix) # assumes A on block diagonal form
11+
size(A,1) == 1 && return [1]
12+
setdiff(findall(diag(A, -1) .== 0), complex_indices(A).+1)
13+
end
14+
15+
function complex_indices(D::AbstractVector)
16+
complex_eigs = imag.(D) .!= 0
17+
findall(complex_eigs)
18+
end
19+
20+
function cdf2rdf(E::Eigen)
21+
# 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) in zip(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+
function modal_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+
function schur_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)
137+
"""
138+
function hess_form(sys)
139+
F = hessenberg(sys.A)
140+
Q = Matrix(F.Q)
141+
A = F.H
142+
B = Q'sys.B
143+
C = sys.C*Q
144+
D = sys.D
145+
ss(A,B,C,D, sys.timeevol), Q, F
146+
end

‎test/runtests.jl‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -101,4 +101,9 @@ using Test
101101
include("test_high_precision_hinf.jl")
102102
end
103103

104+
@testset "canonical_forms" begin
105+
@info "Testing canonical_forms"
106+
include("test_canonical_forms.jl")
107+
end
108+
104109
end

‎test/test_canonical_forms.jl‎

Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
using RobustAndOptimalControl, ControlSystems
2+
import RobustAndOptimalControl: cdf2rdf, blockdiagonalize
3+
4+
function isblockdiagonal(A)
5+
complex_inds = findall(diag(A, -1) .!= 0)
6+
diag(A, -1) ≈ -diag(A, 1) || (return false)
7+
A = A - diagm(diag(A)) # remove main diagonal
8+
A = A - diagm(1=>diag(A, 1))
9+
A = A - diagm(-1=>diag(A, -1))
10+
all(iszero, A)
11+
end
12+
13+
14+
@testset "modal form" begin
15+
@info "Testing modal form"
16+
17+
X = randn(5,5) # odd number to enforce at least one real eigval
18+
E = eigen(X)
19+
D,V = E
20+
Db, Vb = cdf2rdf(E)
21+
@test isblockdiagonal(Db)
22+
@test Vb*Db ≈ X*Vb
23+
@test sort(real(D)) ≈ sort(diag(Db)) # real values on diagonal
24+
ivals = [diag(Db, -1); diag(Db, 1)] # imag values on 1/-1 diagonals
25+
@test all(v ∈ ivals for v in imag(D))
26+
27+
Xb, T = blockdiagonalize(X)
28+
@test T\X*T ≈ Xb
29+
@test isblockdiagonal(Xb)
30+
31+
sys = ssrand(1,1,5) # odd number to enforce at least one real eigval
32+
sysm, T = modal_form(sys)
33+
@test tf(sys) ≈ tf(sysm)
34+
@test sysm ≈ similarity_transform(sys, T)
35+
36+
complex_inds = findall(diag(sysm.A, -1) .!= 0)
37+
@test all(sysm.C[i] > 0 for i ∈ complex_inds) # test coefficient convention
38+
39+
40+
# test balanced property
41+
sys = ssrand(1,1,1, proper=true)
42+
sysm, T = modal_form(sys)
43+
@test tf(sys) ≈ tf(sysm)
44+
@test sysm ≈ similarity_transform(sys, T)
45+
@test abs(sysm.C[1]) ≈ abs(sysm.B[1])
46+
47+
# test C1 property
48+
sys = ssrand(1,1,1, proper=true)
49+
sysm, T = modal_form(sys, C1=true)
50+
@test tf(sys) ≈ tf(sysm)
51+
@test sysm ≈ similarity_transform(sys, T)
52+
@test sysm.C[1] ≈ 1
53+
54+
55+
sys = ssrand(2,3,5) # test that it works for MIMO
56+
sysm, T = modal_form(sys)
57+
@test tf(sys) ≈ tf(sysm)
58+
@test isblockdiagonal(sysm.A)
59+
60+
@test sysm ≈ similarity_transform(sys, T)
61+
62+
63+
# test with repeated eigenvalues
64+
X = ss(tf(1,[1,1,1])^2).A
65+
E = eigen(X)
66+
Db, Vb = cdf2rdf(E)
67+
@test isblockdiagonal(Db)
68+
@test Vb*Db ≈ X*Vb
69+
@test sort(real(E.values)) ≈ sort(diag(Db)) # real values on diagonal
70+
ivals = [diag(Db, -1); diag(Db, 1)] # imag values on 1/-1 diagonals
71+
@test all(v ∈ ivals for v in imag(E.values))
72+
end
73+
74+
75+
76+
77+
@testset "schur form" begin
78+
@info "Testing schur form"
79+
80+
sys = ssrand(1,1,5) # odd number to enforce at least one real eigval
81+
sysm, T = schur_form(sys)
82+
@test tf(sys) ≈ tf(sysm)
83+
84+
@test sysm ≈ similarity_transform(sys, T)
85+
86+
87+
sys = ssrand(2,3,5) # test that it works for MIMO
88+
sysm, _ = schur_form(sys)
89+
@test tf(sys) ≈ tf(sysm)
90+
@test all(d->all(iszero, sysm.A[diagind(sysm.A, d)]), -5:-2)
91+
end
92+
93+
94+
95+
@testset "hess form" begin
96+
@info "Testing hess form"
97+
98+
sys = ssrand(1,1,5) # odd number to enforce at least one real eigval
99+
sysm, T = hess_form(sys)
100+
@test tf(sys) ≈ tf(sysm)
101+
102+
@test sysm ≈ similarity_transform(sys, T)
103+
104+
105+
sys = ssrand(2,3,5) # test that it works for MIMO
106+
sysm, _ = hess_form(sys)
107+
@test tf(sys) ≈ tf(sysm)
108+
@test all(d->all(iszero, sysm.A[diagind(sysm.A, d)]), -5:-2)
109+
end
110+
111+

0 commit comments

Comments
 (0)