Merge pull request #70 from JuliaControl/ap

handle analysis points in addition to symbols

Author: baggepinnen

docs/src/batch_linearization.md

@@ -152,10 +152,10 @@ code = SymbolicControlSystems.print_c_array(stdout, Cs_disc[1:7:end], xs[1:7:end
 The generated code starts by defining the interpolation vector `xs`, this variable is called `Cs_interp_vect` in the generated code. The code then defines all the ``A`` matrices as a 3-dimensional array, followed by a function that performs the interpolation `interpolate_Cs_A`. This function takes the output array as the first argument, a pointer to the 3D array with interpolation matrices, the interpolation vector as well as the interpolation variable `t`, in this document called ``v``. The same code is then repeated for the matrices ``B,C,D`` as well if they require interpolation (if they are all the same, no interpolation code is written). 
 
 ## Linearize around a trajectory
-We can linearize around a trajectory obtained from `solve` using the function [`trajectory_ss`](@ref). We provide it with a vector of time points along the trajectory at which to linearize, and in this case we specify the inputs and outputs to linearize between as analysis points `:r` and `:y`.
+We can linearize around a trajectory obtained from `solve` using the function [`trajectory_ss`](@ref). We provide it with a vector of time points along the trajectory at which to linearize, and in this case we specify the inputs and outputs to linearize between as analysis points `r` and `y`.
 ```@example BATCHLIN
 timepoints = 0:0.01:8
-Ps2, ssys = trajectory_ss(closed_loop, :r, :y, sol; t=timepoints)
+Ps2, ssys = trajectory_ss(closed_loop, closed_loop.r, closed_loop.y, sol; t=timepoints)
 bodeplot(Ps2, w, legend=false)
 ```
 

src/ode_system.jl

@@ -1,3 +1,5 @@
+using ModelingToolkit: AnalysisPoint
+const AP = Union{Symbol, AnalysisPoint}
 import ModelingToolkitStandardLibrary.Blocks as Blocks
 conn = ModelingToolkit.connect
 t = Blocks.t
@@ -168,42 +170,38 @@ function RobustAndOptimalControl.named_ss(
     kwargs...,
 )
 
-    if isa(inputs,  Symbol)
-        nu = 1
-    else
-        inputs = map(inputs) do inp
-            if inp isa ODESystem
-                @variables u(t)
-                if u ∈ Set(unknowns(inp))
-                    inp.u
-                else
-                    error("Input $(inp.name) is an ODESystem and not a variable")
-                end
+    inputs = vcat(inputs)
+    outputs = vcat(outputs)
+
+    inputs = map(inputs) do inp
+        if inp isa ODESystem
+            @variables u(t)
+            if u ∈ Set(unknowns(inp))
+                inp.u
             else
-                inp
+                error("Input $(inp.name) is an ODESystem and not a variable")
             end
+        else
+            inp
         end
-        nu = length(inputs)
     end
-    if isa(outputs,  Symbol)
-        ny = 1
-    else
-        outputs = map(outputs) do out
-            if out isa ODESystem
-                @variables u(t)
-                if u ∈ Set(unknowns(out))
-                    out.u
-                else
-                    error("Outut $(out.name) is an ODESystem and not a variable")
-                end
+    nu = length(inputs)
+
+    outputs = map(outputs) do out
+        if out isa ODESystem
+            @variables u(t)
+            if u ∈ Set(unknowns(out))
+                out.u
             else
-                out
+                error("Outut $(out.name) is an ODESystem and not a variable")
             end
+        else
+            out
         end
-        ny = length(outputs)
     end
+    ny = length(outputs)
     matrices, ssys = ModelingToolkit.linearize(sys, inputs, outputs; kwargs...)
-    symstr(x) = Symbol(string(x))
+    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
@@ -276,25 +274,22 @@ function named_sensitivity_function(
     kwargs...,
 )
 
-    if isa(inputs,  Symbol)
-        nu = 1
-    else
-        inputs = map(inputs) do inp
-            if inp isa ODESystem
-                @variables u(t)
-                if u ∈ Set(unknowns(inp))
-                    inp.u
-                else
-                    error("Input $(inp.name) is an ODESystem and not a variable")
-                end
+    inputs = vcat(inputs)
+    inputs = map(inputs) do inp
+        if inp isa ODESystem
+            @variables u(t)
+            if u ∈ Set(unknowns(inp))
+                inp.u
             else
-                inp
+                error("Input $(inp.name) is an ODESystem and not a variable")
             end
+        else
+            inp
         end
-        nu = length(inputs)
     end
+    nu = length(inputs)
     matrices, ssys = fun(sys, inputs, args...; kwargs...)
-    symstr(x) = Symbol(string(x))
+    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

test/test_ODESystem.jl

@@ -266,14 +266,14 @@ eqs = [connect(r.output, F.input)
     connect(F.output, sys_inner.add.input1)]
 sys_outer = ODESystem(eqs, t, systems = [F, sys_inner, r], name = :outer)
 
-matrices, _ = Blocks.get_sensitivity(sys_outer, [:inner_plant_input, :inner_plant_output])
+matrices, _ = Blocks.get_sensitivity(sys_outer, [sys_outer.inner.plant_input, sys_outer.inner.plant_output])
 S = ss(matrices...)
 
-Sn = get_named_sensitivity(sys_outer, [:inner_plant_input, :inner_plant_output])
+Sn = get_named_sensitivity(sys_outer, [sys_outer.inner.plant_input, sys_outer.inner.plant_output])
 
 @test S == Sn.sys
 
-@test Sn.u == Sn.y == [:inner_plant_input, :inner_plant_output]
+@test Sn.u == Sn.y == [:outer₊inner₊plant_input, :outer₊inner₊plant_output] == [:outer₊inner₊plant_input, :outer₊inner₊plant_output]
 
 
 ## Test connector names

test/test_batchlin.jl

@@ -72,6 +72,6 @@ sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8)
 
 time = 0:0.1:8
 inputs, outputs = [duffing.u.u], [duffing.y.u]
-Ps2, ssys = trajectory_ss(closed_loop, :r, :y, sol; t=time)
+Ps2, ssys = trajectory_ss(closed_loop, closed_loop.r, closed_loop.y, sol; t=time)
 @test length(Ps2) == length(time)
 # bodeplot(Ps2)