From dbce83a3d8204ee0e68459355441a1842bf836f8 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 27 Aug 2025 14:21:12 +0000 Subject: [PATCH 1/5] add scaling capability to model reduction --- src/descriptor.jl | 38 +++++++++++++++++++++++++++++--------- test/test_reduction.jl | 33 ++++++++++++++++++++++++++++++++- 2 files changed, 61 insertions(+), 10 deletions(-) diff --git a/src/descriptor.jl b/src/descriptor.jl index ebd4cd8d..66719c14 100644 --- a/src/descriptor.jl +++ b/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 diff --git a/test/test_reduction.jl b/test/test_reduction.jl index 7d465afe..cc2167b9 100644 --- a/test/test_reduction.jl +++ b/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) \ No newline at end of file +# 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 \ No newline at end of file From 9293132ad6ec82337d2c03a8ded2e1a4de2b384b Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 27 Aug 2025 14:35:27 +0000 Subject: [PATCH 2/5] inverse scaling --- src/descriptor.jl | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/src/descriptor.jl b/src/descriptor.jl index 66719c14..fb29191c 100644 --- a/src/descriptor.jl +++ b/src/descriptor.jl @@ -91,11 +91,12 @@ function baltrunc2(sys::LTISystem; residual=false, n=missing, scaleY=1.0, scaleU # 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) + # Divide by scaling factors to normalize to ~[-1,1] + sys_scaled = ss(A, B / scaleU, C / scaleY, scaleY \ D / scaleU, Ts) sysr, hs = DescriptorSystems.gbalmr(dss(sys_scaled); matchdc=residual, ord=n, kwargs...) - # Inverse scaling after reduction + # Multiply by scaling factors to restore original units Ar, Br, Cr, Dr = ssdata(ss(sysr)) - sys_final = ss(Ar, Br/scaleU, Cr/scaleY, Dr/(scaleY*scaleU), Ts) + sys_final = ss(Ar, Br * scaleU, scaleY * Cr, scaleY * Dr * scaleU, Ts) sys_final, hs end @@ -116,7 +117,8 @@ function baltrunc_coprime(sys, info=nothing; residual=false, n=missing, factoriz # 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) + # Divide by scaling factors to normalize to ~[-1,1] + sys_scaled = ss(A, B / scaleU, scaleY \ C, scaleY \ D / scaleU, Ts) if info !== nothing && hasproperty(info, :NM) @unpack N, M, NM = info else @@ -139,8 +141,8 @@ function baltrunc_coprime(sys, info=nothing; residual=false, n=missing, factoriz Br = BN - BM * (DMi * DN) Dr = (DMi * DN) - # Inverse scaling after reduction - sys_final = ss(Ar, Br/scaleU, Cr/scaleY, Dr/(scaleY*scaleU), Ts) + # Multiply by scaling factors to restore original units + sys_final = ss(Ar, Br * scaleU, scaleY * Cr, scaleY * Dr * scaleU, Ts) sys_final, hs, (; NM, N, M, NMr) end @@ -156,7 +158,8 @@ function baltrunc_unstab(sys::LTISystem, info=nothing; residual=false, n=missing # 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) + # Divide by scaling factors to normalize to ~[-1,1] + sys_scaled = ss(A, B / scaleU, scaleY \ C, scaleY \ D / scaleU, Ts) if info !== nothing && hasproperty(info, :stab) @unpack stab, unstab = info else @@ -167,9 +170,9 @@ function baltrunc_unstab(sys::LTISystem, info=nothing; residual=false, n=missing 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...) - # Inverse scaling after reduction + # Multiply by scaling factors to restore original units Ar, Br, Cr, Dr = ssdata(ss(sysr + unstab)) - sys_final = ss(Ar, Br/scaleU, Cr/scaleY, Dr/(scaleY*scaleU), Ts) + sys_final = ss(Ar, Br * scaleU, scaleY * Cr, scaleY * Dr * scaleU, Ts) sys_final, hs, (; stab, unstab) end From 3b38e7d5ae27c89611964c3b1cc154eed154d95e Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Thu, 28 Aug 2025 04:40:02 +0000 Subject: [PATCH 3/5] timeevol --- src/descriptor.jl | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/src/descriptor.jl b/src/descriptor.jl index fb29191c..48a1af6d 100644 --- a/src/descriptor.jl +++ b/src/descriptor.jl @@ -89,14 +89,13 @@ $(@doc(DescriptorSystems.gbalmr)) """ 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) # Divide by scaling factors to normalize to ~[-1,1] - sys_scaled = ss(A, B / scaleU, C / scaleY, scaleY \ D / scaleU, Ts) + sys_scaled = ss(A, B / scaleU, C / scaleY, scaleY \ D / scaleU, sys.timeevol) sysr, hs = DescriptorSystems.gbalmr(dss(sys_scaled); matchdc=residual, ord=n, kwargs...) # Multiply by scaling factors to restore original units Ar, Br, Cr, Dr = ssdata(ss(sysr)) - sys_final = ss(Ar, Br * scaleU, scaleY * Cr, scaleY * Dr * scaleU, Ts) + sys_final = ss(Ar, Br * scaleU, scaleY * Cr, scaleY * Dr * scaleU, sys.timeevol) sys_final, hs end @@ -115,10 +114,9 @@ $(@doc(DescriptorSystems.gbalmr)) """ 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) # Divide by scaling factors to normalize to ~[-1,1] - sys_scaled = ss(A, B / scaleU, scaleY \ C, scaleY \ D / scaleU, Ts) + sys_scaled = ss(A, B / scaleU, scaleY \ C, scaleY \ D / scaleU, sys.timeevol) if info !== nothing && hasproperty(info, :NM) @unpack N, M, NM = info else @@ -142,7 +140,7 @@ function baltrunc_coprime(sys, info=nothing; residual=false, n=missing, factoriz Dr = (DMi * DN) # Multiply by scaling factors to restore original units - sys_final = ss(Ar, Br * scaleU, scaleY * Cr, scaleY * Dr * scaleU, Ts) + sys_final = ss(Ar, Br * scaleU, scaleY * Cr, scaleY * Dr * scaleU, sys.timeevol) sys_final, hs, (; NM, N, M, NMr) end @@ -156,10 +154,9 @@ See `baltrunc2` for other keyword arguments. """ 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) # Divide by scaling factors to normalize to ~[-1,1] - sys_scaled = ss(A, B / scaleU, scaleY \ C, scaleY \ D / scaleU, Ts) + sys_scaled = ss(A, B / scaleU, scaleY \ C, scaleY \ D / scaleU, sys.timeevol) if info !== nothing && hasproperty(info, :stab) @unpack stab, unstab = info else @@ -172,7 +169,7 @@ function baltrunc_unstab(sys::LTISystem, info=nothing; residual=false, n=missing sysr, hs = DescriptorSystems.gbalmr(stab; matchdc=residual, ord=n-nx_unstab, kwargs...) # Multiply by scaling factors to restore original units Ar, Br, Cr, Dr = ssdata(ss(sysr + unstab)) - sys_final = ss(Ar, Br * scaleU, scaleY * Cr, scaleY * Dr * scaleU, Ts) + sys_final = ss(Ar, Br * scaleU, scaleY * Cr, scaleY * Dr * scaleU, sys.timeevol) sys_final, hs, (; stab, unstab) end From 1c20bf286d86d0f45f638090ae9f95527430ace9 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Thu, 28 Aug 2025 06:33:59 +0000 Subject: [PATCH 4/5] better tests --- test/test_reduction.jl | 39 +++++++++++++++++++++------------------ 1 file changed, 21 insertions(+), 18 deletions(-) diff --git a/test/test_reduction.jl b/test/test_reduction.jl index cc2167b9..9e76e690 100644 --- a/test/test_reduction.jl +++ b/test/test_reduction.jl @@ -303,32 +303,35 @@ controller_reduction_plot(info.Gs,info.Ks, method=:cr) # ncfmargin(P, W1*Ksr) -# === Scaling tests for model reduction === + +# === Scaling invariance 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) + +sys_siso = ssrand(1,1,9,proper=true) # 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 +# baltrunc2: compare frequency response with and without scaling +sysr_noscale, _ = baltrunc2(sys_siso; n=3) +sysr_scale, _ = baltrunc2(sys_siso; n=3, scaleY=scaleY, scaleU=scaleU) +@test sysr_noscale.nx == sysr_scale.nx == 3 +@test hinfnorm2(minreal(sysr_noscale - sysr_scale))[1] < 1e-5 + + +# baltrunc_coprime: compare frequency response with and without scaling +sysr_coprime_noscale, _, _ = baltrunc_coprime(sys_siso; n=3) +sysr_coprime_scale, _, _ = baltrunc_coprime(sys_siso; n=3, scaleY=scaleY, scaleU=scaleU) +@test sysr_coprime_noscale.nx == sysr_coprime_scale.nx == 3 +@test hinfnorm2(minreal(sysr_coprime_noscale - sysr_coprime_scale))[1] < 1e-5 -# 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 +# baltrunc_unstab: compare frequency response with and without scaling 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 \ No newline at end of file +sysr_unstab_noscale, _, _ = baltrunc_unstab(sys_unstab; n=3) +sysr_unstab_scale, _, _ = baltrunc_unstab(sys_unstab; n=3, scaleY=scaleY, scaleU=scaleU) +@test sysr_unstab_noscale.nx == sysr_unstab_scale.nx == 3 +@test hinfnorm2(minreal(sysr_unstab_noscale - sysr_unstab_scale))[1] < 1e-5 From fedb41a7cd4f45abea74c89873a597880f21ad99 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Thu, 28 Aug 2025 06:58:10 +0000 Subject: [PATCH 5/5] up test --- test/test_reduction.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/test/test_reduction.jl b/test/test_reduction.jl index 9e76e690..fd24168f 100644 --- a/test/test_reduction.jl +++ b/test/test_reduction.jl @@ -330,8 +330,7 @@ sysr_coprime_scale, _, _ = baltrunc_coprime(sys_siso; n=3, scaleY=scaleY, scaleU # baltrunc_unstab: compare frequency response with and without scaling -sys_unstab = ss([1.0], [1.0], [1.0], [0.0]) # Unstable pole -sysr_unstab_noscale, _, _ = baltrunc_unstab(sys_unstab; n=3) -sysr_unstab_scale, _, _ = baltrunc_unstab(sys_unstab; n=3, scaleY=scaleY, scaleU=scaleU) +sysr_unstab_noscale, _, _ = baltrunc_unstab(sys_siso; n=3) +sysr_unstab_scale, _, _ = baltrunc_unstab(sys_siso; n=3, scaleY=scaleY, scaleU=scaleU) @test sysr_unstab_noscale.nx == sysr_unstab_scale.nx == 3 @test hinfnorm2(minreal(sysr_unstab_noscale - sysr_unstab_scale))[1] < 1e-5