From 9f1cc2795a79fff46588969ec8e4793934a667cd Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Thu, 6 Aug 2026 12:00:19 +0000 Subject: [PATCH] Default loop-opened signals in trajectory_ss to their trajectory values MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Opening a loop turns the opened signal into a parameter whose value must be supplied in the operating point. Previously the tutorial held these at 0, which linearizes at the wrong point for systems nonlinear in the opened signal — visible as identical controller linearizations along the whole trajectory in the batch-linearization tutorial. trajectory_ss now maps each opened signal to itself as a symbolic operating- point value, which ModelingToolkit.LinearizationOpPoint resolves from the loop-closed solution at each time point (requires the corresponding MTK support for symbolic op values). User-supplied op entries take precedence. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01SWPXe76HwiugWWyLP7BYrL --- Project.toml | 2 +- docs/src/batch_linearization.md | 20 +++++++++------ src/ode_system.jl | 45 +++++++++++++++++++++++++++++---- test/test_batchlin.jl | 8 ++++++ 4 files changed, 61 insertions(+), 14 deletions(-) diff --git a/Project.toml b/Project.toml index 9cd1753..477fb38 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ControlSystemsMTK" uuid = "687d7614-c7e5-45fc-bfc3-9ee385575c88" authors = ["Fredrik Bagge Carlson"] -version = "2.8.0" +version = "2.9.0" [deps] ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e" diff --git a/docs/src/batch_linearization.md b/docs/src/batch_linearization.md index 4116995..a11e2b3 100644 --- a/docs/src/batch_linearization.md +++ b/docs/src/batch_linearization.md @@ -167,7 +167,7 @@ Not how the closed-loop system changes very little along the trajectory, this is Internally, [`trajectory_ss`](@ref) works very much the same as [`batch_ss`](@ref), but constructs operating points automatically along the trajectory using `ModelingToolkit.LinearizationOpPoint`. The operating points are extracted from the differential states and parameters of the solution. We specify the inputs and outputs as analysis points to properly define the linearization interface. -We can replicate the figure above by linearizing the plant and the controller individually, by providing the `loop_openings` argument. Opening a loop breaks the corresponding connection and turns the opened signal into a free input. Its operating-point value is not implied by the solution, so it must be supplied explicitly through the `op` argument (a `Dict` mapping the opened signal to its value); here we hold each opened signal at `0`. When linearizing the plant, we disconnect the controller output by passing `loop_openings=[closed_loop.u]`, and when linearizing the controller, we have various options for disconnecting the plant: +We can replicate the figure above by linearizing the plant and the controller individually, by providing the `loop_openings` argument. Opening a loop breaks the corresponding connection and turns the opened signal into a free input. Its operating-point value is not implied by the equations of the system anymore, but since we are linearizing around a trajectory of the loop-closed system, the natural value is the one the signal has in the solution at each time point. This is what `trajectory_ss` does by default: each opened signal is linearized around its own value in the loop-closed solution. To linearize around some other value, pass it through the `op` argument, e.g., `op = Dict(opened_signal => 0)`. When linearizing the plant, we disconnect the controller output by passing `loop_openings=[closed_loop.u]`, and when linearizing the controller, we have various options for disconnecting the plant: - Break the connection from plant output to controller input by passing `loop_openings=[closed_loop.y]` - Break the connection between the controller and the plant input by passing `loop_openings=[closed_loop.u]` - Break the connection `y` as well as the scheduling variable `v` (which is another form of feedback) by passing `loop_openings=[closed_loop.y, closed_loop.v]` @@ -176,21 +176,21 @@ We will explore these options below, starting with the first option, breaking th ```@example BATCHLIN kwargs = (; adaptive=false, legend=false) plants, _ = trajectory_ss(closed_loop, closed_loop.u, closed_loop.y, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.u]); -controllersy, ssy, ops3, resolved_ops3 = trajectory_ss(closed_loop, closed_loop.r, closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.y], op=Dict(fb.input2.u => 0)); +controllersy, ssy, ops3, resolved_ops3 = trajectory_ss(closed_loop, closed_loop.r, closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.y]); closed_loopsy = feedback.(plants .* controllersy) bodeplot(closed_loopsy, w; title="Loop open at y", kwargs...) ``` -When we open the loop at `u` instead, the plant input is held at `0` while the scheduling variable `v` remains connected, so the controller is linearized at the scheduling value along the trajectory and is fully isolated from the plant: +When we open the loop at `u` instead, the scheduling variable `v` remains connected, so the controller is linearized at the scheduling value along the trajectory and is fully isolated from the plant: ```@example BATCHLIN -controllersu, ssu = trajectory_ss(closed_loop, closed_loop.r, closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.u], op=Dict(duffing.u.u => 0)); +controllersu, ssu = trajectory_ss(closed_loop, closed_loop.r, closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.u]); closed_loopsu = feedback.(plants .* controllersu) bodeplot(closed_loopsu, w; title="Loop open at u", kwargs...) ``` -If we instead break the scheduling feedback `v` in addition to `y` (holding both at `0`), the controller is isolated from the plant but its scheduling is pinned to `0` rather than following the trajectory, so the result differs from opening at `u`: +We may also break the scheduling feedback `v` in addition to `y`. Since the opened signals follow their values in the loop-closed solution, the controller is scheduled along the trajectory also in this case, and the result coincides with opening at `u`: ```@example BATCHLIN -controllersv, ssv = trajectory_ss(closed_loop, closed_loop.r, closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.y, closed_loop.v], op=Dict(fb.input2.u => 0, Cgs.scheduling_input.u => 0)); +controllersv, ssv = trajectory_ss(closed_loop, closed_loop.r, closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.y, closed_loop.v]); closed_loopsv = feedback.(plants .* controllersv) bodeplot(closed_loopsv, w; title="Loop open at v and y", kwargs...) @@ -200,7 +200,7 @@ We have thus far treated the controller as a SISO system, but we could also view ```@example BATCHLIN plants_mimo, _ = trajectory_ss(closed_loop, closed_loop.u, [closed_loop.y, closed_loop.v], sol; t=timepoints, verbose=true, loop_openings=[closed_loop.u]); -controllers_mimo, ssm = trajectory_ss(closed_loop, [closed_loop.r, closed_loop.v], closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.u], op=Dict(duffing.u.u => 0)); +controllers_mimo, ssm = trajectory_ss(closed_loop, [closed_loop.r, closed_loop.v], closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.u]); closed_loops_mimo = feedback.(controllers_mimo .* plants_mimo) # Look at complementary sensitivity function in the input, since this is a SISO system bodeplot(closed_loops_mimo, w; title="Loop open at MIMO", kwargs...) @@ -218,7 +218,11 @@ plot( bodeplot(controllersv, w, legend=false, plotphase=false, title="Loop open at v and y"), ) ``` -Opening at `u` keeps the scheduling connection `v` intact, so we obtain the gain-scheduled controller evaluated at the scheduling value along the trajectory. Opening additionally at `v` and holding it at `0` instead pins the scheduling to `0`, which is why "Loop open at v and y" differs from "Loop open at u". +Opening at `u` keeps the scheduling connection `v` intact, so we obtain the gain-scheduled controller evaluated at the scheduling value along the trajectory. Opening at `v` and `y` yields the same controllers, since the opened scheduling signal follows its value in the loop-closed solution: +```@example BATCHLIN +using Test +@test all(isapprox(freqresp(cv, w), freqresp(cu, w), rtol=1e-6) for (cv, cu) in zip(controllersv, controllersu)) +``` If we only open at `y`, the scheduling feedback through `v` remains in place, so the controller linearizations _still contain the closed loop through the scheduling connection_ `v`. We can verify this by looking at what variables are present in the input-output map ```@example BATCHLIN diff --git a/src/ode_system.jl b/src/ode_system.jl index 93e8f63..5c65cbf 100644 --- a/src/ode_system.jl +++ b/src/ode_system.jl @@ -535,11 +535,12 @@ Operating points are extracted from the solution automatically using `ModelingTo - `outputs`: A vector of variables or analysis points. - `sol`: An ODE solution object. - `t`: Time points along the solution trajectory at which to linearize. The returned array of `StateSpace` objects will be of the same length as `t`. -- `op`: A `Dict` of additional operating-point values that are not available from `sol`. This is required when using `loop_openings`: opening a loop turns the opened signal into a parameter whose value is not implied by the solution, and it must be supplied here (typically set to `0`), e.g. `op = Dict(sys.opened_signal => 0)`. The values are merged into the solution-derived operating point at every time point. -- `kwargs`: Are sent to the linearization functions (e.g., `loop_openings`). +- `loop_openings`: A list of analysis points whose connections are broken during the linearization. Each opened signal becomes a parameter whose value is by default taken from the loop-closed solution `sol` at each time point, such that the linearization is performed around the trajectory also for the opened signals. To linearize around some other value for an opened signal, supply it through `op`. +- `op`: A `Dict` of operating-point values overriding those obtained from `sol`, e.g. `op = Dict(sys.opened_signal => 0)`. Symbolic values are resolved from `sol` at every time point. The values are merged into the solution-derived operating point at every time point. +- `kwargs`: Are sent to the linearization functions. - `named`: If `true`, the returned systems will be of type `NamedStateSpace`, otherwise they will be of type `StateSpace`. """ -function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), op = Dict(), allow_input_derivatives = false, verbose = true, named = true, kwargs...) +function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), op = Dict(), loop_openings = [], allow_input_derivatives = false, verbose = true, named = true, kwargs...) maximum(t) > maximum(sol.t) && @warn("The maximum time in `t`: $(maximum(t)), is larger than the maximum time in `sol.t`: $(maximum(sol.t)).") minimum(t) < minimum(sol.t) && @warn("The minimum time in `t`: $(minimum(t)), is smaller than the minimum time in `sol.t`: $(minimum(sol.t)).") @@ -547,9 +548,11 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), op = Dict output_names = reduce(vcat, ap.input.u for ap in vcat(outputs)) # Use LinearizationOpPoint to let MTK extract operating points from the solution. - # `op` supplies values not available from `sol` (e.g. loop-opening parameters). + # `op` supplies values not available from `sol` (e.g. loop-opening parameters); + # opened signals default to their own value in the loop-closed solution. + op = merge(default_loop_opening_op(sys, loop_openings), Dict(Symbolics.unwrap(k) => v for (k, v) in pairs(op))) oppoint = ModelingToolkit.LinearizationOpPoint(sol, t; op) - lins, ssys, resolved_ops = linearize(sys, inputs, outputs; op = oppoint, allow_input_derivatives, DEFAULT_LINEARIZE_KWARGS..., kwargs...) + lins, ssys, resolved_ops = linearize(sys, inputs, outputs; op = oppoint, allow_input_derivatives, loop_openings, DEFAULT_LINEARIZE_KWARGS..., kwargs...) named_linsystems = map(lins) do l if named @@ -563,6 +566,38 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), op = Dict (; linsystems = named_linsystems, ssys, ops = resolved_ops, resolved_ops) end +""" + default_loop_opening_op(sys, loop_openings) + +Construct default operating-point values for the signals opened by `loop_openings`. +Opening a loop turns the opened signal into a parameter whose value is not implied by the +rest of the system. Each opened signal is mapped to itself as a symbolic value, which +`ModelingToolkit.LinearizationOpPoint` resolves from the loop-closed solution at each time +point — the linearization is thus performed around the trajectory also for the opened +signals. +""" +function default_loop_opening_op(sys, loop_openings) + defop = Dict{Any, Any}() + for ap in ModelingToolkit.canonicalize_ap(sys, collect(loop_openings)) + ap isa ModelingToolkit.AnalysisPoint || continue + ap.outputs === nothing && continue # AP specified by name only, no default possible + for out in ap.outputs + v = strip_root_namespace(sys, ModelingToolkit.ap_var(out)) + defop[v] = v + end + end + defop +end + +# Analysis points obtained through `getproperty` on the root system have their variables +# namespaced with the root system name, while operating points (and solution indexing) use +# root-namespace-stripped variables. +function strip_root_namespace(sys, v) + parts = ModelingToolkit.namespace_hierarchy(Symbolics.getname(v)) + (length(parts) > 1 && parts[1] == nameof(sys)) || return v + Symbolics.rename(v, Symbol(join(parts[2:end], Symbolics.NAMESPACE_SEPARATOR))) +end + "_max_100(t) = length(t) > 100 ? range(extrema(t)..., 100) : t" _max_100(t) = length(t) > 100 ? range(extrema(t)..., 100) : t diff --git a/test/test_batchlin.jl b/test/test_batchlin.jl index e7cc035..7b9dbbf 100644 --- a/test/test_batchlin.jl +++ b/test/test_batchlin.jl @@ -79,3 +79,11 @@ op = Dict(u.u => 0) Ps2, ssys = trajectory_ss(closed_loop, closed_loop.r, closed_loop.y, sol; t=time) @test length(Ps2) == length(time) # bodeplot(Ps2) + +## Loop openings: opened signals default to their value in the loop-closed solution +controllers, _ = trajectory_ss(closed_loop, closed_loop.r, closed_loop.y, sol; t=time, loop_openings=[closed_loop.y]) +@test length(controllers) == length(time) +# The default is equivalent to explicitly mapping the opened signal to its driving signal +controllers2, _ = trajectory_ss(closed_loop, closed_loop.r, closed_loop.y, sol; t=time, loop_openings=[closed_loop.y], op=Dict(fb.input2.u => duffing.y.u)) +w_test = exp10.(LinRange(-2, 2, 30)) +@test all(isapprox(freqresp(c1, w_test), freqresp(c2, w_test), rtol=1e-10) for (c1, c2) in zip(controllers, controllers2))