Use LinearizationOpPoint for trajectory_ss

Project.toml

@@ -1,13 +1,14 @@
 name = "ControlSystemsMTK"
 uuid = "687d7614-c7e5-45fc-bfc3-9ee385575c88"
+version = "2.7.0"
 authors = ["Fredrik Bagge Carlson"]
-version = "2.6.0"
 
 [deps]
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+ModelingToolkitBase = "7771a370-6774-4173-bd38-47e70ca0b839"
 ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
 MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
 RobustAndOptimalControl = "21fd56a4-db03-40ee-82ee-a87907bee541"
@@ -17,7 +18,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 [compat]
 ControlSystemsBase = "1.0.1"
 DataInterpolations = "5, 6, 7, 8"
-ModelingToolkit = "11"
+ModelingToolkit = "11.7"
 ModelingToolkitStandardLibrary = "2"
 MonteCarloMeasurements = "1.1"
 RobustAndOptimalControl = "0.4.14"
@@ -34,3 +35,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["Test", "ControlSystems", "GenericLinearAlgebra", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqNonlinearSolve"]
+
+[sources]
+ModelingToolkit = {url = "https://github.com/SciML/ModelingToolkit.jl", rev = "as/linearize-op"}
+ModelingToolkitBase = {url = "https://github.com/SciML/ModelingToolkit.jl", rev = "as/linearize-op", subdir = "lib/ModelingToolkitBase"}

docs/Project.toml

@@ -4,6 +4,7 @@ ControlSystemsMTK = "687d7614-c7e5-45fc-bfc3-9ee385575c88"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+ModelingToolkitBase = "7771a370-6774-4173-bd38-47e70ca0b839"
 ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
 MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
 OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
@@ -12,3 +13,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RobustAndOptimalControl = "21fd56a4-db03-40ee-82ee-a87907bee541"
 SymbolicControlSystems = "886cb795-8fd3-4b11-92f6-8071e46f71c5"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[sources]
+ModelingToolkit = {url = "https://github.com/SciML/ModelingToolkit.jl", rev = "as/linearize-op"}
+ModelingToolkitBase = {url = "https://github.com/SciML/ModelingToolkit.jl", rev = "as/linearize-op", subdir = "lib/ModelingToolkitBase"}

docs/src/api.md

@@ -19,14 +19,13 @@ get_named_sensitivity
 get_named_comp_sensitivity
 get_named_looptransfer
 ModelingToolkit.linearize_symbolic
-ModelingToolkitStandardLibrary.Blocks.get_sensitivity
-ModelingToolkitStandardLibrary.Blocks.get_comp_sensitivity
-ModelingToolkitStandardLibrary.Blocks.get_looptransfer
+ModelingToolkit.get_sensitivity
+ModelingToolkit.get_comp_sensitivity
+ModelingToolkit.get_looptransfer
 ModelingToolkitStandardLibrary.Blocks.StateSpace
 RobustAndOptimalControl.ss2particles
 ControlSystemsBase.StateSpace
 SymbolicControlSystems.ccode
 SymbolicControlSystems.print_c_array
-ModelingToolkit.reorder_states
-ControlSystemsMTK.fuzz
+ModelingToolkit.reorder_unknowns
 ```

docs/src/batch_linearization.md

@@ -123,7 +123,7 @@ for C in Cs
         connect(Ci.output, duffing.u)
     ]
     @named closed_loop = System(eqs, t, systems=[duffing, Ci, fb, ref, F])
-    prob = ODEProblem(structural_simplify(closed_loop), [F.x => 0, F.xd => 0], (0.0, 8.0))
+    prob = ODEProblem(mtkcompile(closed_loop), [F.x => 0, F.xd => 0], (0.0, 8.0))
     sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8)
     plot!(sol, idxs=[duffing.y.u, duffing.u.u], layout=2, lab="")
 end
@@ -137,8 +137,8 @@ eqs = [
     connect(duffing.y, :v, Cgs.scheduling_input) # Don't forget to connect the scheduling variable!
 ]
 @named closed_loop = System(eqs, t, systems=[duffing, Cgs, fb, ref, F])
-prob = ODEProblem(structural_simplify(closed_loop), [F.xd => 0], (0.0, 8.0))
-sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8, initializealg=SciMLBase.NoInit(), dtmax=0.01)
+prob = ODEProblem(mtkcompile(closed_loop), [F.x => 0, F.xd => 0], (0.0, 8.0))
+sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8, dtmax=0.01)
 plot!(sol, idxs=[duffing.y.u, duffing.u.u], l=(2, :red), lab="Gain scheduled")
 plot!(sol, idxs=F.output.u, l=(1, :black, :dash, 0.5), lab="Ref")
 ```
@@ -164,10 +164,10 @@ bodeplot(Ps2, w, legend=false)
 ```
 Not how the closed-loop system changes very little along the trajectory, this is a good indication that the gain-scheduled controller is able to make the system appear linear.
 
-Internally, [`trajectory_ss`](@ref) works very much the same as [`batch_ss`](@ref), but constructs operating points automatically along the trajectory. This requires that the solution contains the states of the simplified system, accessible through the `idxs` argument like `sol(t, idxs=x)`. By linearizing the same system as we simulated, we ensure that this condition holds, doing so requires that we specify the inputs and outputs as analysis points rather than as variables.
+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. When linearizing the plant, we disconnect the controller input by passing `loop_openings=[closed_loop.u]`, and when linearizing the controller, we have various options for disconnecting the 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 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:
 - 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]);
+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));
 
 closed_loopsy = feedback.(plants .* controllersy)
 bodeplot(closed_loopsy, w; title="Loop open at y", kwargs...)
 ```
-When we open the loop at `u`, we get a slightly different result:
+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:
 ```@example BATCHLIN
-controllersu, ssu = trajectory_ss(closed_loop, closed_loop.r, closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.u]);
+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));
 
 closed_loopsu = feedback.(plants .* controllersu)
 bodeplot(closed_loopsu, w; title="Loop open at u", kwargs...)
 ```
-In this case, all static gains are 1. A similar result is obtained by explicitly breaking the scheduling feedback `v` in addition to an opening of `y`:
+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`:
 ```@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]);
