From 3a5755fb56a8f4f1f00952d1fb28fd0e4040f2d7 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 13 Apr 2026 07:17:27 +0200
Subject: [PATCH 01/12] Use LinearizationOpPoint for trajectory_ss, remove
 robust_sol_getindex

Replace manual operating point extraction in `trajectory_ss` with
MTK's `LinearizationOpPoint` API (SciML/ModelingToolkit.jl#4443).
This removes the fragile `robust_sol_getindex` helper and simplifies
the implementation.

- Use `_build_op_from_solution` to extract differential states + parameters
- Supplement with linearization system unknowns from the solution
- Remove `robust_sol_getindex` (no longer needed)
- Bump version to 2.7.0, require MTK >= 11.7
- Update docs narrative for trajectory_ss

Closes SciML/ModelingToolkit.jl#4159

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
---
 Project.toml                    |   4 +-
 docs/src/batch_linearization.md |   4 +-
 src/ode_system.jl               | 114 +++++++++++---------------------
 3 files changed, 43 insertions(+), 79 deletions(-)

diff --git a/Project.toml b/Project.toml
index 0ed5462..0fbf7b7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 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"
@@ -17,7 +17,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"
diff --git a/docs/src/batch_linearization.md b/docs/src/batch_linearization.md
index 0cded72..4da6462 100644
--- a/docs/src/batch_linearization.md
+++ b/docs/src/batch_linearization.md
@@ -164,7 +164,7 @@ 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:
@@ -221,7 +221,7 @@ plot(
 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))
+@test all(isapprox.(sminreal.(controllersv), sminreal.(controllersu), atol=1e-10))
 ```
 
 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
diff --git a/src/ode_system.jl b/src/ode_system.jl
index a8fad8f..a56e7ce 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -301,7 +301,7 @@ function named_sensitivity_function(
         end
     end
     nu = length(inputs)
-    matrices, ssys = fun(sys, inputs, args...; kwargs...)
+    matrices, ssys, xpt = fun(sys, inputs, args...; kwargs...)
     symstr(x) = Symbol(x isa AnalysisPoint ? x.name : string(x))
     unames = symstr.(inputs)
     fm(x) = convert(Matrix{Float64}, x)
@@ -314,12 +314,16 @@ function named_sensitivity_function(
         lsys = ss(matrices...)
     end
     x_names = get_x_names(lsys, ssys; descriptor, simple_infeigs, balance)
+    u0 = [xpt.p[ModelingToolkit.parameter_index(ssys, i)] for i in ModelingToolkit.inputs(ssys)]
+    xu = (; x=xpt.x, u = u0)
+    extra = Dict(:operating_point => xu)
     nsys = named_ss(
         lsys;
         x = x_names,
         u = unames,
         y = unames, #Symbol.("out_" .* string.(inputs)),
         name = string(Base.nameof(sys)),
+        extra,
     )
     RobustAndOptimalControl.set_extra!(nsys, :ssys, ssys)
     nsys
@@ -512,65 +516,69 @@ end
 
 Linearize `sys` around the trajectory `sol` at times `t`. Returns a vector of `StateSpace` objects and the simplified system.
 
+Operating points are extracted from the solution automatically using [`ModelingToolkit.LinearizationOpPoint`](@ref).
+
 # Arguments:
 - `inputs`: A vector of variables or analysis points.
 - `outputs`: A vector of variables or analysis points.
-- `sol`: An ODE solution object. This solution must contain the states of the simplified system, accessible through the `idxs` argument like `sol(t, idxs=x)`.
+- `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`.
 - `fuzzer`: A function that takes an operating point dictionary and returns an array of "fuzzed" operating points. This is useful for adding noise/uncertainty to the operating points along the trajectory. See [`ControlSystemsMTK.fuzz`](@ref) for such a function.
-- `verbose`: If `true`, print warnings for variables that are not found in `sol`.
-- `kwargs`: Are sent to the linearization functions.
+- `kwargs`: Are sent to the linearization functions (e.g., `loop_openings`).
 - `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), allow_input_derivatives = false, fuzzer = nothing, 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)).")
