add scaling capability to model reduction

Author: baggepinnen

src/descriptor.jl

@@ -87,9 +87,16 @@ Compute the a balanced truncation of order `n` and the hankel singular values
 For keyword arguments, see the docstring of `DescriptorSystems.gbalmr`, reproduced below
 $(@doc(DescriptorSystems.gbalmr))
 """
-function baltrunc2(sys::LTISystem; residual=false, n=missing, kwargs...)
-    sysr, hs = DescriptorSystems.gbalmr(dss(sys); matchdc=residual, ord=n, kwargs...)
-    ss(sysr), hs
+function baltrunc2(sys::LTISystem; residual=false, n=missing, scaleY=1.0, scaleU=1.0, kwargs...)
+    # Apply scaling if needed
+    Ts = isdiscrete(sys) ? sys.Ts : 0
+    A, B, C, D = ssdata(sys)
+    sys_scaled = ss(A, B*scaleU, scaleY*C, scaleY*D*scaleU, Ts)
+    sysr, hs = DescriptorSystems.gbalmr(dss(sys_scaled); matchdc=residual, ord=n, kwargs...)
+    # Inverse scaling after reduction
+    Ar, Br, Cr, Dr = ssdata(ss(sysr))
+    sys_final = ss(Ar, Br/scaleU, Cr/scaleY, Dr/(scaleY*scaleU), Ts)
+    sys_final, hs
 end
 
 """
@@ -105,11 +112,15 @@ Coprime-factor reduction performs a coprime factorization of the model into \$P(
 - `kwargs`: Are passed to `DescriptorSystems.gbalmr`, the docstring of which is reproduced below:
 $(@doc(DescriptorSystems.gbalmr))
 """
-function baltrunc_coprime(sys, info=nothing; residual=false, n=missing, factorization::F = DescriptorSystems.gnlcf, kwargs...) where F
+function baltrunc_coprime(sys, info=nothing; residual=false, n=missing, factorization::F = DescriptorSystems.gnlcf, scaleY=1.0, scaleU=1.0, kwargs...) where F
+    # Apply scaling if needed
+    Ts = isdiscrete(sys) ? sys.Ts : 0
+    A, B, C, D = ssdata(sys)
+    sys_scaled = ss(A, B*scaleU, scaleY*C, scaleY*D*scaleU, Ts)
     if info !== nothing && hasproperty(info, :NM)
         @unpack N, M, NM = info
     else
-        N,M = factorization(dss(sys))
+        N,M = factorization(dss(sys_scaled))
         A,E,B,C,D = DescriptorSystems.dssdata(N)
         NM = DescriptorSystems.dss(A,E,[B M.B],C,[D M.D])
     end
@@ -128,7 +139,9 @@ function baltrunc_coprime(sys, info=nothing; residual=false, n=missing, factoriz
     Br = BN  - BM * (DMi * DN)
     Dr = (DMi * DN)
 
-    ss(Ar,Br,Cr,Dr,sys.timeevol), hs, (; NM, N, M, NMr)
+    # Inverse scaling after reduction
+    sys_final = ss(Ar, Br/scaleU, Cr/scaleY, Dr/(scaleY*scaleU), Ts)
+    sys_final, hs, (; NM, N, M, NMr)
 end
 
 
@@ -139,18 +152,25 @@ Balanced truncation for unstable models. An additive decomposition of sys into `
 
 See `baltrunc2` for other keyword arguments.
 """
-function baltrunc_unstab(sys::LTISystem, info=nothing; residual=false, n=missing, kwargs...)
+function baltrunc_unstab(sys::LTISystem, info=nothing; residual=false, n=missing, scaleY=1.0, scaleU=1.0, kwargs...)
+    # Apply scaling if needed
+    Ts = isdiscrete(sys) ? sys.Ts : 0
+    A, B, C, D = ssdata(sys)
+    sys_scaled = ss(A, B*scaleU, scaleY*C, scaleY*D*scaleU, Ts)
     if info !== nothing && hasproperty(info, :stab)
         @unpack stab, unstab = info
     else
-        stab, unstab = DescriptorSystems.gsdec(dss(sys); job="stable", kwargs...)
+        stab, unstab = DescriptorSystems.gsdec(dss(sys_scaled); job="stable", kwargs...)
     end
     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, (; stab, unstab)
+    # Inverse scaling after reduction
+    Ar, Br, Cr, Dr = ssdata(ss(sysr + unstab))
+    sys_final = ss(Ar, Br/scaleU, Cr/scaleY, Dr/(scaleY*scaleU), Ts)
+    sys_final, hs, (; stab, unstab)
 end
 
 

test/test_reduction.jl

@@ -300,4 +300,35 @@ Ksr, hs, infor = baltrunc_coprime(info.Ks; n)
 controller_reduction_plot(info.Gs,info.Ks)
 controller_reduction_plot(info.Gs,info.Ks, method=:cr)
 
-# ncfmargin(P, W1*Ksr)
+# ncfmargin(P, W1*Ksr)
+
+
+# === Scaling tests for model reduction ===
+using RobustAndOptimalControl: baltrunc2, baltrunc_coprime, baltrunc_unstab
+
+# Simple SISO system
+As = [-1.0]
+Bs = [2.0]
+Cs = [3.0]
+Ds = [0.5]
+sys_siso = ss(As, Bs, Cs, Ds)
+
+# Non-unit scaling factors
+scaleY = 5.0
+scaleU = 0.2
+
+# baltrunc2 scaling test
+sysr, _ = baltrunc2(sys_siso; n=1, scaleY, scaleU)
+@test sysr.nx == 1
+@test abs(dcgain(sys_siso) - dcgain(sysr)) < 1e-8
+
+# baltrunc_coprime scaling test
+sysr_coprime, _, _ = baltrunc_coprime(sys_siso; n=1, scaleY, scaleU)
+@test sysr_coprime.nx == 1
+@test abs(dcgain(sys_siso) - dcgain(sysr_coprime)) < 1e-8
+
+# baltrunc_unstab scaling test
+sys_unstab = ss([1.0], [1.0], [1.0], [0.0]) # Unstable pole
+sysr_unstab, _, _ = baltrunc_unstab(sys_unstab; n=1, scaleY, scaleU)
+@test sysr_unstab.nx == 1
+@test abs(dcgain(sys_unstab) - dcgain(sysr_unstab)) < 1e-8