+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));
 
 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]);
+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));
 
 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,20 +218,21 @@ plot(
     bodeplot(controllersv, w, legend=false, plotphase=false, title="Loop open at v and y"),
 )
 ```
-if we open at both `y` and `v` or we open at `u`, we get controllers for the different values of the scheduling variable, and the corresponding measurement feedback (which is the same as the scheduling variable in this case).
-```@example BATCHLIN
-using Test
-@test all(sminreal.(controllersv) .== sminreal.(controllersu))
-```
+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".
 
-However, if we only open at `y` we get controller linearizations that _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
+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
 sminreal(controllersy[end]).x
 ```
 notice how the state of the plant is included in the controller, this is an indication that we didn't fully isolate the controller during the linearizaiton. If we do the same thing for the controller with the loop opened at `u`, we see that the state of the plant is not included in the controller:
 ```@example BATCHLIN
 sminreal(controllersu[end]).x
 ```
+Opening at `u` thus yields a controller with fewer states than opening only at `y`, since the latter remains coupled to the plant through the scheduling feedback:
+```@example BATCHLIN
+using Test
+@test sminreal(controllersu[end]).nx < sminreal(controllersy[end]).nx
+```
 The call to `sminreal` is important here, it removes the states that are not needed to represent the input-output map of the system. The state of the full model, including the plant state, is present in the linearization before this call. 
 
 

docs/src/index.md

@@ -178,7 +178,11 @@ model = SystemModel() |> complete
 ### Numeric linearization
 We can linearize this model numerically using `named_ss`, this produces a `NamedStateSpace{Continuous, Float64}`
 ```@example LINEAIZE_SYMBOLIC
-op = Dict(model.inertia1.flange_b.phi => 0.0, model.torque.tau.u => 0)
+op = Dict(
+    model.inertia1.flange_b.phi => 0.0,
+    model.inertia2.flange_a.phi => 0.0,
+    model.torque.tau.u => 0,
+)
 lsys = named_ss(model, [model.torque.tau.u], [model.inertia1.phi, model.inertia2.phi]; op)
 ```
 ### Symbolic linearization
@@ -192,7 +196,7 @@ symbolic_sys = ss(mats.A, mats.B, mats.C, mats.D)
 That's pretty cool, but even nicer is to generate some code for this symbolic system. Below, we use `build_function` to generate a function that takes a numeric vector `x` representing the values of the state, and a vector of parameters, and returns a `StaticStateSpace{Continuous, Float64}`. We pass the keyword argument `force_SA=true` to `build_function` to get an allocation-free function. 
 
 ```@example LINEAIZE_SYMBOLIC
-defs = ModelingToolkit.defaults(simplified_sys)
+defs = ModelingToolkit.initial_conditions(simplified_sys)
 defs = merge(Dict(unknowns(model) .=> 0), defs)
 x = ModelingToolkit.get_u0(simplified_sys, defs) # Extract the default state and parameter values
 pars = ModelingToolkit.get_p(simplified_sys, defs, split=false)