Count linearization channels per element of array-valued signals

src/ode_system.jl

@@ -58,6 +58,60 @@ end
 
 symstr(x) = Symbol(x isa AnalysisPoint ? x.name : string(x))
 
+# Number of scalar channels of the signal `x`, a variable or an analysis point (an
+# `AnalysisPoint` or its `Symbol` name) of `sys`. The width of an analysis point is the width
+# of the variable that drives it.
+function signal_width(sys, x)
+    if x isa Union{Symbol, AnalysisPoint}
+        ap = only(ModelingToolkit.canonicalize_ap(sys, x))
+        ap isa Symbol && (ap = AnalysisPoint(ap))
+        _, (x,) = ModelingToolkit.apply_transformation(ModelingToolkit.GetInput(ap), sys)
+    end
+    x = Symbolics.unwrap(x)
+    return Symbolics.symbolic_type(x) == Symbolics.ArraySymbolic() ? length(x) : 1
+end
+
+# One name per scalar channel. A signal named `s` of width `n > 1` is expanded to the channel
+# names `s1, ..., sn`, following `RobustAndOptimalControl.expand_symbol`.
+function channel_names(names, widths)
+    out = Symbol[]
+    sources = Dict{Symbol, Vector{Symbol}}()
+    for (name, w) in zip(names, widths)
+        for c in RobustAndOptimalControl.expand_symbol(name, w)
+            push!(out, c)
+            push!(get!(sources, c, Symbol[]), name)
+        end
+    end
+    if !allunique(out)
+        clashes = join(("`$c` (from $(join(("`$n`" for n in sources[c]), " and ")))" for c in unique(out) if length(sources[c]) > 1), ", ")
+        throw(ArgumentError("The channel names $clashes are not unique. An array-valued signal `s` with `n` elements has the channels `s1, ..., sn`; rename the signals so that these names do not coincide with the names of other signals."))
+    end
+    return out
+end
+
+"""
+    signal_channel_names(sys, signals)
+
+Return the names of the scalar channels of `signals`, in the order in which [`named_ss`](@ref)
+and the `get_named_*` functions name the inputs and outputs of the linearized system. Each
+entry of `signals` is a variable or an analysis point (an `AnalysisPoint` or its name) of
+`sys`. A signal `s` with one element has the channel `s`, and an array-valued signal `s` with
+`n` elements has the channels `s1, ..., sn`. An `ArgumentError` is thrown if two channel
+names coincide.
+"""
+signal_channel_names(sys, signals) =
+    channel_names(symstr.(signals), [signal_width(sys, s) for s in signals])
+
+# Operating-point values of the inputs of `ssys`, one entry per scalar channel
+function input_operating_point(ssys, xpt)
+    u0 = Any[]
+    for i in ModelingToolkit.inputs(ssys)
+        v = xpt.p[ModelingToolkit.parameter_index(ssys, i)]
+        v isa AbstractArray ? append!(u0, vec(v)) : push!(u0, v)
+    end
+    return identity.(u0)
+end
+
 """
     RobustAndOptimalControl.named_ss(sys::ModelingToolkit.AbstractSystem, inputs, outputs; descriptor=true, simple_infeigs=true, balance=false, kwargs...)
 
@@ -96,8 +150,6 @@ function RobustAndOptimalControl.named_ss(
             inp
         end
     end
-    nu = length(inputs)
-
     outputs = map(outputs) do out
         if out isa System
             @variables u(t)
@@ -110,8 +162,12 @@ function RobustAndOptimalControl.named_ss(
             out
         end
     end
+    # The linearization has one channel per element of an array-valued input or output, so
+    # the names and the number of inputs are counted per channel
+    unames = signal_channel_names(sys, inputs)
+    ynames = signal_channel_names(sys, outputs)
+    nu = length(unames)
     matrices, ssys, xpt = ModelingToolkit.linearize(sys, inputs, outputs; DEFAULT_LINEARIZE_KWARGS..., kwargs...)
-    unames = symstr.(inputs)
     if nu > 0 && size(matrices.B, 2) == 2nu
         # This indicates that input derivatives are present
         # Derivative column nu + i is the derivative of input i (see the MTK docs for
@@ -124,18 +180,15 @@ function RobustAndOptimalControl.named_ss(
     else
         lsys = ss(matrices...)
     end
-    pind = [ModelingToolkit.parameter_index(ssys, i) for i in ModelingToolkit.inputs(ssys)]
-    x0 = xpt.x
-    u0 = [xpt.p[pi] for pi in pind] 
-    xu = (; x = x0, u = u0)
+    xu = (; x = xpt.x, u = input_operating_point(ssys, xpt))
     extra = Dict(:operating_point => xu)
     # If simple_infeigs=false, the system might have been reduced and the state names might not match the original system.
     x_names = get_x_names(lsys, ssys; descriptor, simple_infeigs, balance)
     nsys = named_ss(
         lsys;
         x = x_names,
         u = unames,
-        y = symstr.(outputs),
+        y = ynames,
         name = string(Base.nameof(sys)),
         extra,
     )
@@ -157,8 +210,9 @@ function RobustAndOptimalControl.named_ss(
     # the input variables are recovered from the simplified system instead, which
     # is what the operating-point extraction below already relies on.
     inputs = ModelingToolkit.inputs(ssys)
-    nu = length(inputs)
-    unames = symstr.(inputs)
+    unames = signal_channel_names(ssys, inputs)
+    ynames = signal_channel_names(ssys, outputs)
+    nu = length(unames)
     if nu > 0 && size(matrices.B, 2) == 2nu
         # This indicates that input derivatives are present
         # Derivative column nu + i is the derivative of input i (see the MTK docs for
@@ -171,18 +225,15 @@ function RobustAndOptimalControl.named_ss(
     else
         lsys = ss(matrices...)
     end
-    pind = [ModelingToolkit.parameter_index(ssys, i) for i in ModelingToolkit.inputs(ssys)]
-    x0 = xpt.x
-    u0 = [xpt.p[pi] for pi in pind] 
-    xu = (; x = x0, u = u0)
+    xu = (; x = xpt.x, u = input_operating_point(ssys, xpt))
     extra = Dict(:operating_point => xu)
     # If simple_infeigs=false, the system might have been reduced and the state names might not match the original system.
     x_names = get_x_names(lsys, ssys; descriptor, simple_infeigs, balance)
     nsys = named_ss(
         lsys;
         x = x_names,
         u = unames,
-        y = symstr.(outputs),
+        y = ynames,
         name = string(Base.nameof(sys)),
         extra,
     )
@@ -311,10 +362,9 @@ function named_sensitivity_function(
             inp
         end
     end
-    nu = length(inputs)
+    unames = signal_channel_names(sys, inputs)
+    nu = length(unames)
     matrices, ssys, xpt = fun(sys, inputs, args...; DEFAULT_LINEARIZE_KWARGS..., kwargs...)
-    symstr(x) = Symbol(x isa AnalysisPoint ? x.name : string(x))
-    unames = symstr.(inputs)
     fm(x) = convert(Matrix{Float64}, x)
     if nu > 0 && size(matrices.B, 2) == 2nu
         # This indicates that input derivatives are present
@@ -329,8 +379,7 @@ function named_sensitivity_function(
         lsys = ss(matrices...)
     end
     x_names = get_x_names(lsys, ssys; descriptor, simple_infeigs, balance)
-    u0 = [xpt.p[ModelingToolkit.parameter_index(ssys, i)] for i in ModelingToolkit.inputs(ssys)]
-    xu = (; x=xpt.x, u = u0)
+    xu = (; x = xpt.x, u = input_operating_point(ssys, xpt))
     extra = Dict(:operating_point => xu)
     nsys = named_ss(
         lsys;

test/test_ODESystem.jl

@@ -405,4 +405,58 @@ w_test = [1e-5, 1, 100]
 @test freqresp(G, w_test) ≈ freqresp(G3, w_test)
 @test freqresp(G, w_test) ≈ freqresp(G4, w_test)
 
-# bodeplot([G, G2, G3, G4, G5], exp10.(LinRange(-2, 2, 200)))
+# bodeplot([G, G2, G3, G4, G5], exp10.(LinRange(-2, 2, 200)))
+## Analysis points on array-valued connections
+# Each element of an array-valued analysis point is one channel of the linearization. The
+# channel count must not be taken from the number of analysis points: with a 2-vector input
+# and `nu = 1`, the second column of `B` was treated as the derivative of the first.
+@testset "Array-valued analysis points" begin
+    tv = ModelingToolkit.t_nounits
+    Dv = ModelingToolkit.D_nounits
+    @component function VecPlant(; name)
+        @variables x(tv)[1:2] = zeros(2) u(tv)[1:2] [input = true] y(tv)[1:2] [output = true]
+        eqs = [Dv(x[1]) ~ -x[1] + u[1], Dv(x[2]) ~ -2x[2] + 3u[2], y[1] ~ x[1], y[2] ~ x[2]]
+        System(eqs, tv; name)
+    end
+    @component function VecGain(; name)
+        @variables u(tv)[1:2] [input = true] y(tv)[1:2] [output = true]
+        @parameters k[1:2] = [1.0, 5.0]
+        System([y[1] ~ -k[1] * u[1], y[2] ~ -k[2] * u[2]], tv; name)
+    end
+    @named vplant = VecPlant()
+    @named vgain = VecGain()
+    vsys = System(
+        [connect(vplant.y, :y, vgain.u), connect(vgain.y, :u, vplant.u)], tv;
+        systems = [vplant, vgain], name = :vsys
+    )
+    Pref = ss([-1.0 0; 0 -2], [1.0 0; 0 3], I(2), 0)
+    K = diagm([1.0, 5.0])
+    # From the input analysis point to the output analysis point of the closed loop
+    Gref = feedback(Pref, K)
+    for allow_input_derivatives in (false, true)
+        Gv = named_ss(vsys, [:u], [:y]; allow_input_derivatives)
+        @test Gv.u == [:u1, :u2]
+        @test Gv.y == [:y1, :y2]
+        @test length(Gv.extra[:operating_point].u) == 2
+        @test tf(minreal(Gv.sys)) ≈ tf(Gref)
+        # The channels of an analysis point are selected by its name
+        @test size(Gv[:y, :u]) == (2, 2)
+    end
+    # Array-valued input analysis point with scalar outputs
+    Gv = named_ss(vsys, [:u], [vplant.y[1], vplant.y[2]])
+    @test Gv.u == [:u1, :u2]
+    @test tf(minreal(Gv.sys)) ≈ tf(Gref)
+
+    Sv = get_named_sensitivity(vsys, :u)
+    @test Sv.u == Sv.y == [:u1, :u2]
+    # With the negative-feedback gain K, the sensitivity at the plant input is (I + K P)^-1
+    @test tf(minreal(Sv.sys)) ≈ tf(feedback(ss(I(2)), K * Pref))
+
+    @test ControlSystemsMTK.signal_channel_names(vsys, [:u, :y]) == [:u1, :u2, :y1, :y2]
+    Gs = named_ss(vsys, [:u], [vplant.y[1]])
+    @test ControlSystemsMTK.signal_channel_names(vsys, [:u, vplant.y[1]]) == [Gs.u; Gs.y]
+    # The channels of an array-valued signal `x` coincide with the name of a scalar signal `x1`
+    @test_throws ArgumentError ControlSystemsMTK.channel_names([:x, :x1], [2, 1])
+    @test_throws r"`x1` \(from `x` and `x1`\)" ControlSystemsMTK.channel_names([:x, :x1], [2, 1])
+    @test ControlSystemsMTK.channel_names([:x, :x3], [2, 1]) == [:x1, :x2, :x3]
+end