add balanced truncation on coprime factors

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"

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

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]
+    DMi = pinv(D[:, nu+1:end])
+
+    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
 
 
 """

test/test_descriptor.jl

@@ -23,27 +23,39 @@ end >= 94
 
 
 ##
-# sys = ssrand(2,3,4, stable=false)
-# # N,M = coprime_baltrunc(sys, n=3)
-# sysr, _ = coprime_baltrunc(sys, n=3)
-
-# @test sysr.nx == 3
-
-# 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