Default loop-opened signals in trajectory_ss to their trajectory values

Project.toml

@@ -1,7 +1,7 @@
 name = "ControlSystemsMTK"
 uuid = "687d7614-c7e5-45fc-bfc3-9ee385575c88"
 authors = ["Fredrik Bagge Carlson"]
-version = "2.8.1"
+version = "2.9.0"
 
 [deps]
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"

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

src/ode_system.jl

@@ -538,21 +538,24 @@ 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)).")
 
     input_names = reduce(vcat, getproperty.(ap.outputs, :u) for ap in vcat(inputs))
     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
@@ -566,6 +569,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
 

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))