-    # NOTE: we call linearization_funciton twice :( The first call is to get x=unknowns(ssys), the second call provides the operating points.
-    # lin_fun, ssys = linearization_function(sys, inputs, outputs; warn_initialize_determined = false, kwargs...)
-    lin_fun, ssys = linearization_function(sys, inputs, outputs; warn_empty_op = false, warn_initialize_determined = false, kwargs...)
-    x = unknowns(ssys)
 
-    # TODO: The value of the output (or input) of the input analysis points should be mapped to the perturbation vars
-    perturbation_vars = ModelingToolkit.inputs(ssys)
-    # original_inputs = reduce(vcat, unnamespace(ap) for ap in vcat(inputs)) # assuming all inputs are analysis points for now
+    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))
 
-    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 extract operating points from the solution.
+    # _build_op_from_solution gives differential states + parameters; we supplement
+    # with all unknowns of the linearization system to avoid initialization issues.
+    _extract_base_op(ti) = ModelingToolkit._build_op_from_solution(ModelingToolkit.LinearizationOpPoint(sol, ti))
+    tc = collect(t)
 
-    op_nothing = Dict(unknowns(sys) .=> nothing) # Remove all defaults present in the original system
+    # First pass: get the linearization system's unknowns
+    lin_fun0, ssys = linearization_function(sys, inputs, outputs; warn_initialize_determined=false, kwargs...)
+    lin_unknowns = unknowns(ssys)
     defs = ModelingToolkit.initial_conditions(sys)
-    ops = map(t) do ti
-        opsol = Dict(x => robust_sol_getindex(sol, ti, x, defs; verbose) for x in x)
-        # When the new behavior of Break is introduced, speficy the value for all inupts in ssys by `for x in [x; perturbation_vars]` on the line above
-        # opsolu = Dict(new_u => robust_sol_getindex(sol, ti, u, defs; verbose) for (new_u, u) in zip(perturbation_vars, original_inputs))
-        merge(op_nothing, opsol)
+
+    ops = map(tc) do ti
+        op = _extract_base_op(ti)
+        for x in lin_unknowns
+            haskey(op, x) && continue
+            try
+                op[x] = sol(ti, idxs=x)
+            catch
+                val = get(defs, x, nothing)
+                val !== nothing && (op[x] = val)
+            end
+        end
+        op
     end
+
     if fuzzer !== nothing
         opsv = map(ops) do op
             fuzzer(op)
         end
         ops = reduce(vcat, opsv)
-        t = repeat(t, inner = length(ops) ÷ length(t))
+        tc = repeat(tc, inner = length(ops) ÷ length(tc))
     end
-    lin_fun, ssys = linearization_function(sys, inputs, outputs; op=ops[1], kwargs...)#, initialization_abstol=1e-1, initialization_reltol=1e-1, kwargs...) # initializealg=ModelingToolkit.SciMLBase.NoInit()
-    # Main.lin_fun = lin_fun
-    # Main.op1 = ops[1]
-    # Main.ops = ops 
-    # equations(lin_fun.prob.f.initialization_data.initializeprob.f.sys)
-    # observed(lin_fun.prob.f.initialization_data.initializeprob.f.sys)
-    lins_ops = map(zip(ops, t)) do (op, t)
-        linearize(ssys, lin_fun; op, t, allow_input_derivatives)
-        # linearize(sys, inputs, outputs; op, t, allow_input_derivatives) # useful for debugging
+
+    lin_fun, ssys = linearization_function(sys, inputs, outputs; op=ops[1], t=tc[1], initialize=false, kwargs...)
+    lins_ops = map(zip(ops, tc)) do (op, ti)
+        linearize(ssys, lin_fun; op, t=ti, allow_input_derivatives)
     end
     lins = first.(lins_ops)
     resolved_ops = last.(lins_ops)
+
     named_linsystems = map(lins) do l
         if named
