From 49215a093340d70b135d79121f6456dce59d0d50 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Tue, 3 Jun 2025 09:24:38 +0200
Subject: [PATCH 01/20] update to MTKv10

---
 Project.toml           |   4 +-
 src/ode_system.jl      | 142 +++++++----------------------------------
 test/test_ODESystem.jl |  67 +++++++------------
 3 files changed, 46 insertions(+), 167 deletions(-)

diff --git a/Project.toml b/Project.toml
index fe78a39..59bc7e4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ControlSystemsMTK"
 uuid = "687d7614-c7e5-45fc-bfc3-9ee385575c88"
 authors = ["Fredrik Bagge Carlson"]
-version = "2.3.2"
+version = "2.4.0"
 
 [deps]
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
@@ -17,7 +17,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 [compat]
 ControlSystemsBase = "1.0.1"
 DataInterpolations = "3, 4, 5, 6, 7"
-ModelingToolkit = "9.61"
+ModelingToolkit = "10.0"
 ModelingToolkitStandardLibrary = "2"
 MonteCarloMeasurements = "1.1"
 RobustAndOptimalControl = "0.4.14"
diff --git a/src/ode_system.jl b/src/ode_system.jl
index 346b184..9e6c801 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -3,12 +3,12 @@ const AP = Union{Symbol, AnalysisPoint}
 import ModelingToolkitStandardLibrary.Blocks as Blocks
 conn = ModelingToolkit.connect
 t = Blocks.t
-ModelingToolkit.ODESystem(sys::LTISystem; kwargs...) = ODESystem(ss(sys); kwargs...)
+ModelingToolkit.System(sys::LTISystem; kwargs...) = System(ss(sys); kwargs...)
 
 """
-    ModelingToolkit.ODESystem(sys::AbstractStateSpace; name::Symbol, x0 = zeros(sys.nx), x_names, u_names, y_names)
+    ModelingToolkit.System(sys::AbstractStateSpace; name::Symbol, x0 = zeros(sys.nx), x_names, u_names, y_names)
 
-Create an ODESystem from `sys::StateSpace`. 
+Create an System from `sys::StateSpace`. 
 
 # Arguments:
 - `sys`: An instance of `StateSpace` or `NamedStateSpace`.
@@ -19,7 +19,7 @@ The arguments below are automatically set if the system is a `NamedStateSpace`.
 - `u_names`: A vector of symbols with input names. 
 - `y_names`: A vector of symbols with output names. 
 """
