diff --git a/Project.toml b/Project.toml
index 0ed5462..244970d 100644
--- a/Project.toml
+++ b/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"}
diff --git a/docs/Project.toml b/docs/Project.toml
index 48e803d..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"
@@ -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"}
diff --git a/docs/src/api.md b/docs/src/api.md
index 8197657..74846f1 100644
--- a/docs/src/api.md
+++ b/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
 ```
\ No newline at end of file
diff --git a/docs/src/batch_linearization.md b/docs/src/batch_linearization.md
index 0cded72..4116995 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")
 ```
@@ -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,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(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
 ```
@@ -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/docs/src/index.md b/docs/src/index.md
index c1285e3..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
@@ -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)
diff --git a/src/ode_system.jl b/src/ode_system.jl
index a8fad8f..169eeb4 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
@@ -341,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`.
@@ -376,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).
@@ -488,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...)
@@ -508,162 +512,48 @@ 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.
 
+Operating points are extracted from the solution automatically using `ModelingToolkit.LinearizationOpPoint`.
+
 # 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.
+- `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, fuzzer = nothing, 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)).")
-    # 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)) 
-
-    op_nothing = Dict(unknowns(sys) .=> nothing) # Remove all defaults present in the original system
-    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)
-    end
-    if fuzzer !== nothing
-        opsv = map(ops) do op
-            fuzzer(op)
-        end
-        ops = reduce(vcat, opsv)
-        t = repeat(t, inner = length(ops) ÷ length(t))
-    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
-    end
-    lins = first.(lins_ops)
-    resolved_ops = last.(lins_ops)
+
+    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).
+    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
-            # 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
     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}
-
-"""
-    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)
 
 """
@@ -712,10 +602,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 = [
diff --git a/test/test_ODESystem.jl b/test/test_ODESystem.jl
index 6987218..a828205 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
@@ -195,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