-            # Convert to a NamedStateSpace with the same names as the original system
             ynames = allunique(output_names) ? symstr.(output_names) : [Symbol(string(nameof(sys))*"_y$i") for i in 1:length(output_names)]
             unames = allunique(input_names) ? symstr.(input_names) : [Symbol(string(nameof(sys))*"_u$i") for i in 1:length(input_names)]
             nsys = named_ss(ss(l.A, l.B, l.C, l.D); name = string(Base.nameof(sys)), x = symstr.(unknowns(ssys)), u = unames, y = ynames)
-            # RobustAndOptimalControl.merge_nonunique_outputs(RobustAndOptimalControl.merge_nonunique_inputs(nsys))
         else
             ss(l.A, l.B, l.C, l.D)
         end
@@ -620,50 +628,6 @@ end
 MonteCarloMeasurements.vecindex(p::Symbolics.BasicSymbolic,i) = p
 issymbolic(x) = x isa Union{Symbolics.Num, Symbolics.BasicSymbolic}
 
-"""
-    robust_sol_getindex(sol, ti, x, defs; verbose = true)
-
-Extract symbolic variable `x` from ode solution `sol` at time `ti`. This operation may fail
-- If the variable is a dummy derivative that is not present in the solution. In this case, the value is reconstructed by derivative interpolation.
-- The var is not present at all, in this case, the default value in `defs` is returned.
-
-# Arguments:
-- `sol`: An ODESolution
-- `ti`: Time point
-- `defs`: A Dict with default values. 
-- `verbose`: Print a warning if the variable is not found in the solution.
-"""
-function robust_sol_getindex(sol, ti, x, defs; verbose = true)
-    try
-        return sol(ti, idxs=x)
-    catch
-        n = string((x))
-        if occursin("ˍt(", n)
-            n = split(n, "ˍt(")[1]
-            sp = split(n, '₊')
-            varname = sp[end]
-            local var
-            let t = Symbolics.arguments(Symbolics.unwrap(x))[1]
-                @variables var(t)
-            end
-            ModelingToolkit.@set! var.val.f.name = Symbol(varname)
-            namespaces = sp[1:end-1]
-            if !isempty(namespaces)
-                for ns in reverse(namespaces)
-                    var = ModelingToolkit.renamespace(Symbol(ns), var)
-                end
-            end
-            out = sol(ti, Val{1}, idxs=[Num(var)])[]
-            verbose && println("Could not find variable $x in solution, returning $(out) obtained through interpolation of $var.")
-            return out
-        end
-
-        val = get(defs, x, 0.0)
-        verbose && println("Could not find variable $x in solution, returning $val.")
-        return val
-    end
-end
-
 maybe_interp(interpolator, x, t) = allequal(x) ? x[1] : interpolator(x, t)
 
 """

From 0040cdf4b7059cdbe8b41532c00bba93676595d5 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 13 Apr 2026 07:24:53 +0200
Subject: [PATCH 02/12] Remove broken cross-reference to ModelingToolkit
 docstring

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
---
 src/ode_system.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/ode_system.jl b/src/ode_system.jl
index a56e7ce..a6e929c 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -516,7 +516,7 @@ end
 
 Linearize `sys` around the trajectory `sol` at times `t`. Returns a vector of `StateSpace` objects and the simplified system.
 
-Operating points are extracted from the solution automatically using [`ModelingToolkit.LinearizationOpPoint`](@ref).
+Operating points are extracted from the solution automatically using `ModelingToolkit.LinearizationOpPoint`.
 
 # Arguments:
 - `inputs`: A vector of variables or analysis points.

From 1fb8a8b8aacf3d5c0402de912ba084659cf21ec4 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 13 Apr 2026 08:30:28 +0200
Subject: [PATCH 03/12] Remove fuzz functionality, simplify trajectory_ss to
 use LinearizationOpPoint directly
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

trajectory_ss now delegates entirely to MTK's linearize with
LinearizationOpPoint(sol, t) — no manual op construction needed.

Removed ControlSystemsMTK.fuzz and all references to it.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
---
 docs/src/api.md   |  1 -
 src/ode_system.jl | 97 ++++++-----------------------------------------
 2 files changed, 11 insertions(+), 87 deletions(-)

diff --git a/docs/src/api.md b/docs/src/api.md
index 8197657..c051bb5 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -28,5 +28,4 @@ ControlSystemsBase.StateSpace
 SymbolicControlSystems.ccode
 SymbolicControlSystems.print_c_array
 ModelingToolkit.reorder_states
-ControlSystemsMTK.fuzz
 ```
\ No newline at end of file
diff --git a/src/ode_system.jl b/src/ode_system.jl
index a6e929c..c356195 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -492,7 +492,7 @@ centers, radii = fit_complex_perturbations(P * C, w; relative = false, nominal =
 nyquistcircles!(w, centers, radii, ylims = (-4, 1), xlims = (-3, 4))
 ```
 
-See also [`trajectory_ss`](@ref) and [`fuzz`](@ref).
+See also [`trajectory_ss`](@ref).
 """
 function batch_ss(args...; kwargs...)
     lins, ssys, resolved_ops = batch_linearize(args...; kwargs...)
