From 0a540abfa2c9a0f5ddbb0eedc0de9894f4c081b6 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Mon, 21 Feb 2022 14:13:48 +0100 Subject: [PATCH 1/3] add balanced truncation on coprime factors --- Project.toml | 2 +- src/RobustAndOptimalControl.jl | 2 +- src/descriptor.jl | 48 ++++++++++++++++++++++------------ test/test_descriptor.jl | 15 +++++++---- 4 files changed, 43 insertions(+), 24 deletions(-) diff --git a/Project.toml b/Project.toml index 5f9c1fd9..320b2546 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "RobustAndOptimalControl" uuid = "21fd56a4-db03-40ee-82ee-a87907bee541" authors = ["Fredrik Bagge Carlson", "Marcus Greiff"] -version = "0.3.3" +version = "0.3.4" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" diff --git a/src/RobustAndOptimalControl.jl b/src/RobustAndOptimalControl.jl index cd620ad3..d15b9d56 100644 --- a/src/RobustAndOptimalControl.jl +++ b/src/RobustAndOptimalControl.jl @@ -24,7 +24,7 @@ using ChainRulesCore export show_construction, vec2sys include("utils.jl") -export dss, hinfnorm2, h2norm, hankelnorm, nugap, νgap, baltrunc2, stab_unstab, baltrunc_unstab +export dss, hinfnorm2, h2norm, hankelnorm, nugap, νgap, baltrunc2, stab_unstab, baltrunc_unstab, baltrunc_coprime include("descriptor.jl") export ExtendedStateSpace, system_mapping, performance_mapping, noise_mapping, ssdata_e, partition, ss diff --git a/src/descriptor.jl b/src/descriptor.jl index ec3e20a0..8efe7f1b 100644 --- a/src/descriptor.jl +++ b/src/descriptor.jl @@ -86,23 +86,37 @@ function baltrunc2(sys::LTISystem; residual=false, n=missing, kwargs...) ss(sysr), hs end -# function coprime_baltrunc(sys; residual=false, n=missing, kwargs...) -# N,M = DescriptorSystems.glcf(dss(sys), mindeg=true, mininf=true, fast=false) -# nu = size(N.B, 2) -# # A,E,B,C,D = DescriptorSystems.dssdata(N) -# # NM = DescriptorSystems.dss(A,E,[B M.B],C,[D M.D]) -# # Nr1, hs = DescriptorSystems.gbalmr(N; matchdc=residual, ord=n-size(M.A, 1), kwargs...) -# NM = [N M] -# sysr, hs = DescriptorSystems.gbalmr(NM; matchdc=residual, ord=n, kwargs...) - -# A,E,B,C,D = DescriptorSystems.dssdata(sysr) - -# Nr = DescriptorSystems.dss(A,E,B[:, 1:nu],C,D[:, 1:nu]) -# Mr = DescriptorSystems.dss(A,E,B[:, nu+1:end],C,D[:, nu+1:end]) -# sysr = Mr \ Nr -# ss(sysr), hs -# # Nr, Mr -# end +""" + baltrunc_coprime(sys; residual = false, n = missing, factorization::F = DescriptorSystems.glcf, kwargs...) + +Compute a balanced truncation of the left coprime factorization of `sys`. +See [`baltrunc2`](@ref) for additional keyword-argument help. + +# Arguments: +- `factorization`: The function to perform the coprime factorization. A normalized factorization may be used by passing `RobustAndOptimalControl.DescriptorSystems.gnlcf`. +- `kwargs`: Are passed to `DescriptorSystems.gbalmr` +""" +function baltrunc_coprime(sys; residual=false, n=missing, factorization::F = DescriptorSystems.glcf, kwargs...) where F + N,M = factorization(dss(sys)) + nu = size(N.B, 2) + A,E,B,C,D = DescriptorSystems.dssdata(N) + NM = DescriptorSystems.dss(A,E,[B M.B],C,[D M.D]) + sysr, hs = DescriptorSystems.gbalmr(NM; matchdc=residual, ord=n, kwargs...) + + A,E,B,C,D = DescriptorSystems.dssdata(DescriptorSystems.dss2ss(sysr)[1]) + + BN = B[:, 1:nu] + DN = D[:, 1:nu] + BM = B[:, nu+1:end] + DM = D[:, nu+1:end] + + Ar = A - BM * (DM \ C) + Cr = (DM \ C) + Br = BN - BM * (DM \ DN) + Dr = (DM \ DN) + + ss(Ar,Br,Cr,Dr,sys.timeevol), hs +end """ diff --git a/test/test_descriptor.jl b/test/test_descriptor.jl index 3b77a19d..f746def5 100644 --- a/test/test_descriptor.jl +++ b/test/test_descriptor.jl @@ -23,11 +23,15 @@ end >= 94 ## -# sys = ssrand(2,3,4, stable=false) -# # N,M = coprime_baltrunc(sys, n=3) -# sysr, _ = coprime_baltrunc(sys, n=3) +sys = ssrand(2,3,40, stable=true) +sysus = ssrand(2,3,2, stable=true) +sysus.A .*= -1 +sys = sys + sysus -# @test sysr.nx == 3 +sysr, hs = baltrunc_coprime(sys, n=20, factorization = RobustAndOptimalControl.DescriptorSystems.glcf) + +@test sysr.nx == 20 +@test linfnorm(sysr - sys)[1] < 3e-3 # bodeplot([sys, sysr]) @@ -46,4 +50,5 @@ sys = sys + sysus sysr, hs = baltrunc_unstab(sys, n=20) @test sysr.nx == 20 @test linfnorm(sysr - sys)[1] < 1e-3 -# bodeplot([sys, sysr]) \ No newline at end of file +# bodeplot([sys, sysr]) + From 6e6e89a031b6a53902de4c117699d635f1188f90 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Mon, 21 Feb 2022 14:16:03 +0100 Subject: [PATCH 2/3] add test for number of unstable poles --- test/test_descriptor.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/test/test_descriptor.jl b/test/test_descriptor.jl index f746def5..7bac2bf5 100644 --- a/test/test_descriptor.jl +++ b/test/test_descriptor.jl @@ -28,11 +28,14 @@ sysus = ssrand(2,3,2, stable=true) sysus.A .*= -1 sys = sys + sysus -sysr, hs = baltrunc_coprime(sys, n=20, factorization = RobustAndOptimalControl.DescriptorSystems.glcf) +sysr, hs = RobustAndOptimalControl.baltrunc_coprime(sys, n=20, factorization = RobustAndOptimalControl.DescriptorSystems.glcf) -@test sysr.nx == 20 +@test sysr.nx <= 20 @test linfnorm(sysr - sys)[1] < 3e-3 +e = poles(sysr) +@test count(e->real(e)>0, e) == 2 # test that the two unstable poles were preserved + # bodeplot([sys, sysr]) ## stab_unstab @@ -48,7 +51,7 @@ sysus = ssrand(2,3,2, stable=true) sysus.A .*= -1 sys = sys + sysus sysr, hs = baltrunc_unstab(sys, n=20) -@test sysr.nx == 20 +@test sysr.nx <= 20 @test linfnorm(sysr - sys)[1] < 1e-3 # bodeplot([sys, sysr]) From 4f4667c4f5fc9c87ba61282f4564d5e964e47371 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Mon, 21 Feb 2022 15:12:11 +0100 Subject: [PATCH 3/3] prefactor matrix --- src/descriptor.jl | 10 +++---- test/test_descriptor.jl | 66 ++++++++++++++++++++++------------------- 2 files changed, 40 insertions(+), 36 deletions(-) diff --git a/src/descriptor.jl b/src/descriptor.jl index 8efe7f1b..0586e2b7 100644 --- a/src/descriptor.jl +++ b/src/descriptor.jl @@ -108,12 +108,12 @@ function baltrunc_coprime(sys; residual=false, n=missing, factorization::F = Des BN = B[:, 1:nu] DN = D[:, 1:nu] BM = B[:, nu+1:end] - DM = D[:, nu+1:end] + DMi = pinv(D[:, nu+1:end]) - Ar = A - BM * (DM \ C) - Cr = (DM \ C) - Br = BN - BM * (DM \ DN) - Dr = (DM \ DN) + Ar = A - BM * (DMi * C) + Cr = (DMi * C) + Br = BN - BM * (DMi * DN) + Dr = (DMi * DN) ss(Ar,Br,Cr,Dr,sys.timeevol), hs end diff --git a/test/test_descriptor.jl b/test/test_descriptor.jl index 7bac2bf5..d1d3d199 100644 --- a/test/test_descriptor.jl +++ b/test/test_descriptor.jl @@ -23,35 +23,39 @@ end >= 94 ## -sys = ssrand(2,3,40, stable=true) -sysus = ssrand(2,3,2, stable=true) -sysus.A .*= -1 -sys = sys + sysus - -sysr, hs = RobustAndOptimalControl.baltrunc_coprime(sys, n=20, factorization = RobustAndOptimalControl.DescriptorSystems.glcf) - -@test sysr.nx <= 20 -@test linfnorm(sysr - sys)[1] < 3e-3 - -e = poles(sysr) -@test count(e->real(e)>0, e) == 2 # test that the two unstable poles were preserved - -# bodeplot([sys, sysr]) - -## stab_unstab -sys = ssrand(2,3,40, stable=false) -stab, unstab = stab_unstab(sys) -@test all(real(poles(stab)) .< 0) -@test all(real(poles(unstab)) .>= 0) -@test linfnorm(stab + unstab - sys)[1] < 1e-8 - -## baltrunc_unstab -sys = ssrand(2,3,40, stable=true) -sysus = ssrand(2,3,2, stable=true) -sysus.A .*= -1 -sys = sys + sysus -sysr, hs = baltrunc_unstab(sys, n=20) -@test sysr.nx <= 20 -@test linfnorm(sysr - sys)[1] < 1e-3 -# bodeplot([sys, sysr]) +@testset "unstable baltrunc" begin + @info "Testing unstable baltrunc" + for proper = [true, false] + sys = ssrand(2,3,40; stable=true, proper) + sysus = ssrand(2,3,2; stable=true, proper) + sysus.A .*= -1 + sys = sys + sysus + + sysr, hs = RobustAndOptimalControl.baltrunc_coprime(sys, n=20, factorization = RobustAndOptimalControl.DescriptorSystems.glcf) + + @test sysr.nx <= 20 + @test linfnorm(sysr - sys)[1] < 9e-3 + + e = poles(sysr) + @test count(e->real(e)>0, e) == 2 # test that the two unstable poles were preserved + # bodeplot([sys, sysr]) + end + ## stab_unstab + sys = ssrand(2,3,40, stable=false) + stab, unstab = stab_unstab(sys) + @test all(real(poles(stab)) .< 0) + @test all(real(poles(unstab)) .>= 0) + @test linfnorm(stab + unstab - sys)[1] < 1e-8 + + ## baltrunc_unstab + sys = ssrand(2,3,40, stable=true) + sysus = ssrand(2,3,2, stable=true) + sysus.A .*= -1 + sys = sys + sysus + sysr, hs = baltrunc_unstab(sys, n=20) + @test sysr.nx <= 20 + @test linfnorm(sysr - sys)[1] < 1e-3 + # bodeplot([sys, sysr]) + +end \ No newline at end of file