add utility to convert uncertain model to multimodel with many outputs

Author: baggepinnen

Project.toml

@@ -18,6 +18,7 @@ Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
@@ -35,6 +36,7 @@ MonteCarloMeasurements = "1.0"
 Optim = "1.5"
 Plots = "1"
 RecipesBase = "1"
+SparseArrays = "1"
 Statistics = "1"
 UnPack = "1.0"
 julia = "1.7"

src/RobustAndOptimalControl.jl

@@ -37,7 +37,7 @@ include("ExtendedStateSpace.jl")
 export modal_form, schur_form, hess_form, frequency_separation
 include("canonical.jl")
 
-export δ, δr, δc, δss, nominal, UncertainSS, uss, blocksort, sys_from_particles, ss2particles
+export δ, δr, δc, δss, nominal, UncertainSS, uss, blocksort, sys_from_particles, mo_sys_from_particles, ss2particles
 include("uncertainty_interface.jl")
 
 export makeweight, neglected_delay, gain_and_delay_uncertainty, neglected_lag, fit_complex_perturbations

src/named_systems2.jl

@@ -277,9 +277,9 @@ function Base.show(io::IO, G::NamedStateSpace)
     end
     print(io, "Named")
     show(io, G.sys)
-    print(io, "\nWith state  names: "); println(io, join(G.x, ' '))
-    print(io, "     input  names: "); println(io, join(G.u, ' '))
-    print(io, "     output names: "); println(io, join(G.y, ' '))
+    length(G.x) < 50 && (print(io, "\nWith state  names: "); println(io, join(G.x, ' ')))
+    length(G.u) < 50 && (print(io, "     input  names: "); println(io, join(G.u, ' ')))
+    length(G.y) < 50 && (print(io, "     output names: "); println(io, join(G.y, ' ')))
 end
 
 function Base.:-(s1::NamedStateSpace{T,S}) where {T <: CS.TimeEvolution, S}

src/uncertainty_interface.jl

@@ -1,3 +1,4 @@
+using SparseArrays
 abstract type UncertainElement{C, F} end
 
 struct δ{C, F<:Real} <: UncertainElement{C, F}
@@ -477,11 +478,11 @@ using MonteCarloMeasurements: vecindex
 
 Return the `i`th system from a system `P` with `Particles` coefficients.
 
-If called without an index, return a vector of systems, one for each possibly `i`.
+If called without an index, return a vector of systems, one for each possible `i`.
 
 This function is used to convert from an uncertain representation using `Particles` to a "multi-model" representation using multiple `StateSpace` models.
 
-See also [`ss2particles`](@ref) and `MonteCarloMeasurements.nominal`.
+See also [`ss2particles`](@ref), [`mo_sys_from_particles`](@ref) and `MonteCarloMeasurements.nominal`.
 """
 function sys_from_particles(P, i)
     A,B,C,D = ssdata(P)
@@ -502,6 +503,36 @@ function sys_from_particles(P)
     end
 end
 
+"""
+    mo_sys_from_particles(P::AbstractStateSpace; sparse = nparticles(P.A) > 500)
+
+Converts a state-space model with `Particles` coefficients to a single state-space model with the same number of inputs as `P`, but with `nparticles(P.A)` times the output and state dimensions.
+
+If `sparse` is true, the resulting model will be a `HeteroStateSpace` with sparse `A`, `C`, and `D` matrices.
+"""
+function mo_sys_from_particles(P::AbstractStateSpace; sparse = nparticles(P.A) > 500)
+    systems = sys_from_particles(P)
+    A   = blockdiag(s.A for s in systems)
+    B   = reduce(vcat, s.B for s in systems)
+    C   = blockdiag(s.C for s in systems)
+    D   = reduce(vcat, s.D for s in systems)
+    if sparse
+        if count(!iszero, D) / length(D) < 0.1
+            D = SparseArrays.sparse(D)
+        end
+        Pmo = HeteroStateSpace(SparseArrays.sparse(A), B, SparseArrays.sparse(C), D, ControlSystemsBase.timeevol(P))
+    else
+        Pmo = ss(A, B, C, D, ControlSystemsBase.timeevol(P))
+    end
+    suffix_i(x, i) = Symbol.(string.(x) .* string("_", i))
+    P isa NamedStateSpace || (return Pmo)
+    named_ss(Pmo, P.name;
+        x = reduce(vcat, suffix_i(P.x, i) for i in eachindex(systems)),
+        u = P.u,
+        y = reduce(vcat, suffix_i(P.y, i) for i in eachindex(systems)),
+    )
+end
+
 
 function ControlSystemsBase.lsim(sys::DelayLtiSystem{T,S}, u, t::AbstractArray{<:Real}, args...; x0=fill(zero(T), nstates(sys)), kwargs...) where {T, S<:AbstractParticles}
     syss = sys_from_particles(sys)

test/test_uncertainty.jl

@@ -1,4 +1,4 @@
-using RobustAndOptimalControl, ControlSystemsBase, MonteCarloMeasurements
+using RobustAndOptimalControl, ControlSystemsBase, MonteCarloMeasurements, Test
 
 d = δr()
 @test d.val == 0
@@ -348,4 +348,43 @@ nyquistplot(P*C, w[1:10:end], points=true, xlims=(-3.5, 2.5), ylims=(-5, 1.5), M
 
 P = ss(-1, 1, 1, 0)
 Pd = P * delay(0.1 .. 0.3)
-@test Pd.P.P ≈ (tf(P) * delay(0.1 .. 0.3)).P.P
+@test Pd.P.P ≈ (tf(P) * delay(0.1 .. 0.3)).P.P
+
+
+## Convert a named particle system to a single system with multiple outputs
+# Tests that output and state names are set in a reasonable way
+N = 20
+A = randn(2,2) .+ 0.1Particles(N)
+B = randn(2,3) .+ 0.1Particles(N)
+C = randn(2,2) .+ 0.1Particles(N)
+D = 0
+P = named_ss(ss(A,B,C,D), "P")
+
+MOP = RobustAndOptimalControl.mo_sys_from_particles(P);
+
+@test MOP.ny == P.ny*N
+@test MOP.nx == P.nx*N
+@test MOP.nu == P.nu
+
+@test MOP.x[1] == :Px1_1
+@test MOP.x[2] == :Px2_1
+@test MOP.x[3] == :Px1_2
+@test MOP.u == P.u
+
+# res = step(MOP, 0:0.1:1, method=:tustin)
+
+N = 2000
+A = randn(2,2) .+ 0.1Particles(N)
+B = randn(2,3) .+ 0.1Particles(N)
+C = randn(2,2) .+ 0.1Particles(N)
+D = 0
+P = named_ss(ss(A,B,C,D), "P")
+
+MOP = RobustAndOptimalControl.mo_sys_from_particles(P);
+
+using SparseArrays
+@test MOP.ny == P.ny*N
+@test MOP.nx == P.nx*N
+@test MOP.nu == P.nu
+
+@test MOP.A isa SparseMatrixCSC