@@ -512,7 +512,7 @@ end
 # end
 
 """
-    linsystems, ssys = trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), fuzzer=nothing, verbose = true, kwargs...)
+    linsystems, ssys = trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), kwargs...)
 
 Linearize `sys` around the trajectory `sol` at times `t`. Returns a vector of `StateSpace` objects and the simplified system.
 
@@ -523,56 +523,20 @@ 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`.
-- `fuzzer`: A function that takes an operating point dictionary and returns an array of "fuzzed" operating points. This is useful for adding noise/uncertainty to the operating points along the trajectory. See [`ControlSystemsMTK.fuzz`](@ref) for such a function.
 - `kwargs`: Are sent to the linearization functions (e.g., `loop_openings`).
 - `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), allow_input_derivatives = false, fuzzer = nothing, verbose = true, named = true, kwargs...)
+function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), 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 extract operating points from the solution.
-    # _build_op_from_solution gives differential states + parameters; we supplement
-    # with all unknowns of the linearization system to avoid initialization issues.
-    _extract_base_op(ti) = ModelingToolkit._build_op_from_solution(ModelingToolkit.LinearizationOpPoint(sol, ti))
+    # Use LinearizationOpPoint to let MTK extract operating points from the solution
     tc = collect(t)
-
-    # First pass: get the linearization system's unknowns
-    lin_fun0, ssys = linearization_function(sys, inputs, outputs; warn_initialize_determined=false, kwargs...)
-    lin_unknowns = unknowns(ssys)
-    defs = ModelingToolkit.initial_conditions(sys)
-
-    ops = map(tc) do ti
-        op = _extract_base_op(ti)
-        for x in lin_unknowns
-            haskey(op, x) && continue
-            try
-                op[x] = sol(ti, idxs=x)
-            catch
-                val = get(defs, x, nothing)
-                val !== nothing && (op[x] = val)
-            end
-        end
-        op
-    end
-
-    if fuzzer !== nothing
-        opsv = map(ops) do op
-            fuzzer(op)
-        end
-        ops = reduce(vcat, opsv)
-        tc = repeat(tc, inner = length(ops) ÷ length(tc))
-    end
-
-    lin_fun, ssys = linearization_function(sys, inputs, outputs; op=ops[1], t=tc[1], initialize=false, kwargs...)
-    lins_ops = map(zip(ops, tc)) do (op, ti)
-        linearize(ssys, lin_fun; op, t=ti, allow_input_derivatives)
-    end
-    lins = first.(lins_ops)
-    resolved_ops = last.(lins_ops)
+    op = ModelingToolkit.LinearizationOpPoint(sol, tc)
+    lins, ssys, resolved_ops = linearize(sys, inputs, outputs; op, allow_input_derivatives, kwargs...)
 
     named_linsystems = map(lins) do l
         if named
@@ -583,51 +547,12 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_inp
             ss(l.A, l.B, l.C, l.D)
         end
     end
-    (; linsystems = named_linsystems, ssys, ops, resolved_ops)
+    (; linsystems = named_linsystems, ssys, ops = resolved_ops, resolved_ops)
 end
 
 "_max_100(t) = length(t) > 100 ? range(extrema(t)..., 100) : t"
 _max_100(t) = length(t) > 100 ? range(extrema(t)..., 100) : t
 
-"""
-    fuzz(op, p; N = 10, parameters = true, variables = true)
-
-"Fuzz" an operating point `op::Dict` by changing each non-zero value to an uncertain number with multiplicative uncertainty `p`, represented by `N` samples, i.e., `p = 0.1` means that the value is multiplied by a `N` numbers between 0.9 and 1.1.
-
-`parameters` and `variables` indicate whether to fuzz parameters and state variables, respectively.
-
-This function modifies all variables the same way. For more fine-grained control, load the `MonteCarloMeasurements` package and use the `Particles` type directly, followed by `MonteCarloMeasurements.particle_dict2dict_vec(op)`, i.e., the following makes `uncertain_var` uncertain with a 10% uncertainty:
-```julia
-using MonteCarloMeasurements
-op = ModelingToolkit.defaults(sys)
-op[uncertain_var] = op[uncertain_var] * Particles(10, Uniform(0.9, 1.1))
-ops = MonteCarloMeasurements.particle_dict2dict_vec(op)
-batch_ss(model, inputs, outputs, ops)
-```
-If you have more than one uncertain parameter, it's important to use the same number of particles for all of them (10 in the example above).
-
-To make use of this function in [`trajectory_ss`](@ref), pass something like
-```
-fuzzer = op -> ControlSystemsMTK.fuzz(op, 0.02; N=10)
-```
-to fuzz each operating point 10 times with a 2% uncertainty. The resulting number of operating points will increase by 10x.
-"""
-function fuzz(op, p; N=10, parameters = true, variables = true)
-    op = map(collect(keys(op))) do key
-        par = ModelingToolkit.isparameter(key)
-        val = op[key]
-        par && !parameters && return (key => val)
-        !par && !variables && return (key => val)
-        aval = abs(val)
-        uval = issymbolic(val) ? val : iszero(val) ? 0.0 : Particles(N, MonteCarloMeasurements.Uniform(val-p*aval, val+p*aval))
-        key => uval
-    end |> Dict
-    MonteCarloMeasurements.particle_dict2dict_vec(op)
-end
-
-MonteCarloMeasurements.vecindex(p::Symbolics.BasicSymbolic,i) = p
-issymbolic(x) = x isa Union{Symbolics.Num, Symbolics.BasicSymbolic}
-
 maybe_interp(interpolator, x, t) = allequal(x) ? x[1] : interpolator(x, t)
 
 """
