From f6dfcaa396379b1ade5e5e1bfa3a611db2c62fef Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Tue, 28 Jan 2025 06:32:51 +0100 Subject: [PATCH 1/2] handle analysis points in addition to symbols --- src/ode_system.jl | 73 ++++++++++++++++++++++------------------------- 1 file changed, 34 insertions(+), 39 deletions(-) diff --git a/src/ode_system.jl b/src/ode_system.jl index 563e55d..04acb5a 100644 --- a/src/ode_system.jl +++ b/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,40 +170,36 @@ 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)) unames = symstr.(inputs) @@ -276,23 +274,20 @@ 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)) unames = symstr.(inputs) From e4a0cde18027074bdf040a7560e6bef6346711c8 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Tue, 28 Jan 2025 07:07:47 +0100 Subject: [PATCH 2/2] use ap access style from MTK --- docs/src/batch_linearization.md | 4 ++-- src/ode_system.jl | 4 ++-- test/test_ODESystem.jl | 6 +++--- test/test_batchlin.jl | 2 +- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/src/batch_linearization.md b/docs/src/batch_linearization.md index bf68d2c..57bcb5e 100644 --- a/docs/src/batch_linearization.md +++ b/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) ``` diff --git a/src/ode_system.jl b/src/ode_system.jl index 04acb5a..6d8995f 100644 --- a/src/ode_system.jl +++ b/src/ode_system.jl @@ -201,7 +201,7 @@ function RobustAndOptimalControl.named_ss( 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 @@ -289,7 +289,7 @@ function named_sensitivity_function( 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 diff --git a/test/test_ODESystem.jl b/test/test_ODESystem.jl index 7c0546f..b2fbaf8 100644 --- a/test/test_ODESystem.jl +++ b/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 diff --git a/test/test_batchlin.jl b/test/test_batchlin.jl index 94b2eb7..411f89b 100644 --- a/test/test_batchlin.jl +++ b/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)