-function ModelingToolkit.ODESystem(
+function ModelingToolkit.System(
     sys::AbstractStateSpace;
     name::Symbol,
     x0 = zeros(sys.nx),
@@ -39,112 +39,14 @@ function ModelingToolkit.ODESystem(
         [uc[i].u ~ input.u[i] for i in 1:length(uc)];
         [yc[i].u ~ output.u[i] for i in 1:length(yc)];
     ]
-    extend(ODESystem(eqs, t; name, systems), ssblock)
-end
-
-"""
-    sconnect(input::Function, sys::T; name)
-    sconnect(input::Num,      sys::T; name)
-
-Connect a function `input(t)` to `sys.input`
-
-# Examples:
-```julia
-sconnect(sin, sys)   # Connect a funciton, assumed to be a function of time
-sconnect(sin(t), sys) # Connect a Num
-```
-"""
-function sconnect(
-    input::Union{Function, Num},
-    sys::T;
-    name = Symbol("$(sys.name) with input"),
-) where {T<:ModelingToolkit.AbstractTimeDependentSystem}
-    @named output = Blocks.RealOutput()
-    T(
-        [
-            sys.input.u ~ (input isa Num ? input : input(t))
-            output.u ~ sys.output.u
-        ],
-        t;
-        systems = [sys, output],
-        name,
-    )
-end
-
-"""
-    sconnect(sys1::T, sys2::T; name)
-
-Connect systems in series, equivalent to `sys2*sys1` or `series(sys1, sys2)` in ControlSystems.jl terminology
-"""
-function sconnect(
-    sys1::T,
-    sys2::T;
-    name = Symbol("$(sys1.name)*$(sys2.name)"),
-) where {T<:ModelingToolkit.AbstractTimeDependentSystem}
-    @named output = Blocks.RealOutput() # TODO: missing size
-    @named input = Blocks.RealInput() # TODO: missing size
-    T(
-        [
-            conn(input, sys2.input)
-            conn(output, sys1.output)
-            conn(sys2.output, sys1.input)
-        ],
-        t;
-        name,
-        systems = [sys1, sys2, output, input],
-    )
-end
-
-"""
-    G = ControlSystemsBase.feedback(loopgain::T; name)
-
-Form the feedback-interconnection
-\$G = L/(1+L)\$
-
-The system `G` will be a new system with `input` and `output` connectors.
-"""
-function ControlSystemsBase.feedback(
-    loopgain::T;
-    name = Symbol("feedback $(loopgain.name)"),
-) where {T<:ModelingToolkit.AbstractTimeDependentSystem}
-    add = Blocks.Add(k1 = 1, k2 = -1, name = :feedback)
-    @named input = Blocks.RealInput()
-    @named output = Blocks.RealOutput()
-    T(
-        [
-            input.u ~ add.input1.u
-            output.u ~ loopgain.output.u
-            conn(loopgain.output, add.input2)
-            conn(add.output, loopgain.input)
-        ],
-        t;
-        systems = [input, output, loopgain, add],
-        name,
-    )
-end
-
-function Base.:(*)(s1::T, s2::T) where {T<:ModelingToolkit.AbstractTimeDependentSystem}
-    name = Symbol(string(s1.name) * "_" * string(s2.name))
-    @named input = Blocks.RealInput()
-    @named output = Blocks.RealOutput()
-    eqs = [
-        conn(s1.input, s2.output)
-        output.u ~ s1.output.u
-    ]
-    systems = [output, s1, s2]
-    if any(s.name == :input for s in s2.systems)
-        push!(eqs, input.u ~ s2.input.u)
-        push!(systems, input)
-    end
-    T(eqs, t; systems, name)
+    extend(System(eqs, t; name, systems), ssblock)
 end
 
 
 numeric(x::Num) = x.val
 
-
 function ControlSystemsBase.ss(
-    sys::ModelingToolkit.AbstractTimeDependentSystem,
+    sys::ModelingToolkit.AbstractSystem,
     inputs,
     outputs;
     kwargs...
@@ -154,9 +56,9 @@ end
 
 
 """
-    RobustAndOptimalControl.named_ss(sys::ModelingToolkit.AbstractTimeDependentSystem, inputs, outputs; descriptor=true, kwargs...)
+    RobustAndOptimalControl.named_ss(sys::ModelingToolkit.AbstractSystem, inputs, outputs; descriptor=true, kwargs...)
 
-Convert an `ODESystem` to a `NamedStateSpace` using linearization. `inputs, outputs` are vectors of variables determining the inputs and outputs respectively. See docstring of `ModelingToolkit.linearize` for more info on `kwargs`.
+Convert an `System` to a `NamedStateSpace` using linearization. `inputs, outputs` are vectors of variables determining the inputs and outputs respectively. See docstring of `ModelingToolkit.linearize` for more info on `kwargs`.
 
 If `descriptor = true` (default), this method automatically converts systems that MTK has failed to produce a proper form for into a proper linear statespace system using the method described here:
 https://juliacontrol.github.io/ControlSystemsMTK.jl/dev/#Internals:-Transformation-of-non-proper-models-to-proper-statespace-form
@@ -166,7 +68,7 @@ If `descriptor = false`, the system is instead converted to a statespace realiza
 See also [`ModelingToolkit.linearize`](@ref) which is the lower-level function called internally. The functions [`get_named_sensitivity`](@ref), [`get_named_comp_sensitivity`](@ref), [`get_named_looptransfer`](@ref) similarily provide convenient ways to compute sensitivity functions while retaining signal names in the same way as `named_ss`. The corresponding lower-level functions `get_sensitivity`, `get_comp_sensitivity` and `get_looptransfer` are available in ModelingToolkitStandardLibrary.Blocks and are documented in [MTKstdlib: Linear analysis](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/linear_analysis/).
 """
 function RobustAndOptimalControl.named_ss(
-    sys::ModelingToolkit.AbstractTimeDependentSystem,
+    sys::ModelingToolkit.AbstractSystem,
     inputs,
     outputs;
     descriptor = true,
@@ -177,12 +79,12 @@ function RobustAndOptimalControl.named_ss(
     outputs = vcat(outputs)
 
     inputs = map(inputs) do inp
-        if inp isa ODESystem
+        if inp isa System
             @variables u(t)
             if u ∈ Set(unknowns(inp))
                 inp.u
             else
-                error("Input $(inp.name) is an ODESystem and not a variable")
+                error("Input $(inp.name) is an System and not a variable")
             end
         else
             inp
@@ -191,12 +93,12 @@ function RobustAndOptimalControl.named_ss(
     nu = length(inputs)
 
     outputs = map(outputs) do out
-        if out isa ODESystem
+        if out isa System
             @variables u(t)
             if u ∈ Set(unknowns(out))
                 out.u
             else
-                error("Outut $(out.name) is an ODESystem and not a variable")
+                error("Outut $(out.name) is an System and not a variable")
             end
         else
             out
@@ -286,19 +188,19 @@ get_named_looptransfer
 
 function named_sensitivity_function(
     fun,
-    sys::ModelingToolkit.AbstractTimeDependentSystem,
+    sys::ModelingToolkit.AbstractSystem,
     inputs, args...;
     kwargs...,
 )
 
     inputs = vcat(inputs)
     inputs = map(inputs) do inp
-        if inp isa ODESystem
+        if inp isa System
             @variables u(t)
             if u ∈ Set(unknowns(inp))
                 inp.u
             else
-                error("Input $(inp.name) is an ODESystem and not a variable")
+                error("Input $(inp.name) is an System and not a variable")
             end
         else
             inp
@@ -347,7 +249,7 @@ if isdefined(ModelingToolkit, :get_disturbance_system)
 end
 
 """
-    build_quadratic_cost_matrix(linear_sys, ssys::ODESystem, costs::Vector{Pair})
+    build_quadratic_cost_matrix(linear_sys, ssys::System, costs::Vector{Pair})
 
 For a system that has been linearized, assemble a quadratic cost matrix (for LQR or Kalman filtering) that penalizes states or outputs of simplified system `ssys` according to the vector of pairs `costs`.
 
@@ -362,7 +264,7 @@ The second problem above, the ordering of the states, can be worked around using
 - `ssys`: Output of [`linearize`](@ref).
 - `costs`: A vector of pairs
 """
-function build_quadratic_cost_matrix(matrices::NamedTuple, ssys::ODESystem, costs::AbstractVector{<:Pair})
+function build_quadratic_cost_matrix(matrices::NamedTuple, ssys::System, costs::AbstractVector{<:Pair})
     x = ModelingToolkit.unknowns(ssys)
     y = ModelingToolkit.outputs(ssys)
     # y = getproperty.(ModelingToolkit.observed(ssys), :lhs)
@@ -382,7 +284,7 @@ function build_quadratic_cost_matrix(matrices::NamedTuple, ssys::ODESystem, cost
 end
 
 """
-    build_quadratic_cost_matrix(sys::ODESystem, inputs::Vector, costs::Vector{Pair}; kwargs...)
+    build_quadratic_cost_matrix(sys::System, inputs::Vector, costs::Vector{Pair}; kwargs...)
 
 Assemble a quadratic cost matrix (for LQR or Kalman filtering) that penalizes states or outputs of system `sys` according to the vector of pairs `costs`.
 
@@ -397,7 +299,7 @@ The second problem above, the ordering of the states, can be worked around using
 - `inputs`: A vector of variables that are to be considered controlled inputs for the LQR controller.
 - `costs`: A vector of pairs.
 """
-function build_quadratic_cost_matrix(sys::ODESystem, inputs::AbstractVector, costs::AbstractVector{<:Pair}; kwargs...)
+function build_quadratic_cost_matrix(sys::System, inputs::AbstractVector, costs::AbstractVector{<:Pair}; kwargs...)
     matrices, ssys = ModelingToolkit.linearize(sys, inputs, first.(costs); kwargs...)
     x = ModelingToolkit.unknowns(ssys)
     y = ModelingToolkit.outputs(ssys)
@@ -455,7 +357,7 @@ eqs = [D(x) ~ v
        y ~ x]
 
 
-@named duffing = ODESystem(eqs, t)
+@named duffing = System(eqs, t)
 
 bounds = getbounds(duffing, unknowns(duffing))
 sample_within_bounds((l, u)) = (u - l) * rand() + l
@@ -699,7 +601,7 @@ function GainScheduledStateSpace(systems, vt; interpolator, x = zeros(systems[1]
          for i in 1:nx];
         collect(output.u .~ C * x .+ D * (input.u .- u0) .+ y0)
     ]
-    compose(ODESystem(eqs, t, name = name), [input, output, scheduling_input])
+    compose(System(eqs, t, name = name), [input, output, scheduling_input])
 end
 
 "LPVStateSpace is equivalent to GainScheduledStateSpace, see the docs for GainScheduledStateSpace."
diff --git a/test/test_ODESystem.jl b/test/test_ODESystem.jl
index 3e426f0..b160c00 100644
--- a/test/test_ODESystem.jl
+++ b/test/test_ODESystem.jl
@@ -11,34 +11,12 @@ connect = ModelingToolkit.connect
 P0 = tf(1.0, [1, 1]) |> ss
 C0 = pid(1, 1) * tf(1, [0.01, 1]) |> ss
 
-@named P = ODESystem(P0)
-@test P isa ODESystem
+@named P = System(P0)
+@test P isa System
 # @test length(ModelingToolkit.outputs(P)) == P0.ny
 # @test length(ModelingToolkit.inputs(P)) == P0.nu
 # @named nonlinear_P = sconnect(x->sign(x)*sqrt(abs(x)), P) # apply input-nonlinearity
-@named C        = ODESystem(C0)
-@named loopgain = sconnect(C, P)
-@named ref      = Blocks.Sine(frequency = 1)
-fb0             = feedback(loopgain, name = :fb) * ref
-fb              = structural_simplify(fb0)
-
-# @test length(unknowns(P)) == 3 # 1 + u + y
-# @test length(unknowns(C)) == 4 # 2 + u + y
-
-x0 = Pair[
-    collect(loopgain.P.x) .=> 1.0;
-]
-
-prob = ODEProblem(fb, x0, (0.0, 10.0))
-sol = solve(prob, Rodas5())
-@test Int(sol.retcode) == 1
-isinteractive() && plot(sol)
-
-fb              = structural_simplify(sconnect(sin, P))
-prob = ODEProblem(fb, x0, (0.0, 10.0))
-sol2 = solve(prob, Rodas5())
-@test Int(sol2.retcode) == 1
-isinteractive() && plot(sol2)
+@named C        = System(C0)
 
 Pc = complete(P)
 op = Dict(Pc.input.u => 0.0)
@@ -65,7 +43,7 @@ P1 = ss(Pc, [Pc.input.u], [Pc.output.u]; op)
 
 
 
-# === Go the other way, ODESystem -> StateSpace ================================
+# === Go the other way, System -> StateSpace ================================
 x = unknowns(P) # I haven't figured out a good way to access states, so this is a bit manual and ugly
 @unpack input, output = P
 P02_named = named_ss(P, [input.u], [output.u]; op)
@@ -82,8 +60,8 @@ x = unknowns(C)
 @test ss(C02) == C0
 
 
-## Back again for a complete round trip, test that ODESystem get correct names
-@named P2 = ODESystem(P02_named)
+## Back again for a complete round trip, test that System get correct names
+@named P2 = System(P02_named)
 # @test Set(unknowns(P)) ⊆ Set(unknowns(P2))
 # @test Set(ModelingToolkit.inputs(P)) ⊆ Set(ModelingToolkit.inputs(P2))
 # @test Set(ModelingToolkit.outputs(P)) ⊆ Set(ModelingToolkit.outputs(P2))
@@ -128,13 +106,13 @@ Cp0 = tx * pid(2, 0) * tf(1, [0.01, 1])  #|> ss
 P0i = deepcopy(P0)
 P0i.D .= 1e-8
 
-@named RF = ODESystem(tf(1, [0.001, 1])) # ref filter
-@named RFv = ODESystem(tf(1, [0.001, 1]))
-@named P = ODESystem(P0) # system model
-@named Fv = ODESystem(inv(P0i[2, 1])) # feedforward vel
-@named Fp = ODESystem(inv(P0i[1, 1])) # feedforward pos
-@named Cv = ODESystem(Cv0) # vel controller
-@named Cp = ODESystem(Cp0) # pos controller
+@named RF = System(tf(1, [0.001, 1])) # ref filter
+@named RFv = System(tf(1, [0.001, 1]))
+@named P = System(P0) # system model
+@named Fv = System(inv(P0i[2, 1])) # feedforward vel
+@named Fp = System(inv(P0i[1, 1])) # feedforward pos
+@named Cv = System(Cv0) # vel controller
+@named Cp = System(Cp0) # pos controller
 fb =
     let ref0 = Blocks.Sine(amplitude = 0.2, frequency = 1, name = :r),
         disturbance = Blocks.Step(height = 100, start_time = 10, name = :d)
@@ -142,7 +120,7 @@ fb =
         @named input = Blocks.RealInput()
         @named output = Blocks.RealOutput()
         Dₜ = Differential(t)
-        fb = ODESystem(
+        fb = System(
             [
                 conn(ref0.output, RF.input) # filter position reference
                 expand_derivatives(Dₜ(ref0.output.u) ~ RFv.input.u) # Filter vel reference
@@ -158,16 +136,15 @@ fb =
             name = :feedback,
         )
     end
-simplified_sys = structural_simplify(fb)
+simplified_sys = mtkcompile(fb)
 
 
 x0 = Pair[
     P.x[1] => 0.0
     P.x[3] => 0.0
 ]
-p = Pair[]
 
-prob = ODEProblem(simplified_sys, x0, (0.0, 20.0), p, jac = true)
+prob = ODEProblem(simplified_sys, x0, (0.0, 20.0), jac = true)
 sol = solve(prob, Rodas5(), saveat = 0:0.01:20)
 if isinteractive()
     @show sol.retcode
@@ -211,9 +188,9 @@ function SystemModel(u=nothing; name=:model)
     ]
     if u !== nothing 
         push!(eqs, connect(u.output, :u, torque.tau))
-        return @named model = ODESystem(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper, u])
+        return @named model = System(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper, u])
     end
-    ODESystem(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper], name)
+    System(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper], name)
 end
 
 model = SystemModel() |> complete
@@ -261,14 +238,14 @@ eqs = [connect(P.output, :plant_output, add.input2)
     connect(add.output, C.input)
     connect(C.output, :plant_input, P.input)]
 
-sys_inner = ODESystem(eqs, t, systems = [P, C, add], name = :inner)
+sys_inner = System(eqs, t, systems = [P, C, add], name = :inner)
 
 @named r = Blocks.Constant(k = 1)
 @named F = Blocks.FirstOrder(k = 1, T = 3)
 
 eqs = [connect(r.output, F.input)
     connect(F.output, sys_inner.add.input1)]
-sys_outer = ODESystem(eqs, t, systems = [F, sys_inner, r], name = :outer)
+sys_outer = System(eqs, t, systems = [F, sys_inner, r], name = :outer)
 
 matrices, _ = Blocks.get_sensitivity(sys_outer, [sys_outer.inner.plant_input, sys_outer.inner.plant_output])
 S = ss(matrices...)
@@ -282,7 +259,7 @@ Sn = get_named_sensitivity(sys_outer, [sys_outer.inner.plant_input, sys_outer.in
 
 ## Test connector names
 P = named_ss(ssrand(1,1,1), u=:jörgen, y=:solis)
-@named Pode = ODESystem(P)
+@named Pode = System(P)
 ModelingToolkit.isconnector(Pode.jörgen)
 ModelingToolkit.isconnector(Pode.solis)
 
@@ -313,7 +290,7 @@ eqs = [connect(link1.TX1, cart.flange)
        connect(cart.flange, force.flange)
        connect(link1.TY1, fixed.flange)]
 
-@named model = ODESystem(eqs, t, [], []; systems = [link1, cart, force, fixed])
+@named model = System(eqs, t, [], []; systems = [link1, cart, force, fixed])
 lin_outputs = [cart.s, cart.v, link1.A, link1.dA]
 lin_inputs = [force.f.u]
 

From 3a86b09ea5ae639ac629dbf24ba90b014f543aeb Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Tue, 3 Jun 2025 09:27:55 +0200
Subject: [PATCH 02/20] bump cache version

---
 .github/workflows/CI.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 9b34ab6..31ac343 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -21,7 +21,7 @@ jobs:
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
+      - uses: actions/cache@v2
         env:
           cache-name: cache-artifacts
         with:

From 5a0656a68a401137239b206bf749e1f7c238ede1 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Tue, 3 Jun 2025 09:29:58 +0200
Subject: [PATCH 03/20] up CI script

---
 .github/workflows/CI.yml | 39 ++++++++++++++++++++++-----------------
 1 file changed, 22 insertions(+), 17 deletions(-)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 31ac343..027b603 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -1,40 +1,45 @@
 name: CI
 on:
-  - push
-  - pull_request
+  push:
+    branches:
+      - main
+    tags: ['*']
+  pull_request:
+  workflow_dispatch:
+concurrency:
+  # Skip intermediate builds: always.
+  # Cancel intermediate builds: only if it is a pull request build.
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
 jobs:
   test:
-    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
     runs-on: ${{ matrix.os }}
+    timeout-minutes: 60
+    permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created
+      actions: write
+      contents: read
     strategy:
       fail-fast: false
       matrix:
         version:
-          - '1'
+          - '1.11'
         os:
           - ubuntu-latest
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v2
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-
-            ${{ runner.os }}-
+      - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v4
         with:
-          file: lcov.info
+          files: lcov.info
           token: ${{ secrets.CODECOV_TOKEN }}
+          fail_ci_if_error: false

From aba32639ed1981cfa874d5bd96765dc1b9561403 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Tue, 3 Jun 2025 09:33:04 +0200
Subject: [PATCH 04/20] remove unused export

---
 src/ControlSystemsMTK.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/ControlSystemsMTK.jl b/src/ControlSystemsMTK.jl
index 418cc0f..4b9ad66 100644
--- a/src/ControlSystemsMTK.jl
+++ b/src/ControlSystemsMTK.jl
@@ -29,7 +29,7 @@ using UnPack
 
 # using SymbolicControlSystems
 
-export sconnect, feedback, ODESystem, unknowns, observed, named_ss
+export feedback, ODESystem, unknowns, observed, named_ss
 export batch_ss, trajectory_ss, GainScheduledStateSpace
 export build_quadratic_cost_matrix
 

From e6188ea94649e04e83a9f21b0c7f1d8df4ca13ab Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Tue, 3 Jun 2025 10:06:00 +0200
Subject: [PATCH 05/20] update other files

---
 docs/src/batch_linearization.md |  8 ++++----
 docs/src/index.md               | 14 +++++++-------
 src/ControlSystemsMTK.jl        |  6 +++---
 test/runtests.jl                |  4 ++--
 test/test_batchlin.jl           |  8 ++++----
 5 files changed, 20 insertions(+), 20 deletions(-)

diff --git a/docs/src/batch_linearization.md b/docs/src/batch_linearization.md
index 30feb3c..199ccf5 100644
--- a/docs/src/batch_linearization.md
+++ b/docs/src/batch_linearization.md
@@ -25,7 +25,7 @@ eqs = [D(x) ~ v
        y.u ~ x]
 
 
-@named duffing = ODESystem(eqs, t, systems=[y, u])
+@named duffing = System(eqs, t, systems=[y, u])
 ```
 
 ## Batch linearization
@@ -111,13 +111,13 @@ plot(layout=2)
 
 # Simulate each individual controller
 for C in Cs
-    @named Ci = ODESystem(C)
+    @named Ci = System(C)
     eqs = [
         closed_loop_eqs
         connect(fb.output, Ci.input)
         connect(Ci.output, duffing.u)
     ]
-    @named closed_loop = ODESystem(eqs, t, systems=[duffing, Ci, fb, ref, F])
+    @named closed_loop = System(eqs, t, systems=[duffing, Ci, 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)
     plot!(sol, idxs=[duffing.y.u, duffing.u.u], layout=2, lab="")
@@ -131,7 +131,7 @@ eqs = [
     connect(Cgs.output, duffing.u)
     connect(duffing.y, Cgs.scheduling_input) # Don't forget to connect the scheduling variable!
 ]
-@named closed_loop = ODESystem(eqs, t, systems=[duffing, Cgs, fb, ref, F])
+@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())
 plot!(sol, idxs=[duffing.y.u, duffing.u.u], l=(2, :red), lab="Gain scheduled")
diff --git a/docs/src/index.md b/docs/src/index.md
index 3dfd0ad..a346fa8 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -21,7 +21,7 @@ pkg> add ControlSystemsMTK
 
 
 ## From ControlSystems to ModelingToolkit
-Simply calling `ODESystem(sys)` converts a `StateSpace` object from ControlSystems into the corresponding [`ModelingToolkitStandardLibrary.Blocks.StateSpace`](http://mtkstdlib.sciml.ai/dev/API/blocks/#ModelingToolkitStandardLibrary.Blocks.StateSpace). If `sys` is a [named statespace object](https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/#Named-systems), the names of inputs and outputs will be retained in the `ODESystem` as connectors, that is, if `my_input` is an input variable in the named statespace object, `my_input` will be a connector of type `RealInput` in the resulting ODESystem. Names of state variables are currently ignored.
+Simply calling `System(sys)` converts a `StateSpace` object from ControlSystems into the corresponding [`ModelingToolkitStandardLibrary.Blocks.StateSpace`](http://mtkstdlib.sciml.ai/dev/API/blocks/#ModelingToolkitStandardLibrary.Blocks.StateSpace). If `sys` is a [named statespace object](https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/#Named-systems), the names of inputs and outputs will be retained in the `System` as connectors, that is, if `my_input` is an input variable in the named statespace object, `my_input` will be a connector of type `RealInput` in the resulting System. Names of state variables are currently ignored.
 
 ### Example:
 
@@ -41,7 +41,7 @@ D =
 
 Continuous-time state-space model
 
-julia> @named P = ODESystem(P0)
+julia> @named P = System(P0)
 Model P with 2 equations
 States (3):
   x[1](t) [defaults to 0.0]
@@ -65,15 +65,15 @@ using ControlSystemsMTK, ControlSystemsBase, ModelingToolkit, RobustAndOptimalCo
 P = named_ss(DemoSystems.double_mass_model(outputs = [1,3]), u=:torque, y=[:motor_angle, :load_angle])
 ```
 
-When we convert this system to an ODESystem, we get a system with connectors `P.torque` and `P.motor_angle`, in addition to the standard connectors `P.input` and `P.output`:
+When we convert this system to a `System`, we get a system with connectors `P.torque` and `P.motor_angle`, in addition to the standard connectors `P.input` and `P.output`:
 ```@example CONNECT
-@named P_ode = ODESystem(P)
+@named P_ode = System(P)
 ```
 Here, `P.torque` is equal to `P.input`, so you may choose to connect to either of them. However, since the output is multivariable, the connector `P.output` represents both outputs, while `P.motor_angle` and `P.load_angle` represent the individual scalar outputs.
 
 
 ## From ModelingToolkit to ControlSystems
-An `ODESystem` can be converted to a named statespace object from [RobustAndOptimalControl.jl](https://github.com/JuliaControl/RobustAndOptimalControl.jl) by calling [`named_ss`](@ref)
+A `System` can be converted to a named statespace object from [RobustAndOptimalControl.jl](https://github.com/JuliaControl/RobustAndOptimalControl.jl) by calling [`named_ss`](@ref)
 
 ```julia
 named_ss(ode_sys, inputs, outputs; op)
@@ -167,9 +167,9 @@ function SystemModel(u=nothing; name=:model)
     ]
     if u !== nothing 
         push!(eqs, connect(u.output, :u, torque.tau))
-        return @named model = ODESystem(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper, u])
+        return @named model = System(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper, u])
     end
-    ODESystem(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper], name)
+    System(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper], name)
 end
 
 model = SystemModel() |> complete
diff --git a/src/ControlSystemsMTK.jl b/src/ControlSystemsMTK.jl
index 4b9ad66..1e987e9 100644
--- a/src/ControlSystemsMTK.jl
+++ b/src/ControlSystemsMTK.jl
@@ -3,7 +3,7 @@ using RobustAndOptimalControl: NamedStateSpace
 #=
 Ideas: All connections handled by ModelingToolkit.
 Names: 
-- handled either by named system, or directly in constructor to ODESystem. 
+- handled either by named system, or directly in constructor to System. 
 Functions: 
 - Give me linear system from [u1, u3] to [qm, a]
 If the linearization of a full system produces a named system, one could implement getindex for vectors of names and obtain the desired transfer functions.
@@ -20,7 +20,7 @@ using ModelingToolkit, ControlSystemsBase
 using ControlSystemsBase: ssdata, AbstractStateSpace, Continuous, nstates, noutputs, ninputs
 # using ControlSystemIdentification
 using RobustAndOptimalControl, MonteCarloMeasurements
-import ModelingToolkit: ODESystem, FnType, Symbolics
+import ModelingToolkit: System, FnType, Symbolics
 using ModelingToolkit: unknowns, observed, isdifferential
 using Symbolics
 using Symbolics: jacobian, solve_for
@@ -29,7 +29,7 @@ using UnPack
 
 # using SymbolicControlSystems
 
-export feedback, ODESystem, unknowns, observed, named_ss
+export feedback, System, unknowns, observed, named_ss
 export batch_ss, trajectory_ss, GainScheduledStateSpace
 export build_quadratic_cost_matrix
 
diff --git a/test/runtests.jl b/test/runtests.jl
index a6e5c0a..bea6f84 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,8 +2,8 @@ using ControlSystemsMTK
 using Test
 
 @testset "ControlSystemsMTK.jl" begin
-    @testset "ODESystem" begin
-        @info "Testing ODESystem"        
+    @testset "System" begin
+        @info "Testing System"        
         include("test_ODESystem.jl")
     end
 
diff --git a/test/test_batchlin.jl b/test/test_batchlin.jl
index 11804c4..90dcbec 100644
--- a/test/test_batchlin.jl
+++ b/test/test_batchlin.jl
@@ -19,7 +19,7 @@ eqs = [D(x) ~ v
        y.u ~ x]
 
 
-@named duffing = ODESystem(eqs, t, systems=[y, u])
+@named duffing = System(eqs, t, systems=[y, u])
 
 bounds = getbounds(duffing, unknowns(duffing))
 sample_within_bounds((l, u)) = (u - l) * rand() + l
@@ -39,7 +39,7 @@ Ps, ssys = batch_ss(duffing, inputs, outputs , ops)
 
 using DataInterpolations
 @named Cgs = GainScheduledStateSpace(Ps, xs, interpolator=LinearInterpolation)
-@test Cgs isa ODESystem
+@test Cgs isa System
 # This is tested better in the docs
 
 ## C-code generation
@@ -53,7 +53,7 @@ import ModelingToolkitStandardLibrary.Blocks
 @named fb = Blocks.Add(k2=-1)
 @named ref = Blocks.Square(frequency=1/6, amplitude=0.5, offset=0.5, start_time=1)
 @named F = Blocks.SecondOrder(w=10, d=0.7)
-@named C = ODESystem(pid(1,1,0; state_space=true, Tf=0.01))
+@named C = System(pid(1,1,0; state_space=true, Tf=0.01))
 
 
 closed_loop_eqs = [
@@ -64,7 +64,7 @@ closed_loop_eqs = [
     ModelingToolkit.connect(C.output, duffing.u)
 ]
 
-@named closed_loop = ODESystem(closed_loop_eqs, t, systems=[duffing, C, fb, ref, F])
+@named closed_loop = System(closed_loop_eqs, t, systems=[duffing, C, fb, ref, F])
 
 ssys = structural_simplify(closed_loop)
 prob = ODEProblem(ssys, [F.xd => 0], (0.0, 8.0))

From 7d10603788d6a01eddbcec92ff80e2741a535336 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Wed, 4 Jun 2025 11:07:51 +0200
Subject: [PATCH 06/20] update usage of `get_u0_p`

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

diff --git a/docs/src/index.md b/docs/src/index.md
index a346fa8..c1285e3 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -194,7 +194,8 @@ That's pretty cool, but even nicer is to generate some code for this symbolic sy
 ```@example LINEAIZE_SYMBOLIC
 defs = ModelingToolkit.defaults(simplified_sys)
 defs = merge(Dict(unknowns(model) .=> 0), defs)
-x, pars = ModelingToolkit.get_u0_p(simplified_sys, defs, defs) # Extract the default state and parameter values
+x = ModelingToolkit.get_u0(simplified_sys, defs) # Extract the default state and parameter values
+pars = ModelingToolkit.get_p(simplified_sys, defs, split=false)
 
 fun = Symbolics.build_function(symbolic_sys, unknowns(simplified_sys), ModelingToolkit.parameters(simplified_sys);
     expression = Val{false}, # Generate a compiled function rather than a Julia expression
diff --git a/test/test_ODESystem.jl b/test/test_ODESystem.jl
index b160c00..282a33c 100644
--- a/test/test_ODESystem.jl
+++ b/test/test_ODESystem.jl
@@ -214,7 +214,7 @@ sys = ss((mats...,)[1:4]...)
 
 defs = ModelingToolkit.defaults(ssys)
 defs = merge(Dict(unknowns(model) .=> 0), defs)
-_, p = ModelingToolkit.get_u0_p(ssys, defs, defs)
+p = ModelingToolkit.get_p(ssys, defs, split=false)
 
 
 sympars = ModelingToolkit.parameters(ssys)

From ba9bd65760e37012fb113ccb54d7fe21bc55b7cc Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 9 Jun 2025 09:07:03 +0200
Subject: [PATCH 07/20] more updates

---
 test/test_ODESystem.jl | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/test/test_ODESystem.jl b/test/test_ODESystem.jl
index 282a33c..44be28e 100644
--- a/test/test_ODESystem.jl
+++ b/test/test_ODESystem.jl
@@ -57,7 +57,7 @@ P02 = ss(P02_named)
 # same for controller
 x = unknowns(C)
 @nonamespace C02 = named_ss(C, [C.input], [C.output]; op)
-@test ss(C02) == C0
+@test tf(ss(C02)) ≈ tf(C0)
 
 
 ## Back again for a complete round trip, test that System get correct names
@@ -296,8 +296,12 @@ lin_inputs = [force.f.u]
 
 # => nothing to remove extra defaults
 op = Dict(cart.s => 10, cart.v => 0, link1.A => -pi/2, link1.dA => 0, force.f.u => 0, link1.x1 => nothing, link1.y1 => nothing, link1.x2 => nothing, link1.x_cm => nothing)
+# op = Dict(cart.s => 10, cart.v => 0, link1.A => -pi/2, link1.dA => 0, force.f.u => 0)
+
+guesses = [link1.fy1 => 0.1, cart.f => 0.1]
+
 G = named_ss(model, lin_inputs, lin_outputs; allow_symbolic = true, op,
-    allow_input_derivatives = true, zero_dummy_der = true)
+    allow_input_derivatives = true, zero_dummy_der = true, guesses)
 G = sminreal(G)
 @info "minreal"
 G = minreal(G)

From a7fd9f3ddef27a3fa55ccafebeb01f0dcbfa505b Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 9 Jun 2025 10:48:35 +0200
Subject: [PATCH 08/20] turn off existing defaults

---
 docs/src/batch_linearization.md |  4 ++--
 src/ode_system.jl               | 18 ++++++++++--------
 2 files changed, 12 insertions(+), 10 deletions(-)

diff --git a/docs/src/batch_linearization.md b/docs/src/batch_linearization.md
index 199ccf5..b061a91 100644
--- a/docs/src/batch_linearization.md
+++ b/docs/src/batch_linearization.md
@@ -56,7 +56,7 @@ P = RobustAndOptimalControl.ss2particles(Ps) # convert to a single StateSpace sy
 
 notice how some coefficients are plotted like uncertain numbers `-13.8 ± 4.3`. We can plot such models as well:
 ```@example BATCHLIN
-bodeplot(P, w, legend=:bottomright) # Should look similar to the one above
+bodeplot(P, w, legend=:bottomright, adaptive=false) # Should look similar to the one above
 ```
 
 ## Controller tuning
@@ -154,7 +154,7 @@ The generated code starts by defining the interpolation vector `xs`, this variab
 We can linearize around a trajectory obtained from `solve` using the function [`trajectory_ss`](@ref). We provide it with a vector of time points along the trajectory at which to linearize, and in this case we specify the inputs and outputs to linearize between as analysis points `r` and `y`.
 ```@example BATCHLIN
 timepoints = 0:0.01:8
-Ps2, ssys = trajectory_ss(closed_loop, closed_loop.r, closed_loop.y, sol; t=timepoints, initialize=true, verbose=true)
+Ps2, ssys = trajectory_ss(closed_loop, closed_loop.r, closed_loop.y, sol; t=timepoints, verbose=true)
 bodeplot(Ps2, w, legend=false)
 ```
 
diff --git a/src/ode_system.jl b/src/ode_system.jl
index 9e6c801..7c732d9 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -425,16 +425,18 @@ Linearize `sys` around the trajectory `sol` at times `t`. Returns a vector of `S
 function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_input_derivatives = false, fuzzer = nothing, verbose = 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_initialize_determined = false, kwargs...)
+    lin_fun, ssys = linearization_function(sys, inputs, outputs; warn_empty_op = false, warn_initialize_determined = false, kwargs...)
 
     x = unknowns(ssys)
+    op_nothing = Dict(unknowns(sys) .=> nothing) # Remove all defaults present in the original system
     defs = ModelingToolkit.defaults(sys)
     ops = map(t) do ti
-        Dict(x => robust_sol_getindex(sol, ti, x, defs; verbose) for x in x)
+        opsol = Dict(x => robust_sol_getindex(sol, ti, x, defs; verbose) for x in x)
+        merge(op_nothing, opsol)
     end
+    # @show ops[1]
     if fuzzer !== nothing
         opsv = map(ops) do op
             fuzzer(op)
@@ -442,7 +444,7 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_inp
         ops = reduce(vcat, opsv)
         t = repeat(t, inner = length(ops) ÷ length(t))
     end
-    # lin_fun, ssys = linearization_function(sys, inputs, outputs; op=ops[1], initialize, kwargs...)
+    lin_fun, ssys = linearization_function(sys, inputs, outputs; op=ops[1], kwargs...) # initializealg=ModelingToolkit.SciMLBase.NoInit()
     lins = 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, initialize=false)[1]
@@ -584,10 +586,10 @@ function GainScheduledStateSpace(systems, vt; interpolator, x = zeros(systems[1]
         description = "Scheduling variable of gain-scheduled statespace system $name",
     ]
     
-    @variables A(v)[1:nx, 1:nx] = systems[1].A
-    @variables B(v)[1:nx, 1:nu] = systems[1].B
-    @variables C(v)[1:ny, 1:nx] = systems[1].C
-    @variables D(v)[1:ny, 1:nu] = systems[1].D
+    @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
     A,B,C,D = collect.((A,B,C,D))
 
     eqs = [

From d45c0b4b835a42ffcd96531f66351d9bf6a397d8 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Wed, 18 Jun 2025 08:44:20 +0200
Subject: [PATCH 09/20] handle operating points

---
 Project.toml           |  2 +-
 src/ode_system.jl      | 35 +++++++++++++++++++++++++----------
 test/test_ODESystem.jl |  5 +++--
 3 files changed, 29 insertions(+), 13 deletions(-)

diff --git a/Project.toml b/Project.toml
index 59bc7e4..7a9e495 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,7 +17,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 [compat]
 ControlSystemsBase = "1.0.1"
 DataInterpolations = "3, 4, 5, 6, 7"
-ModelingToolkit = "10.0"
+ModelingToolkit = "10.3"
 ModelingToolkitStandardLibrary = "2"
 MonteCarloMeasurements = "1.1"
 RobustAndOptimalControl = "0.4.14"
diff --git a/src/ode_system.jl b/src/ode_system.jl
index 7c732d9..dfcba71 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -104,7 +104,7 @@ function RobustAndOptimalControl.named_ss(
             out
         end
     end
-    matrices, ssys = ModelingToolkit.linearize(sys, inputs, outputs; kwargs...)
+    matrices, ssys, xpt = ModelingToolkit.linearize(sys, inputs, outputs; kwargs...)
     symstr(x) = Symbol(x isa AnalysisPoint ? x.name : string(x))
     unames = symstr.(inputs)
     if nu > 0 && size(matrices.B, 2) == 2nu
@@ -115,12 +115,19 @@ function RobustAndOptimalControl.named_ss(
     else
         lsys = ss(matrices...)
     end
+    pind = [ModelingToolkit.parameter_index(ssys, i) for i in ModelingToolkit.inputs(ssys)]
+    x0 = xpt.x
+    u0 = [xpt.p[pi] for pi in pind] 
+    xu = (; x = x0, u = u0)
+    extra = Dict(:operating_point => xu)
+
     named_ss(
         lsys;
         x = symstr.(unknowns(ssys)),
         u = unames,
         y = symstr.(outputs),
         name = string(Base.nameof(sys)),
+        extra,
     )
 end
 
@@ -300,7 +307,7 @@ The second problem above, the ordering of the states, can be worked around using
 - `costs`: A vector of pairs.
 """
 function build_quadratic_cost_matrix(sys::System, inputs::AbstractVector, costs::AbstractVector{<:Pair}; kwargs...)
-    matrices, ssys = ModelingToolkit.linearize(sys, inputs, first.(costs); kwargs...)
+    matrices, ssys, extras = ModelingToolkit.linearize(sys, inputs, first.(costs); kwargs...)
     x = ModelingToolkit.unknowns(ssys)
     y = ModelingToolkit.outputs(ssys)
     nx = length(x)
@@ -323,10 +330,12 @@ function batch_linearize(sys, inputs, outputs, ops::AbstractVector{<:AbstractDic
         allow_input_derivatives = false,
         kwargs...)
     lin_fun, ssys = linearization_function(sys, inputs, outputs; op=ops[1], kwargs...)
-    lins = map(ops) do op
+    lins_ops = map(ops) do op
         linearize(ssys, lin_fun; op, t, allow_input_derivatives)
     end
-    lins, ssys
+    lins = first.(lins_ops)
+    resolved_ops = last.(lins_ops)
+    lins, ssys, resolved_ops
 end
 
 """
@@ -404,8 +413,8 @@ nyquistcircles!(w, centers, radii, ylims = (-4, 1), xlims = (-3, 4))
 See also [`trajectory_ss`](@ref) and [`fuzz`](@ref).
 """
 function batch_ss(args...; kwargs...)
-    lins, ssys = batch_linearize(args...; kwargs...)
-    [ss(l...) for l in lins], ssys
+    lins, ssys, resolved_ops = batch_linearize(args...; kwargs...)
+    [ss(l...) for l in lins], ssys, resolved_ops
 end
 
 """
@@ -422,18 +431,22 @@ Linearize `sys` around the trajectory `sol` at times `t`. Returns a vector of `S
 - `verbose`: If `true`, print warnings for variables that are not found in `sol`.
 - `kwargs`: Are sent to the linearization functions.
 """
-function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_input_derivatives = false, fuzzer = nothing, verbose = true, kwargs...)
+function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_input_derivatives = false, fuzzer = nothing, verbose = 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
+    @show perturbation_vars = ModelingToolkit.inputs(ssys)
+    # original_inputs = [ap.input.u for ap in vcat(inputs)] # assuming all inputs are analysis points for now
     op_nothing = Dict(unknowns(sys) .=> nothing) # Remove all defaults present in the original system
     defs = ModelingToolkit.defaults(sys)
     ops = map(t) do ti
         opsol = Dict(x => robust_sol_getindex(sol, ti, x, defs; verbose) for x in x)
+        # 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
     # @show ops[1]
@@ -445,11 +458,13 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_inp
         t = repeat(t, inner = length(ops) ÷ length(t))
     end
     lin_fun, ssys = linearization_function(sys, inputs, outputs; op=ops[1], kwargs...) # initializealg=ModelingToolkit.SciMLBase.NoInit()
-    lins = map(zip(ops, t)) do (op, t)
+    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, initialize=false)[1]
     end
-    (; linsystems = [ss(l...) for l in lins], ssys, ops)
+    lins = first.(lins_ops)
+    resolved_ops = last.(lins_ops)
+    (; linsystems = [ss(l...) for l in lins], ssys, ops, resolved_ops)
 end
 
 "_max_100(t) = length(t) > 100 ? range(extrema(t)..., 100) : t"
diff --git a/test/test_ODESystem.jl b/test/test_ODESystem.jl
index 44be28e..438851f 100644
--- a/test/test_ODESystem.jl
+++ b/test/test_ODESystem.jl
@@ -303,6 +303,7 @@ guesses = [link1.fy1 => 0.1, cart.f => 0.1]
 G = named_ss(model, lin_inputs, lin_outputs; allow_symbolic = true, op,
     allow_input_derivatives = true, zero_dummy_der = true, guesses)
 G = sminreal(G)
+@test 10 ∈ RobustAndOptimalControl.operating_point(G).x
 @info "minreal"
 G = minreal(G)
 @info "poles"
@@ -311,8 +312,8 @@ ps = poles(G)
 @test minimum(abs, ps) < 1e-6
 @test minimum(abs, complex(0, 1.3777260367206716) .- ps) < 1e-10
 
-lsys, syss = linearize(model, lin_inputs, lin_outputs, op = op,
-    allow_input_derivatives = true)
+lsys, syss = linearize(model, lin_inputs, lin_outputs; op = op,
+    allow_input_derivatives = true, guesses)
 lsyss, sysss = ModelingToolkit.linearize_symbolic(model, lin_inputs, lin_outputs;
     allow_input_derivatives = true)
 

From 14cd42b071b3e447962f04630735c6277d809cbe Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Wed, 18 Jun 2025 08:44:33 +0200
Subject: [PATCH 10/20] careful with state-space basis in interpolation

---
 docs/src/batch_linearization.md | 82 ++++++++++++++++++++++++++++-----
 1 file changed, 71 insertions(+), 11 deletions(-)

diff --git a/docs/src/batch_linearization.md b/docs/src/batch_linearization.md
index b061a91..fb87d3d 100644
--- a/docs/src/batch_linearization.md
+++ b/docs/src/batch_linearization.md
@@ -31,7 +31,7 @@ eqs = [D(x) ~ v
 ## Batch linearization
 To perform batch linearization, we create a vector of operating points, and then linearize the model around each of these points. The function [`batch_ss`](@ref) does this for us, and returns a vector of `StateSpace` models, one for each operating point. An operating point is a `Dict` that maps variables in the MTK model to numerical values. In the example below, we simply sample the variables uniformly within their bounds specified when we created the variables (normally, we might want to linearize on stationary points)
 ```@example BATCHLIN
-N = 16 # Number of samples
+N = 5 # Number of samples
 xs = range(getbounds(x)[1], getbounds(x)[2], length=N)
 ops = Dict.(u.u => 0, x .=> xs)
 ```
@@ -61,6 +61,10 @@ bodeplot(P, w, legend=:bottomright, adaptive=false) # Should look similar to the
 
 ## Controller tuning
 Let's also do some controller tuning for the linearized models above. The function `batch_tune` is not really required here, but it shows how we might go about building more sophisticated tools for batch tuning. In this example, we will tune a PID controller using the function [`loopshapingPID`](@ref). Note, this procedure is not limited to tuning a gain-scheduled PID controller, it should work for gain-scheduling of any LTI controller. 
+
+!!! "note" Interpolating between controllers
+    There are multiple ways in which one could interpolate between different controllers. The two most common approaches are to interpolate their outputs, and interpolating their coefficients. When interpolating the coefficients, it is important to ensure that all controllers have the same structure for the interpolation to be meaningful. One may for example interpolate between PID coefficients, or between the coefficients of a state-space model. When interpolating state-space matrices, the systems must all share the same basis, i.e., the state variables must all have the same meaning among the interpolated systems. When converting a transfer function to state-space form, a numerical balancing is performed, this alters the meaning of the state variables and may introduce artificial dynamics to the interpolated system. We thus pass `balance=false` to the function `ss` to avoid this, or pick a form explicitly, e.g., `modal_form`.
+
 ```@example BATCHLIN
 function batch_tune(f, Ps)
     f.(Ps)
@@ -68,7 +72,8 @@ end
 
 Cs = batch_tune(Ps) do P
     C, kp, ki, kd, fig, CF = loopshapingPID(P, 7; Mt=1.2, Tf = 1/100)
-    ss(CF)
+    ss(CF, balance=false)
+    modal_form(ss(CF, balance=true))[1]
 end
 
 P = RobustAndOptimalControl.ss2particles(Ps)
@@ -103,7 +108,7 @@ using DataInterpolations # Required to interpolate between the controllers
 connect    = ModelingToolkit.connect
 
 closed_loop_eqs = [
-    connect(ref.output, F.input)
+    connect(ref.output, :r0, F.input)
     connect(F.output,  :r, fb.input1) # Add an analysis point :r
     connect(duffing.y, :y, fb.input2) # Add an analysis point :y
 ]
@@ -118,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.xd => 0], (0.0, 8.0))
+    prob = ODEProblem(structural_simplify(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
@@ -128,25 +133,25 @@ end
 eqs = [
     closed_loop_eqs
     connect(fb.output, Cgs.input)
-    connect(Cgs.output, duffing.u)
-    connect(duffing.y, Cgs.scheduling_input) # Don't forget to connect the scheduling variable!
+    connect(Cgs.output, :u, duffing.u)
+    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())
+sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8, initializealg=SciMLBase.NoInit(), 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")
 ```
 
-If everything worked as expected, the gain-scheduled controller should perform better than each of the included controllers individually. 
+If everything worked as expected, the gain-scheduled controller should perform reasonably well across the entire scheduling range. 
 
 
 ## C-Code generation
-We can generate C-code to interpolate our controller using the function [`SymbolicControlSystems.print_c_array`](@ref) from [SymbolicControlSystems.jl](https://github.com/JuliaControl/SymbolicControlSystems.jl). If the controller is a standard [`ControlSystemsBase.StateSpace`](@ref) object, a function that filters the input through the controller can be generated by calling [`SymbolicControlSystems.ccode`](@ref). But if the controller is a vector of controllers representing a gain-scheduled controller, a function that creates the interpolated dynamics is written. In the code below, we shorten the vector of controllers to make the generated C-code easier to read by passing `Cs[1:7:end]` and `xs[1:7:end]`
+We can generate C-code to interpolate our controller using the function [`SymbolicControlSystems.print_c_array`](@ref) from [SymbolicControlSystems.jl](https://github.com/JuliaControl/SymbolicControlSystems.jl). If the controller is a standard [`ControlSystemsBase.StateSpace`](@ref) object, a function that filters the input through the controller can be generated by calling [`SymbolicControlSystems.ccode`](@ref). But if the controller is a vector of controllers representing a gain-scheduled controller, a function that creates the interpolated dynamics is written. In the code below, we shorten the vector of controllers to make the generated C-code easier to read by passing `Cs[1:3:end]` and `xs[1:3:end]`
 ```@example BATCHLIN
 using SymbolicControlSystems, ControlSystemsBase
 Cs_disc = c2d.(Cs, 0.05, :tustin) # Discretize the controller before generating code
-code = SymbolicControlSystems.print_c_array(stdout, Cs_disc[1:7:end], xs[1:7:end], "Cs")
+code = SymbolicControlSystems.print_c_array(stdout, Cs_disc[1:3:end], xs[1:3:end], "Cs")
 ```
 The generated code starts by defining the interpolation vector `xs`, this variable is called `Cs_interp_vect` in the generated code. The code then defines all the ``A`` matrices as a 3-dimensional array, followed by a function that performs the interpolation `interpolate_Cs_A`. This function takes the output array as the first argument, a pointer to the 3D array with interpolation matrices, the interpolation vector as well as the interpolation variable `t`, in this document called ``v``. The same code is then repeated for the matrices ``B,C,D`` as well if they require interpolation (if they are all the same, no interpolation code is written). 
 
@@ -154,13 +159,68 @@ The generated code starts by defining the interpolation vector `xs`, this variab
 We can linearize around a trajectory obtained from `solve` using the function [`trajectory_ss`](@ref). We provide it with a vector of time points along the trajectory at which to linearize, and in this case we specify the inputs and outputs to linearize between as analysis points `r` and `y`.
 ```@example BATCHLIN
 timepoints = 0:0.01:8
-Ps2, ssys = trajectory_ss(closed_loop, closed_loop.r, closed_loop.y, sol; t=timepoints, verbose=true)
+Ps2, ssys, ops2, resolved_ops = trajectory_ss(closed_loop, closed_loop.r, closed_loop.y, sol; t=timepoints, verbose=true);
 bodeplot(Ps2, w, legend=false)
 ```
 
 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.
 
 
+Notice how the peak of the transfer function changes along the trajectory. 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:
+- 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]`
+
+We will explore these options below, starting with the first option, breaking the connection `y`:
+```@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 = trajectory_ss(closed_loop, closed_loop.r, closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.y]);
+
+closed_loopsy = feedback.(plants .* controllersy)
+bodeplot(closed_loopsy, w; title="Loop open at y", kwargs...)
+```
+When we open the loop at `u`, we get a slightly different result:
+```@example BATCHLIN
+controllersu, ssu = trajectory_ss(closed_loop, closed_loop.r, closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.u]);
+
+closed_loopsu = feedback.(plants .* controllersu)
+bodeplot(closed_loopsu, w; title="Loop open at u", kwargs...)
+```
+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`:
+```@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]);
+
+closed_loopsv = feedback.(plants .* controllersv)
+bodeplot(closed_loopsv, w; title="Loop open at v and y", kwargs...)
+```
+
+We have thus far treated the controller as a SISO system, but we could also view it as a system with two inputs, the measurement feedback and the scheduling feedback. For completeness, we show below how to derive the corresponding MIMO systems
+
+```@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]);
+
+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...)
+```
+
+
+
+Why the discrepancy? We can understand the difference by comparing the Bode plots of the controllers. 
+
+```@example BATCHLIN
+plot(
+    bodeplot(Cs,           w, legend=false, plotphase=false, title="Designed controllers"),
+    bodeplot(controllersy, w, legend=false, plotphase=false, title="Loop open at y"),
+    bodeplot(controllersu, w, legend=false, plotphase=false, title="Loop open at u"),
+    bodeplot(controllersv, w, legend=false, plotphase=false, title="Loop open at v and y"),
+)
+```
+if we open at both `y` and `v`, we only get the controller designed for the origin. If 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).
+
+However, if we only open at `y` we get controller linearizations that contains the closed loop through the scheduling connection `v`. The easiest way to ensure that the controller is properly disconnected from the plant while taking into account the different scheduling along the trajectory is thus to break at `u`.
+
 ## Summary
 We have seen how to
 - Perform linearization of a nonlinear ModelingToolkit model in multiple different operating points

From 01531475e3a5fd67404d22fe1083655482ce4de9 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Thu, 19 Jun 2025 10:52:33 +0200
Subject: [PATCH 11/20] simplify implementation of `named_sensitivity_function`

---
 docs/src/batch_linearization.md |  2 +-
 src/ode_system.jl               | 24 ++++++------------------
 2 files changed, 7 insertions(+), 19 deletions(-)

diff --git a/docs/src/batch_linearization.md b/docs/src/batch_linearization.md
index fb87d3d..7acc34a 100644
--- a/docs/src/batch_linearization.md
+++ b/docs/src/batch_linearization.md
@@ -175,7 +175,7 @@ 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 = 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]);
 
 closed_loopsy = feedback.(plants .* controllersy)
 bodeplot(closed_loopsy, w; title="Loop open at y", kwargs...)
diff --git a/src/ode_system.jl b/src/ode_system.jl
index dfcba71..6a61e3a 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -197,12 +197,13 @@ function named_sensitivity_function(
     fun,
     sys::ModelingToolkit.AbstractSystem,
     inputs, args...;
+    descriptor = true,
     kwargs...,
 )
 
     inputs = vcat(inputs)
     inputs = map(inputs) do inp
-        if inp isa System
+    if inp isa System
             @variables u(t)
             if u ∈ Set(unknowns(inp))
                 inp.u
@@ -219,23 +220,10 @@ function named_sensitivity_function(
     unames = symstr.(inputs)
     fm(x) = convert(Matrix{Float64}, x)
     if nu > 0 && size(matrices.B, 2) == 2nu
-        nx = size(matrices.A, 1)
-         # This indicates that input derivatives are present
-        duinds = findall(any(!iszero, eachcol(matrices.B[:, nu+1:end])))
-        B̄ = matrices.B[:, duinds .+ nu]
-        ndu = length(duinds)
-        B = matrices.B[:, 1:nu]
-        Iu = duinds .== (1:nu)'
-        E = [I(nx) -B̄; zeros(ndu, nx+ndu)]
-
-        Ae = cat(matrices.A, -I(ndu), dims=(1,2))
-        Be = [B; Iu]
-        Ce = [fm(matrices.C) zeros(ny, ndu)]
-        De = fm(matrices.D[:, 1:nu])
-        dsys = dss(Ae, E, Be, Ce, De)
-        lsys = ss(RobustAndOptimalControl.DescriptorSystems.dss2ss(dsys)[1])
-        # unames = [unames; Symbol.("der_" .* string.(unames))]
-        # sys = ss(matrices...)
+        # This indicates that input derivatives are present
+        duinds = findall(any(!iszero, eachcol(matrices.B[:, nu+1:end]))) .+ nu
+        u2du = (1:nu) .=> duinds # This maps inputs to their derivatives
+        lsys = causal_simplification(matrices, u2du; descriptor)
     else
         lsys = ss(matrices...)
     end

From 94229352c06c7620ce375363bbf63547d0feaf2a Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Thu, 3 Jul 2025 13:04:57 +0200
Subject: [PATCH 12/20] add specification of op for original inputs

does not work, since the original inputs have the wrong namespace
---
 src/ode_system.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/ode_system.jl b/src/ode_system.jl
index 6a61e3a..0a0c286 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -419,7 +419,7 @@ Linearize `sys` around the trajectory `sol` at times `t`. Returns a vector of `S
 - `verbose`: If `true`, print warnings for variables that are not found in `sol`.
 - `kwargs`: Are sent to the linearization functions.
 """
-function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_input_derivatives = false, fuzzer = nothing, verbose = true,kwargs...)
+function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_input_derivatives = false, fuzzer = nothing, verbose = 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.
@@ -428,12 +428,12 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_inp
     x = unknowns(ssys)
 
     # TODO: The value of the output (or input) of the input analysis points should be mapped to the perturbation vars
-    @show perturbation_vars = ModelingToolkit.inputs(ssys)
-    # original_inputs = [ap.input.u for ap in vcat(inputs)] # assuming all inputs are analysis points for now
+    perturbation_vars = ModelingToolkit.inputs(ssys)
+    original_inputs = [ap.input.u for ap in vcat(inputs)] # assuming all inputs are analysis points for now
     op_nothing = Dict(unknowns(sys) .=> nothing) # Remove all defaults present in the original system
     defs = ModelingToolkit.defaults(sys)
     ops = map(t) do ti
-        opsol = Dict(x => robust_sol_getindex(sol, ti, x, defs; verbose) for x in x)
+        opsol = Dict(x => robust_sol_getindex(sol, ti, x, defs; verbose) for x in [x; original_inputs])
         # 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

From af8a3dd2e1cd020ec8eb81ad8c12d5c2ca0a5882 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Thu, 3 Jul 2025 13:35:02 +0200
Subject: [PATCH 13/20] unnamespace

---
 src/ode_system.jl | 17 ++++++++++++++---
 1 file changed, 14 insertions(+), 3 deletions(-)

diff --git a/src/ode_system.jl b/src/ode_system.jl
index 0a0c286..53e7b93 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -405,6 +405,17 @@ function batch_ss(args...; kwargs...)
     [ss(l...) for l in lins], ssys, resolved_ops
 end
 
+function unnamespace(ap)
+    # ap_name = ap.name
+    # new_name = join(ModelingToolkit.namespace_hierarchy(ap_name)[2:end], Symbolics.NAMESPACE_SEPARATOR)
+    # ap2 = ModelingToolkit.@set ap.name = Symbol(new_name)
+    # ap2
+
+    ap_name = ModelingToolkit.SymbolicIndexingInterface.getname(ap.input.u)
+    new_name = join(ModelingToolkit.namespace_hierarchy(ap_name)[2:end], Symbolics.NAMESPACE_SEPARATOR)
+    new_var = Symbolics.rename(ap.input.u, Symbol(new_name))
+end
+
 """
     linsystems, ssys = trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), fuzzer=nothing, verbose = true, kwargs...)
 
@@ -429,7 +440,7 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_inp
 
     # 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 = [ap.input.u for ap in vcat(inputs)] # assuming all inputs are analysis points for now
+    original_inputs = [unnamespace(ap) for ap in vcat(inputs)] # assuming all inputs are analysis points for now
     op_nothing = Dict(unknowns(sys) .=> nothing) # Remove all defaults present in the original system
     defs = ModelingToolkit.defaults(sys)
     ops = map(t) do ti
@@ -437,7 +448,6 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_inp
         # 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
-    # @show ops[1]
     if fuzzer !== nothing
         opsv = map(ops) do op
             fuzzer(op)
@@ -445,8 +455,9 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_inp
         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...) # initializealg=ModelingToolkit.SciMLBase.NoInit()
+    lin_fun, ssys = linearization_function(sys, inputs, outputs; op=ops[1], initialization_abstol=1e-2, initialization_reltol=1e-2, kwargs...) # initializealg=ModelingToolkit.SciMLBase.NoInit()
     lins_ops = map(zip(ops, t)) do (op, t)
+        @show t
         linearize(ssys, lin_fun; op, t, allow_input_derivatives)
         # linearize(sys, inputs, outputs; op, t, allow_input_derivatives, initialize=false)[1]
     end

From 929b9f5a1f54cb1ab2419c49101d086963d1c14a Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 4 Jul 2025 14:16:45 +0200
Subject: [PATCH 14/20] finally working

---
 docs/src/batch_linearization.md | 47 ++++++++++++++++++++++-----------
 src/ode_system.jl               | 31 +++++++++++-----------
 2 files changed, 47 insertions(+), 31 deletions(-)

diff --git a/docs/src/batch_linearization.md b/docs/src/batch_linearization.md
index 7acc34a..9924276 100644
--- a/docs/src/batch_linearization.md
+++ b/docs/src/batch_linearization.md
@@ -115,18 +115,18 @@ closed_loop_eqs = [
 plot(layout=2)
 
 # Simulate each individual controller
-for C in Cs
-    @named Ci = System(C)
-    eqs = [
-        closed_loop_eqs
-        connect(fb.output, Ci.input)
-        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))
-    sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8)
-    plot!(sol, idxs=[duffing.y.u, duffing.u.u], layout=2, lab="")
-end
+# for C in Cs
+#     @named Ci = System(C)
+#     eqs = [
+#         closed_loop_eqs
+#         connect(fb.output, Ci.input)
+#         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))
+#     sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8)
+#     plot!(sol, idxs=[duffing.y.u, duffing.u.u], layout=2, lab="")
+# end
 
 # Simulate gain-scheduled controller
 @named Cgs = GainScheduledStateSpace(Cs, xs, interpolator=LinearInterpolation)
@@ -162,11 +162,12 @@ timepoints = 0:0.01:8
 Ps2, ssys, ops2, resolved_ops = trajectory_ss(closed_loop, closed_loop.r, closed_loop.y, sol; t=timepoints, verbose=true);
 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.
 
 
-Notice how the peak of the transfer function changes along the trajectory. 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. 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:
 - 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]`
@@ -207,7 +208,7 @@ bodeplot(closed_loops_mimo, w; title="Loop open at MIMO", kwargs...)
 
 
 
-Why the discrepancy? We can understand the difference by comparing the Bode plots of the controllers. 
+Why are the results different depending on where we open the loop? We can understand the difference by comparing the Bode plots of the controllers. 
 
 ```@example BATCHLIN
 plot(
@@ -217,9 +218,23 @@ plot(
     bodeplot(controllersv, w, legend=false, plotphase=false, title="Loop open at v and y"),
 )
 ```
-if we open at both `y` and `v`, we only get the controller designed for the origin. If 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).
+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).
+
+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
+```@example BATCHLIN
+Cy = named_ss(controllersy[end], x=Symbol.(unknowns(ssy)))
+sminreal(Cy).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
+Cu = named_ss(controllersu[end], x=Symbol.(unknowns(ssu)))
+sminreal(Cu).x
+```
+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. 
+
+
 
-However, if we only open at `y` we get controller linearizations that contains the closed loop through the scheduling connection `v`. The easiest way to ensure that the controller is properly disconnected from the plant while taking into account the different scheduling along the trajectory is thus to break at `u`.
+The easiest way to ensure that the controller is properly disconnected from the plant while taking into account the different scheduling along the trajectory is thus to break at `u`.
 
 ## Summary
 We have seen how to
diff --git a/src/ode_system.jl b/src/ode_system.jl
index 53e7b93..b84e80a 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -405,16 +405,13 @@ function batch_ss(args...; kwargs...)
     [ss(l...) for l in lins], ssys, resolved_ops
 end
 
-function unnamespace(ap)
-    # ap_name = ap.name
-    # new_name = join(ModelingToolkit.namespace_hierarchy(ap_name)[2:end], Symbolics.NAMESPACE_SEPARATOR)
-    # ap2 = ModelingToolkit.@set ap.name = Symbol(new_name)
-    # ap2
-
-    ap_name = ModelingToolkit.SymbolicIndexingInterface.getname(ap.input.u)
-    new_name = join(ModelingToolkit.namespace_hierarchy(ap_name)[2:end], Symbolics.NAMESPACE_SEPARATOR)
-    new_var = Symbolics.rename(ap.input.u, Symbol(new_name))
-end
+# function unnamespace(ap)
+#     map(ap.outputs) do out
+#         ap_name = ModelingToolkit.SymbolicIndexingInterface.getname(out.u) 
+#         new_name = join(ModelingToolkit.namespace_hierarchy(ap_name)[2:end], Symbolics.NAMESPACE_SEPARATOR)
+#         Symbolics.rename(ap.input.u, Symbol(new_name))
+#     end
+# end
 
 """
     linsystems, ssys = trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), fuzzer=nothing, verbose = true, kwargs...)
@@ -440,11 +437,11 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_inp
 
     # 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 = [unnamespace(ap) for ap in vcat(inputs)] # assuming all inputs are analysis points for now
+    # @show original_inputs = reduce(vcat, unnamespace(ap) for ap in vcat(inputs)) # assuming all inputs are analysis points for now
     op_nothing = Dict(unknowns(sys) .=> nothing) # Remove all defaults present in the original system
     defs = ModelingToolkit.defaults(sys)
     ops = map(t) do ti
-        opsol = Dict(x => robust_sol_getindex(sol, ti, x, defs; verbose) for x in [x; original_inputs])
+        opsol = Dict(x => robust_sol_getindex(sol, ti, x, defs; verbose) for x in x)
         # 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
@@ -455,11 +452,15 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_inp
         ops = reduce(vcat, opsv)
         t = repeat(t, inner = length(ops) ÷ length(t))
     end
-    lin_fun, ssys = linearization_function(sys, inputs, outputs; op=ops[1], initialization_abstol=1e-2, initialization_reltol=1e-2, kwargs...) # initializealg=ModelingToolkit.SciMLBase.NoInit()
+    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)
-        @show t
         linearize(ssys, lin_fun; op, t, allow_input_derivatives)
-        # linearize(sys, inputs, outputs; op, t, allow_input_derivatives, initialize=false)[1]
+        # linearize(sys, inputs, outputs; op, t, allow_input_derivatives) # useful for debugging
     end
     lins = first.(lins_ops)
     resolved_ops = last.(lins_ops)

From 9c1ce5cc91dd48e48906f3e1b209f96fa9a37465 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 4 Jul 2025 14:20:20 +0200
Subject: [PATCH 15/20] add a test

---
 docs/src/batch_linearization.md | 29 ++++++++++++++++-------------
 1 file changed, 16 insertions(+), 13 deletions(-)

diff --git a/docs/src/batch_linearization.md b/docs/src/batch_linearization.md
index 9924276..0cb0676 100644
--- a/docs/src/batch_linearization.md
+++ b/docs/src/batch_linearization.md
@@ -115,18 +115,18 @@ closed_loop_eqs = [
 plot(layout=2)
 
 # Simulate each individual controller
-# for C in Cs
-#     @named Ci = System(C)
-#     eqs = [
-#         closed_loop_eqs
-#         connect(fb.output, Ci.input)
-#         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))
-#     sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8)
-#     plot!(sol, idxs=[duffing.y.u, duffing.u.u], layout=2, lab="")
-# end
+for C in Cs
+    @named Ci = System(C)
+    eqs = [
+        closed_loop_eqs
+        connect(fb.output, Ci.input)
+        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))
+    sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8)
+    plot!(sol, idxs=[duffing.y.u, duffing.u.u], layout=2, lab="")
+end
 
 # Simulate gain-scheduled controller
 @named Cgs = GainScheduledStateSpace(Cs, xs, interpolator=LinearInterpolation)
@@ -219,6 +219,10 @@ 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))
+```
 
 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
 ```@example BATCHLIN
@@ -248,6 +252,5 @@ Batch linearization in multiple different operating points is an intuitive way t
 
 
 ```@example BATCHLIN
-using Test
 @test sol(6.99, idxs=closed_loop.duffing.y.u) ≈ 0.0 atol=0.01
 ```
\ No newline at end of file

From a9fa17ff4730ed9b2cdd091fa2c7656a500f5c43 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 4 Jul 2025 14:34:23 +0200
Subject: [PATCH 16/20] return named systems

---
 docs/src/batch_linearization.md |  6 ++----
 src/ode_system.jl               | 15 ++++++++++++---
 2 files changed, 14 insertions(+), 7 deletions(-)

diff --git a/docs/src/batch_linearization.md b/docs/src/batch_linearization.md
index 0cb0676..0cded72 100644
--- a/docs/src/batch_linearization.md
+++ b/docs/src/batch_linearization.md
@@ -226,13 +226,11 @@ using Test
 
 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
 ```@example BATCHLIN
-Cy = named_ss(controllersy[end], x=Symbol.(unknowns(ssy)))
-sminreal(Cy).x
+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
-Cu = named_ss(controllersu[end], x=Symbol.(unknowns(ssu)))
-sminreal(Cu).x
+sminreal(controllersu[end]).x
 ```
 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 b84e80a..574ac43 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -54,6 +54,8 @@ function ControlSystemsBase.ss(
     named_ss(sys, inputs, outputs; kwargs...).sys # just discard the names
 end
 
+symstr(x) = Symbol(x isa AnalysisPoint ? x.name : string(x))
+
 
 """
     RobustAndOptimalControl.named_ss(sys::ModelingToolkit.AbstractSystem, inputs, outputs; descriptor=true, kwargs...)
@@ -105,7 +107,6 @@ function RobustAndOptimalControl.named_ss(
         end
     end
     matrices, ssys, xpt = ModelingToolkit.linearize(sys, inputs, outputs; kwargs...)
-    symstr(x) = Symbol(x isa AnalysisPoint ? x.name : string(x))
     unames = symstr.(inputs)
     if nu > 0 && size(matrices.B, 2) == 2nu
         # This indicates that input derivatives are present
@@ -437,7 +438,11 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_inp
 
     # 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)
-    # @show original_inputs = reduce(vcat, unnamespace(ap) for ap in vcat(inputs)) # assuming all inputs are analysis points for now
+    # 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.defaults(sys)
     ops = map(t) do ti
@@ -464,7 +469,11 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_inp
     end
     lins = first.(lins_ops)
     resolved_ops = last.(lins_ops)
-    (; linsystems = [ss(l...) for l in lins], ssys, ops, resolved_ops)
+    named_linsystems = map(lins) do l
+        # Convert to a NamedStateSpace with the same names as the original system
+        named_ss(ss(l.A, l.B, l.C, l.D); name = string(Base.nameof(sys)), x = symstr.(unknowns(ssys)), u = symstr.(input_names), y = symstr.(output_names))
+    end
+    (; linsystems = named_linsystems, ssys, ops, resolved_ops)
 end
 
 "_max_100(t) = length(t) > 100 ? range(extrema(t)..., 100) : t"

From b2b3f670c999c0b96998926d3505f71211078ac4 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 4 Jul 2025 14:39:55 +0200
Subject: [PATCH 17/20] add names to batch_ss as well

---
 src/ode_system.jl | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/src/ode_system.jl b/src/ode_system.jl
index 574ac43..130d463 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -346,7 +346,7 @@ unsafe_comparisons(true)
 @variables x(t)=0 [bounds = (-2, 2)]
 @variables v(t)=0
 @variables u(t)=0
-@variables y(t)=0
+@variables y(t)
 
 D = Differential(t)
 
@@ -360,6 +360,7 @@ eqs = [D(x) ~ v
 bounds = getbounds(duffing, unknowns(duffing))
 sample_within_bounds((l, u)) = (u - l) * rand() + l
 # Create a vector of operating points
+N = 10
 ops = map(1:N) do i
     op = Dict(x => sample_within_bounds(bounds[x]) for x in keys(bounds) if isfinite(bounds[x][1]))
 end
@@ -403,7 +404,11 @@ See also [`trajectory_ss`](@ref) and [`fuzz`](@ref).
 """
 function batch_ss(args...; kwargs...)
     lins, ssys, resolved_ops = batch_linearize(args...; kwargs...)
-    [ss(l...) for l in lins], ssys, resolved_ops
+    named_linsystems = map(lins) do l
+        # Convert to a NamedStateSpace with the same names as the original system
+        named_ss(ss(l.A, l.B, l.C, l.D); name = string(Base.nameof(ssys)), x = symstr.(unknowns(ssys)))
+    end
+    named_linsystems, ssys, resolved_ops
 end
 
 # function unnamespace(ap)

From f82ba4cca2d55979f4659e2de94ca8f20243e025 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 4 Jul 2025 18:53:24 +0200
Subject: [PATCH 18/20] getfield -> getproperty

---
 src/ode_system.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/ode_system.jl b/src/ode_system.jl
index 130d463..1d8116d 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -600,10 +600,10 @@ function GainScheduledStateSpace(systems, vt; interpolator, x = zeros(systems[1]
     s1 = first(systems)
     (; nx, nu, ny) = s1
 
-    Aint = [maybe_interp(interpolator, getindex.(getfield.(systems, :A), i, j), vt) for i = 1:nx, j = 1:nx]
-    Bint = [maybe_interp(interpolator, getindex.(getfield.(systems, :B), i, j), vt) for i = 1:nx, j = 1:nu]
-    Cint = [maybe_interp(interpolator, getindex.(getfield.(systems, :C), i, j), vt) for i = 1:ny, j = 1:nx]
-    Dint = [maybe_interp(interpolator, getindex.(getfield.(systems, :D), i, j), vt) for i = 1:ny, j = 1:nu]
+    Aint = [maybe_interp(interpolator, getindex.(getproperty.(systems, :A), i, j), vt) for i = 1:nx, j = 1:nx]
+    Bint = [maybe_interp(interpolator, getindex.(getproperty.(systems, :B), i, j), vt) for i = 1:nx, j = 1:nu]
+    Cint = [maybe_interp(interpolator, getindex.(getproperty.(systems, :C), i, j), vt) for i = 1:ny, j = 1:nx]
+    Dint = [maybe_interp(interpolator, getindex.(getproperty.(systems, :D), i, j), vt) for i = 1:ny, j = 1:nu]
 
     @named input = Blocks.RealInput(nin = nu)
     @named scheduling_input = Blocks.RealInput()

From e20ba3e763c5b4952937969e205eed01a2c1398b Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 4 Jul 2025 19:20:18 +0200
Subject: [PATCH 19/20] .sys

---
 test/test_batchlin.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/test_batchlin.jl b/test/test_batchlin.jl
index 90dcbec..e7cc035 100644
--- a/test/test_batchlin.jl
+++ b/test/test_batchlin.jl
@@ -32,8 +32,8 @@ inputs, outputs = [u.u], [y.u]
 Ps, ssys = batch_ss(duffing, inputs, outputs , ops)
 @test length(Ps) == N
 
-@test Ps[1] == ss(linearize(duffing, inputs, outputs; op=ops[1])[1]...)
-@test Ps[end] == ss(linearize(duffing, inputs, outputs; op=ops[end])[1]...)
+@test Ps[1].sys == ss(linearize(duffing, inputs, outputs; op=ops[1])[1]...)
+@test Ps[end].sys == ss(linearize(duffing, inputs, outputs; op=ops[end])[1]...)
 
 ##
 

From 13de3856e884566184d7f5ce9f9a620922d10650 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 4 Jul 2025 19:34:32 +0200
Subject: [PATCH 20/20] handle non-unique names

---
 src/ode_system.jl | 14 +++++++++++---
 1 file changed, 11 insertions(+), 3 deletions(-)

diff --git a/src/ode_system.jl b/src/ode_system.jl
index 3e59ad2..29b6118 100644
--- a/src/ode_system.jl
+++ b/src/ode_system.jl
@@ -439,8 +439,9 @@ Linearize `sys` around the trajectory `sol` at times `t`. Returns a vector of `S
 - `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.
+- `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, kwargs...)
+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.
@@ -482,8 +483,15 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_inp
     lins = first.(lins_ops)
     resolved_ops = last.(lins_ops)
     named_linsystems = map(lins) do l
-        # Convert to a NamedStateSpace with the same names as the original system
-        named_ss(ss(l.A, l.B, l.C, l.D); name = string(Base.nameof(sys)), x = symstr.(unknowns(ssys)), u = symstr.(input_names), y = symstr.(output_names))
+        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)
 end