Merge pull request #40 from JuliaControl/stabunstab

add stable/unstable decomposition

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.2"
+version = "0.3.3"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -26,7 +26,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 [compat]
 ChainRulesCore = "1"
 ComponentArrays = "0.9, 0.10, 0.11"
-ControlSystems = "0.11.9"
+ControlSystems = "0.12"
 DescriptorSystems = "1.2"
 Distributions = "0.25"
 GenericSchur = "0.5.2"

src/RobustAndOptimalControl.jl

@@ -24,7 +24,7 @@ using ChainRulesCore
 export show_construction, vec2sys
 include("utils.jl")
 
-export dss, hinfnorm2, h2norm, hankelnorm, nugap, νgap, baltrunc2
+export dss, hinfnorm2, h2norm, hankelnorm, nugap, νgap, baltrunc2, stab_unstab, baltrunc_unstab
 include("descriptor.jl")
 
 export ExtendedStateSpace, system_mapping, performance_mapping, noise_mapping, ssdata_e, partition, ss

src/descriptor.jl

@@ -86,6 +86,52 @@ 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_unstab(sys::LTISystem; residual = false, n = missing, kwargs...)
+
+Balanced truncation for unstable models. An additive decomposition of sys into `sys = sys_stable + sys_unstable` is performed after which `sys_stable` is reduced. The order `n` must not be less than the number of unstable poles.
+
+See `baltrunc2` for other keyword arguments.
+"""
+function baltrunc_unstab(sys::LTISystem; residual=false, n=missing, kwargs...)
+    stab, unstab = DescriptorSystems.gsdec(dss(sys); job="stable", kwargs...)
+    nx_unstab = size(unstab.A, 1)
+    if n isa Integer && n < nx_unstab
+        error("The model contains $(nx_unstab) poles outside the stability region, the reduced-order model must be of at least this order.")
+    end
+    sysr, hs = DescriptorSystems.gbalmr(stab; matchdc=residual, ord=n-nx_unstab, kwargs...)
+    ss(sysr + unstab), hs
+end
+
+"""
+    stab, unstab = stab_unstab(sys; kwargs...)
+
+Decompose `sys` into `sys = stab + unstab` where `stab` contains all stable poles and `unstab` contains unstable poles. See $(@doc(DescriptorSystems.gsdec)) for keyword arguments (argument `job` is set to `"stable"` in this function).
+"""
+function stab_unstab(sys; kwargs...)
+    stab, unstab = DescriptorSystems.gsdec(dss(sys); job="stable", kwargs...)
+    ss(stab), ss(unstab)
+end
+
 ##
 
 function Base.:\(G1::AbstractStateSpace, G2::AbstractStateSpace)

test/test_descriptor.jl

@@ -21,3 +21,29 @@ G2 = ssrand(2,3,4, proper=true)
     end
 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])