@@ -676,10 +601,10 @@ function GainScheduledStateSpace(systems, vt; interpolator, x = zeros(systems[1]
         description = "Scheduling variable of gain-scheduled statespace system $name",
     ]
     
-    @variables A(t)[1:nx, 1:nx] = systems[1].A
-    @variables B(t)[1:nx, 1:nu] = systems[1].B
-    @variables C(t)[1:ny, 1:nx] = systems[1].C
-    @variables D(t)[1:ny, 1:nu] = systems[1].D
+    @variables A(t)[1:nx, 1:nx] [guess=systems[1].A]
+    @variables B(t)[1:nx, 1:nu] [guess=systems[1].B]
+    @variables C(t)[1:ny, 1:nx] [guess=systems[1].C]
+    @variables D(t)[1:ny, 1:nu] [guess=systems[1].D]
     A,B,C,D = collect.((A,B,C,D))
 
     eqs = [

From 6412266566452383731480ee482473ad5bd3e912 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 13 Apr 2026 08:38:24 +0200
Subject: [PATCH 04/12] test and docs updates

---
 docs/src/batch_linearization.md | 6 +++---
 docs/src/index.md               | 2 +-
 test/test_ODESystem.jl          | 2 ++
 3 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/docs/src/batch_linearization.md b/docs/src/batch_linearization.md
index 4da6462..c0f9b14 100644
--- a/docs/src/batch_linearization.md
+++ b/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")
 ```
diff --git a/docs/src/index.md b/docs/src/index.md
index c1285e3..1222bc6 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -192,7 +192,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)
diff --git a/test/test_ODESystem.jl b/test/test_ODESystem.jl
index 6987218..907e519 100644
--- a/test/test_ODESystem.jl
+++ b/test/test_ODESystem.jl
@@ -162,6 +162,8 @@ using ModelingToolkitStandardLibrary.Mechanical.Rotational
 using ModelingToolkitStandardLibrary.Blocks: Sine
 using ModelingToolkit: connect
 import ModelingToolkitStandardLibrary.Blocks
+Spring = Rotational.Spring
+Damper = Rotational.Damper
 t = Blocks.t
 
 # Parameters

From c0056174bba5742f6d32627812492d7e58828aea Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 13 Apr 2026 09:34:21 +0200
Subject: [PATCH 05/12] rm collect

---
 src/ode_system.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/ode_system.jl b/src/ode_system.jl
index c356195..c5c12e4 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -534,8 +534,7 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_inp
     output_names = reduce(vcat, ap.input.u for ap in vcat(outputs))
 
     # Use LinearizationOpPoint to let MTK extract operating points from the solution
-    tc = collect(t)
-    op = ModelingToolkit.LinearizationOpPoint(sol, tc)
+    op = ModelingToolkit.LinearizationOpPoint(sol, t)
     lins, ssys, resolved_ops = linearize(sys, inputs, outputs; op, allow_input_derivatives, kwargs...)
 
     named_linsystems = map(lins) do l

From d36903a9c8cfb24489d508959e58882d24ee51c9 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Tue, 28 Apr 2026 10:34:25 +0200
Subject: [PATCH 06/12] Use [sources] to pin MTK to PR branch for CI testing

Pins ModelingToolkit and ModelingToolkitBase to the as/linearize-op
PR branch (SciML/ModelingToolkit.jl#4443) so CI can validate the
LinearizationOpPoint integration.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
---
 Project.toml | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/Project.toml b/Project.toml
index 0fbf7b7..244970d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,7 @@ 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"
@@ -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"}

From ab44e46683db8507d41550e790dbe74273e164c7 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Tue, 28 Apr 2026 11:24:49 +0200
Subject: [PATCH 07/12] Add missing initial condition for inertia2 and pin docs
 MTK to PR branch

- test_ODESystem.jl: add inertia2.flange_a.phi => 0.0 to op (resolves
  cyclic guesses error: the two inertias are independent states coupled
  through a spring/damper, both need explicit initial conditions)
- docs/Project.toml: add [sources] to pin MTK to as/linearize-op so the
  docs build can find LinearizationOpPoint

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
---
 docs/Project.toml      | 4 ++++
 test/test_ODESystem.jl | 6 +++++-
 2 files changed, 9 insertions(+), 1 deletion(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 48e803d..0077e29 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -12,3 +12,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"}
diff --git a/test/test_ODESystem.jl b/test/test_ODESystem.jl
index 907e519..a828205 100644
--- a/test/test_ODESystem.jl
+++ b/test/test_ODESystem.jl
@@ -197,7 +197,11 @@ end
 
 model = SystemModel() |> complete
 
-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)
 @test -1000 ∈ lsys.A
 @test -10 ∈ lsys.A

From 5899811ade2b1d81ade6256bb85d6d53f156b265 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Tue, 28 Apr 2026 11:37:16 +0200
Subject: [PATCH 08/12] Add ModelingToolkitBase to docs deps for [sources] pin

Pkg requires packages in [sources] to also be listed in [deps] or
[extras]. Adding ModelingToolkitBase to [deps] satisfies this.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
---
 docs/Project.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/Project.toml b/docs/Project.toml
index 0077e29..436ca87 100644
--- a/docs/Project.toml
+++ b/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"

From b48de94896973ac2277e97db85b8d01277f8103c Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Tue, 28 Apr 2026 12:45:09 +0200
Subject: [PATCH 09/12] Update docs API references to new MTK module locations

Several functions moved from ModelingToolkitStandardLibrary.Blocks to
ModelingToolkit, and reorder_states was renamed to reorder_unknowns.

- api.md: point @docs blocks at ModelingToolkit.{get_sensitivity,
  get_comp_sensitivity, get_looptransfer, reorder_unknowns}
- ode_system.jl: update prose references to reorder_unknowns

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
---
 docs/src/api.md   | 8 ++++----
 src/ode_system.jl | 4 ++--
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/docs/src/api.md b/docs/src/api.md
index c051bb5..74846f1 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -19,13 +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
+ModelingToolkit.reorder_unknowns
 ```
\ No newline at end of file
diff --git a/src/ode_system.jl b/src/ode_system.jl
index c5c12e4..6bfff37 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -345,7 +345,7 @@ The motivation for this function is that ModelingToolkit does not guarantee
 - Which states are selected as states after simplification.
 - The order of the states.
 
-The second problem above, the ordering of the states, can be worked around using `reorder_states`, but the first problem cannot be solved by trivial reordering. This function thus accepts an array of costs for a user-selected state realization, and assembles the correct cost matrix for the state realization selected by MTK. To do this, the funciton needs the linearization (`linear_sys`) as well as the simplified system, both of which are outputs of [`linearize`](@ref).
+The second problem above, the ordering of the states, can be worked around using `reorder_unknowns`, but the first problem cannot be solved by trivial reordering. This function thus accepts an array of costs for a user-selected state realization, and assembles the correct cost matrix for the state realization selected by MTK. To do this, the funciton needs the linearization (`linear_sys`) as well as the simplified system, both of which are outputs of [`linearize`](@ref).
 
 # Arguments:
 - `linear_sys`: Output of [`linearize`](@ref), an object containing a property called `C`. This can be a [`ControlSystemsBase.StateSpace`](@ref) or a `NamedTuple` with a field `C`.
@@ -380,7 +380,7 @@ The motivation for this function is that ModelingToolkit does not guarantee
 - Which states are selected as states after simplification.
 - The order of the states.
 
-The second problem above, the ordering of the states, can be worked around using `reorder_states`, but the first problem cannot be solved by trivial reordering. This function thus accepts an array of costs for a user-selected state realization, and assembles the correct cost matrix for the state realization selected by MTK. To do this, the funciton performs a linearization between inputs and the cost outputs. The linearization is used to determine the matrix entries belonging to states that are not part of the realization chosen by MTK.
+The second problem above, the ordering of the states, can be worked around using `reorder_unknowns`, but the first problem cannot be solved by trivial reordering. This function thus accepts an array of costs for a user-selected state realization, and assembles the correct cost matrix for the state realization selected by MTK. To do this, the funciton performs a linearization between inputs and the cost outputs. The linearization is used to determine the matrix entries belonging to states that are not part of the realization chosen by MTK.
 
 # Arguments:
 - `sys`: The system to be linearized (not simplified).

From cd36409431ea4e68f27bc3aafe1a6f09de1171f4 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Tue, 28 Apr 2026 14:20:36 +0200
Subject: [PATCH 10/12] Add missing inertia2 initial condition in index.md doc
 example

Same fix as in test_ODESystem.jl: the two inertias are independent
states coupled through a spring/damper, both need explicit initial
conditions to avoid the "Cyclic guesses detected" error.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
---
 docs/src/index.md | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/docs/src/index.md b/docs/src/index.md
index 1222bc6..ec985dd 100644
--- a/docs/src/index.md
+++ b/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

From 625da6a592ead2acf9bdd7c3f3d3304283b8f0a9 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Tue, 30 Jun 2026 16:46:31 +0000
Subject: [PATCH 11/12] Hold loop-opening signals at zero in trajectory_ss /
 batch_linearization docs

`trajectory_ss` now accepts an `op` keyword that is forwarded to
`ModelingToolkit.LinearizationOpPoint(sol, t; op)`, supplying operating-point
values that are not available from the solution. This is required by the
updated MTK linearization, which errors if a loop-opening parameter (the signal
turned into a parameter by `loop_openings`) is not assigned a value.

The batch-linearization tutorial now holds each opened signal at 0 via `op`,
and the text is updated to the new semantics: opening at `u` keeps the
scheduling variable connected (controller linearized along the trajectory and
fully isolated), opening only at `y` leaves the scheduling feedback in place
(controller still coupled to the plant), and opening `y`+`v` pins the
scheduling to 0 (a different controller). The previous
`controllersv == controllersu` test no longer holds and is replaced by a check
that opening at `u` yields a lower-order (isolated) controller than opening at
`y`.

Requires the corresponding MTK change (SciML/ModelingToolkit.jl#4679).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 docs/src/batch_linearization.md | 27 ++++++++++++++-------------
 src/ode_system.jl               | 10 ++++++----
 2 files changed, 20 insertions(+), 17 deletions(-)

diff --git a/docs/src/batch_linearization.md b/docs/src/batch_linearization.md
index c0f9b14..4116995 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. 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,13 +218,9 @@ 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(isapprox.(sminreal.(controllersv), sminreal.(controllersu), atol=1e-10))
-```
+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
 ```
@@ -232,6 +228,11 @@ notice how the state of the plant is included in the controller, this is an indi
 ```@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. 
 
 
diff --git a/src/ode_system.jl b/src/ode_system.jl
index 6bfff37..169eeb4 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -523,19 +523,21 @@ 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`).
 - `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), allow_input_derivatives = false, verbose = true, named = true, kwargs...)
+function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), op = Dict(), 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 = ModelingToolkit.LinearizationOpPoint(sol, t)
-    lins, ssys, resolved_ops = linearize(sys, inputs, outputs; op, allow_input_derivatives, kwargs...)
+    # Use LinearizationOpPoint to let MTK extract operating points from the solution.
+    # `op` supplies values not available from `sol` (e.g. loop-opening parameters).
+    oppoint = ModelingToolkit.LinearizationOpPoint(sol, t; op)
+    lins, ssys, resolved_ops = linearize(sys, inputs, outputs; op = oppoint, allow_input_derivatives, kwargs...)
 
     named_linsystems = map(lins) do l
         if named

From 0ff0dec3a92e8b508eb2bd1367c7d386b0f886c1 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Thu, 23 Jul 2026 10:56:02 +0000
Subject: [PATCH 12/12] Point [sources] to MTK master (as/linearize-op merged
 and deleted)

SciML/ModelingToolkit.jl#4443 (which included the trajectory-linearization
perf and loop-opening fixes from #4679) has merged to master, and the
`as/linearize-op` branch was deleted, so the previous pin no longer resolves.
Pin to `master` until an MTK release >= 11.36 (the first version containing
these changes; latest registered is 11.35.2) is available, after which the
`[sources]` entries can be dropped and the compat bumped.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 Project.toml      | 4 ++--
 docs/Project.toml | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/Project.toml b/Project.toml
index 244970d..afebd73 100644
--- a/Project.toml
+++ b/Project.toml
@@ -37,5 +37,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 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"}
+ModelingToolkit = {url = "https://github.com/SciML/ModelingToolkit.jl", rev = "master"}
+ModelingToolkitBase = {url = "https://github.com/SciML/ModelingToolkit.jl", rev = "master", subdir = "lib/ModelingToolkitBase"}
diff --git a/docs/Project.toml b/docs/Project.toml
index 436ca87..399f2b7 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -15,5 +15,5 @@ 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"}
+ModelingToolkit = {url = "https://github.com/SciML/ModelingToolkit.jl", rev = "master"}
+ModelingToolkitBase = {url = "https://github.com/SciML/ModelingToolkit.jl", rev = "master", subdir = "lib/ModelingToolkitBase"}