add balanced truncation on coprime factors

Author: baggepinnen

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]
+    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
 
 
 """

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])
+# bodeplot([sys, sysr])
+