diff --git a/docs/src/basics/MTKModel_Connector.md b/docs/src/basics/MTKModel_Connector.md
index 8cc106b05d..ac285d0253 100644
--- a/docs/src/basics/MTKModel_Connector.md
+++ b/docs/src/basics/MTKModel_Connector.md
@@ -130,6 +130,7 @@ end
 `@connector`s accepts begin blocks of `@components`, `@equations`, `@extend`, `@parameters`, `@structural_parameters`, `@variables`. These keywords mean the same as described above for `@mtkmodel`.
 
 !!! note
+    
     For more examples of usage, checkout [ModelingToolkitStandardLibrary.jl](https://github.com/SciML/ModelingToolkitStandardLibrary.jl/)
 
 ### What's a `structure` dictionary?
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 1e4531a33b..ab68cc9e7e 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -229,7 +229,8 @@ for prop in [:eqs
     :metadata
     :gui_metadata
     :discrete_subsystems
-    :unknown_states]
+    :unknown_states
+    :split_idxs]
     fname1 = Symbol(:get_, prop)
     fname2 = Symbol(:has_, prop)
     @eval begin
@@ -1273,14 +1274,34 @@ See also [`linearize`](@ref) which provides a higher-level interface.
 function linearization_function(sys::AbstractSystem, inputs,
     outputs; simplify = false,
     initialize = true,
+    op = Dict(),
+    p = DiffEqBase.NullParameters(),
+    zero_dummy_der = false,
     kwargs...)
-    sys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs; simplify,
+    ssys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs;
+        simplify,
         kwargs...)
+    if zero_dummy_der
+        dummyder = setdiff(states(ssys), states(sys))
+        defs = Dict(x => 0.0 for x in dummyder)
+        @set! ssys.defaults = merge(defs, defaults(ssys))
+        op = merge(defs, op)
+    end
+    sys = ssys
+    x0 = merge(defaults(sys), op)
+    u0, p, _ = get_u0_p(sys, x0, p; use_union = false, tofloat = true)
+    p, split_idxs = split_parameters_by_type(p)
+    ps = parameters(sys)
+    if p isa Tuple
+        ps = Base.Fix1(getindex, ps).(split_idxs)
+        ps = (ps...,) #if p is Tuple, ps should be Tuple
+    end
+
     lin_fun = let diff_idxs = diff_idxs,
         alge_idxs = alge_idxs,
         input_idxs = input_idxs,
         sts = states(sys),
-        fun = ODEFunction{true, SciMLBase.FullSpecialize}(sys),
+        fun = ODEFunction{true, SciMLBase.FullSpecialize}(sys, states(sys), ps; p = p),
         h = build_explicit_observed_function(sys, outputs),
         chunk = ForwardDiff.Chunk(input_idxs)
 
@@ -1599,11 +1620,12 @@ function linearize(sys, inputs, outputs; op = Dict(), t = 0.0,
     allow_input_derivatives = false,
     zero_dummy_der = false,
     kwargs...)
-    lin_fun, ssys = linearization_function(sys, inputs, outputs; kwargs...)
-    if zero_dummy_der
-        dummyder = setdiff(states(ssys), states(sys))
-        op = merge(op, Dict(x => 0.0 for x in dummyder))
-    end
+    lin_fun, ssys = linearization_function(sys,
+        inputs,
+        outputs;
+        zero_dummy_der,
+        op,
+        kwargs...)
     linearize(ssys, lin_fun; op, t, allow_input_derivatives), ssys
 end
 
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 8528dd3c12..f11d1283c6 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -354,6 +354,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, dvs = s
     checkbounds = false,
     sparsity = false,
     analytic = nothing,
+    split_idxs = nothing,
     kwargs...) where {iip, specialize}
     f_gen = generate_function(sys, dvs, ps; expression = Val{eval_expression},
         expression_module = eval_module, checkbounds = checkbounds,
@@ -508,6 +509,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, dvs = s
         nothing
     end
 
+    @set! sys.split_idxs = split_idxs
     ODEFunction{iip, specialize}(f;
         sys = sys,
         jac = _jac === nothing ? nothing : _jac,
@@ -765,7 +767,7 @@ Take dictionaries with initial conditions and parameters and convert them to num
 """
 function get_u0_p(sys,
     u0map,
-    parammap;
+    parammap = nothing;
     use_union = true,
     tofloat = true,
     symbolic_u0 = false)
@@ -773,7 +775,9 @@ function get_u0_p(sys,
     ps = parameters(sys)
 
     defs = defaults(sys)
-    defs = mergedefaults(defs, parammap, ps)
+    if parammap !== nothing
+        defs = mergedefaults(defs, parammap, ps)
+    end
     defs = mergedefaults(defs, u0map, dvs)
 
     if symbolic_u0
@@ -835,7 +839,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
     f = constructor(sys, dvs, ps, u0; ddvs = ddvs, tgrad = tgrad, jac = jac,
         checkbounds = checkbounds, p = p,
         linenumbers = linenumbers, parallel = parallel, simplify = simplify,
-        sparse = sparse, eval_expression = eval_expression,
+        sparse = sparse, eval_expression = eval_expression, split_idxs,
         kwargs...)
     implicit_dae ? (f, du0, u0, p) : (f, u0, p)
 end
diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index dcbecfcfb0..3b828702c2 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -139,6 +139,10 @@ struct ODESystem <: AbstractODESystem
     used for ODAEProblem.
     """
     unknown_states::Union{Nothing, Vector{Any}}
+    """
+    split_idxs: a vector of vectors of indices for the split parameters.
+    """
+    split_idxs::Union{Nothing, Vector{Vector{Int}}}
 
     function ODESystem(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad,
         jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults,
@@ -146,8 +150,8 @@ struct ODESystem <: AbstractODESystem
         devents, metadata = nothing, gui_metadata = nothing,
         tearing_state = nothing,
         substitutions = nothing, complete = false,
-        discrete_subsystems = nothing, unknown_states = nothing;
-        checks::Union{Bool, Int} = true)
+        discrete_subsystems = nothing, unknown_states = nothing,
+        split_idxs = nothing; checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckComponents) > 0
             check_variables(dvs, iv)
             check_parameters(ps, iv)
@@ -161,7 +165,7 @@ struct ODESystem <: AbstractODESystem
             ctrl_jac, Wfact, Wfact_t, name, systems, defaults, torn_matching,
             connector_type, preface, cevents, devents, metadata, gui_metadata,
             tearing_state, substitutions, complete, discrete_subsystems,
-            unknown_states)
+            unknown_states, split_idxs)
     end
 end
 
diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index 8e87d145c2..b37aad5d19 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -293,7 +293,7 @@ dprob = DiscreteProblem(js, u₀map, tspan, parammap)
 function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple, Nothing},
     parammap = DiffEqBase.NullParameters();
     checkbounds = false,
-    use_union = false,
+    use_union = true,
     kwargs...)
     dvs = states(sys)
     ps = parameters(sys)
diff --git a/src/utils.jl b/src/utils.jl
index abea65ab21..ab17bd3d8d 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -657,46 +657,40 @@ function promote_to_concrete(vs; tofloat = true, use_union = true)
     end
     T = eltype(vs)
     if Base.isconcretetype(T) && (!tofloat || T === float(T)) # nothing to do
-        vs
+        return vs
     else
         sym_vs = filter(x -> SymbolicUtils.issym(x) || SymbolicUtils.istree(x), vs)
         isempty(sym_vs) || throw_missingvars_in_sys(sym_vs)
-        C = typeof(first(vs))
-        I = Int8
-        has_int = false
-        has_array = false
-        has_bool = false
-        array_T = nothing
+
+        C = nothing
         for v in vs
-            if v isa AbstractArray
-                has_array = true
-                array_T = typeof(v)
+            E = typeof(v)
+            if E <: Number
+                if tofloat
+                    E = float(E)
+                end
             end
-            E = eltype(v)
-            C = promote_type(C, E)
-            if E <: Integer
-                has_int = true
-                I = promote_type(I, E)
+            if C === nothing
+                C = E
             end
-            if E <: Bool
-                has_bool = true
+            if use_union
+                C = Union{C, E}
+            else
+                @assert C==E "`promote_to_concrete` can't make type $E uniform with $C"
+                C = E
             end
         end
-        if tofloat && !has_array
-            C = float(C)
-        elseif has_array || (use_union && has_int && C !== I)
-            if has_array
-                C = Union{C, array_T}
-            end
-            if has_int
-                C = Union{C, I}
-            end
-            if has_bool
-                C = Union{C, Bool}
+
+        y = similar(vs, C)
+        for i in eachindex(vs)
+            if (vs[i] isa Number) & tofloat
+                y[i] = float(vs[i]) #needed because copyto! can't convert Int to Float automatically
+            else
+                y[i] = vs[i]
             end
-            return copyto!(similar(vs, C), vs)
         end
-        convert.(C, vs)
+
+        return y
     end
 end
 
diff --git a/src/variables.jl b/src/variables.jl
index 4cffd3fdeb..854680998e 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -145,15 +145,23 @@ function SciMLBase.process_p_u0_symbolic(prob::Union{SciMLBase.AbstractDEProblem
                                 " Please use `remake` with the `u0` keyword argument as a vector of values, paying attention to state order."))
     end
 
-    # assemble defaults
-    defs = defaults(prob.f.sys)
-    defs = mergedefaults(defs, prob.p, parameters(prob.f.sys))
-    defs = mergedefaults(defs, p, parameters(prob.f.sys))
-    defs = mergedefaults(defs, prob.u0, states(prob.f.sys))
-    defs = mergedefaults(defs, u0, states(prob.f.sys))
-
-    u0 = varmap_to_vars(u0, states(prob.f.sys); defaults = defs, tofloat = true)
-    p = varmap_to_vars(p, parameters(prob.f.sys); defaults = defs)
+    sys = prob.f.sys
+    defs = defaults(sys)
+    ps = parameters(sys)
+    if has_split_idxs(sys) && (split_idxs = get_split_idxs(sys)) !== nothing
+        for (i, idxs) in enumerate(split_idxs)
+            defs = mergedefaults(defs, prob.p[i], ps[idxs])
+        end
+    else
+        # assemble defaults
+        defs = defaults(sys)
+        defs = mergedefaults(defs, prob.p, ps)
+    end
+    defs = mergedefaults(defs, p, ps)
+    sts = states(sys)
+    defs = mergedefaults(defs, prob.u0, sts)
+    defs = mergedefaults(defs, u0, sts)
+    u0, p, defs = get_u0_p(sys, defs)
 
     return p, u0
 end
diff --git a/test/input_output_handling.jl b/test/input_output_handling.jl
index 12f68f9f04..249a94e9d0 100644
--- a/test/input_output_handling.jl
+++ b/test/input_output_handling.jl
@@ -124,8 +124,8 @@ using ModelingToolkitStandardLibrary.Mechanical.Rotational
 t = ModelingToolkitStandardLibrary.Mechanical.Rotational.t
 @named inertia1 = Inertia(; J = 1)
 @named inertia2 = Inertia(; J = 1)
-@named spring = Spring(; c = 10)
-@named damper = Damper(; d = 3)
+@named spring = Rotational.Spring(; c = 10)
+@named damper = Rotational.Damper(; d = 3)
 @named torque = Torque(; use_support = false)
 @variables y(t) = 0
 eqs = [connect(torque.flange, inertia1.flange_a)
diff --git a/test/latexify/10.tex b/test/latexify/10.tex
index 9a8335bd34..09ec74ab72 100644
--- a/test/latexify/10.tex
+++ b/test/latexify/10.tex
@@ -1,5 +1,5 @@
 \begin{align}
-\frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& \frac{\sigma \left(  - x\left( t \right) + y\left( t \right) \right) \frac{\mathrm{d}}{\mathrm{d}t} \left(  - y\left( t \right) + x\left( t \right) \right)}{\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t}} \\
-0 =&  - y\left( t \right) + \frac{1}{10} \sigma \left( \rho - z\left( t \right) \right) x\left( t \right) \\
-\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t} =& \left( y\left( t \right) \right)^{\frac{2}{3}} x\left( t \right) - \beta z\left( t \right)
+\frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& \frac{\left(  - x\left( t \right) + y\left( t \right) \right) \frac{\mathrm{d}}{\mathrm{d}t} \left( x\left( t \right) - y\left( t \right) \right) \sigma}{\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t}} \\
+0 =&  - y\left( t \right) + \frac{1}{10} x\left( t \right) \left(  - z\left( t \right) + \rho \right) \sigma \\
+\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t} =& \left( y\left( t \right) \right)^{\frac{2}{3}} x\left( t \right) - z\left( t \right) \beta
 \end{align}
diff --git a/test/latexify/20.tex b/test/latexify/20.tex
index e32ea81eb7..9de213aafd 100644
--- a/test/latexify/20.tex
+++ b/test/latexify/20.tex
@@ -1,5 +1,5 @@
 \begin{align}
-\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& \left(  - u(t)_1 + u(t)_2 \right) p_3 \\
-0 =&  - u(t)_2 + \frac{1}{10} \left(  - u(t)_1 + p_1 \right) p_2 p_3 u(t)_1 \\
+\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& p_3 \left(  - u(t)_1 + u(t)_2 \right) \\
+0 =&  - u(t)_2 + \frac{1}{10} \left( p_1 - u(t)_1 \right) p_2 p_3 u(t)_1 \\
 \frac{\mathrm{d}}{\mathrm{d}t} u(t)_3 =& u(t)_2^{\frac{2}{3}} u(t)_1 - p_3 u(t)_3
 \end{align}
diff --git a/test/latexify/30.tex b/test/latexify/30.tex
index bc46f41e9b..b83feeba72 100644
--- a/test/latexify/30.tex
+++ b/test/latexify/30.tex
@@ -1,5 +1,5 @@
 \begin{align}
-\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& \left(  - u(t)_1 + u(t)_2 \right) p_3 \\
-\frac{\mathrm{d}}{\mathrm{d}t} u(t)_2 =&  - u(t)_2 + \frac{1}{10} \left(  - u(t)_1 + p_1 \right) p_2 p_3 u(t)_1 \\
+\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& p_3 \left(  - u(t)_1 + u(t)_2 \right) \\
+\frac{\mathrm{d}}{\mathrm{d}t} u(t)_2 =&  - u(t)_2 + \frac{1}{10} \left( p_1 - u(t)_1 \right) p_2 p_3 u(t)_1 \\
 \frac{\mathrm{d}}{\mathrm{d}t} u(t)_3 =& u(t)_2^{\frac{2}{3}} u(t)_1 - p_3 u(t)_3
 \end{align}
diff --git a/test/odesystem.jl b/test/odesystem.jl
index 56d436d76c..39bcda58de 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -23,7 +23,8 @@ eqs = [D(x) ~ σ * (y - x),
 ModelingToolkit.toexpr.(eqs)[1]
 @named de = ODESystem(eqs; defaults = Dict(x => 1))
 subed = substitute(de, [σ => k])
-@test isequal(sort(parameters(subed), by = string), [k, β, ρ])
+ssort(eqs) = sort(eqs, by = string)
+@test isequal(ssort(parameters(subed)), [k, β, ρ])
 @test isequal(equations(subed),
     [D(x) ~ k * (y - x)
         D(y) ~ (ρ - z) * x - y
@@ -241,7 +242,7 @@ p2 = (k₁ => 0.04,
     k₃ => 1e4)
 tspan = (0.0, 100000.0)
 prob1 = ODEProblem(sys, u0, tspan, p)
-@test prob1.f.sys === sys
+@test prob1.f.sys == sys
 prob12 = ODEProblem(sys, u0, tspan, [0.04, 3e7, 1e4])
 prob13 = ODEProblem(sys, u0, tspan, (0.04, 3e7, 1e4))
 prob14 = ODEProblem(sys, u0, tspan, p2)
@@ -348,14 +349,14 @@ eqs = [0 ~ x + z
     D(accumulation_y) ~ y
     D(accumulation_z) ~ z
     D(x) ~ y]
-@test sort(equations(asys), by = string) == eqs
+@test ssort(equations(asys)) == ssort(eqs)
 @variables ac(t)
 asys = add_accumulations(sys, [ac => (x + y)^2])
 eqs = [0 ~ x + z
     0 ~ x - y
     D(ac) ~ (x + y)^2
     D(x) ~ y]
-@test sort(equations(asys), by = string) == eqs
+@test ssort(equations(asys)) == ssort(eqs)
 
 sys2 = ode_order_lowering(sys)
 M = ModelingToolkit.calculate_massmatrix(sys2)
diff --git a/test/runtests.jl b/test/runtests.jl
index e6387e2cd9..b60d82a6ce 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -56,4 +56,6 @@ end
 @safetestset "FuncAffect Test" include("funcaffect.jl")
 @safetestset "Constants Test" include("constants.jl")
 # Reference tests go Last
-@safetestset "Latexify recipes Test" include("latexify.jl")
+if VERSION >= v"1.9"
+    @safetestset "Latexify recipes Test" include("latexify.jl")
+end
diff --git a/test/split_parameters.jl b/test/split_parameters.jl
index ef8f434dca..d37b9e2c13 100644
--- a/test/split_parameters.jl
+++ b/test/split_parameters.jl
@@ -2,6 +2,29 @@ using ModelingToolkit, Test
 using ModelingToolkitStandardLibrary.Blocks
 using OrdinaryDiffEq
 
+x = [1, 2.0, false, [1, 2, 3], Parameter(1.0)]
+
+y = ModelingToolkit.promote_to_concrete(x)
+@test eltype(y) == Union{Float64, Parameter{Float64}, Vector{Int64}}
+
+y = ModelingToolkit.promote_to_concrete(x; tofloat = false)
+@test eltype(y) == Union{Bool, Float64, Int64, Parameter{Float64}, Vector{Int64}}
+
+x = [1, 2.0, false, [1, 2, 3]]
+y = ModelingToolkit.promote_to_concrete(x)
+@test eltype(y) == Union{Float64, Vector{Int64}}
+
+x = Any[1, 2.0, false]
+y = ModelingToolkit.promote_to_concrete(x; tofloat = false)
+@test eltype(y) == Union{Bool, Float64, Int64}
+
+y = ModelingToolkit.promote_to_concrete(x; use_union = false)
+@test eltype(y) == Float64
+
+x = Float16[1.0, 2.0, 3.0]
+y = ModelingToolkit.promote_to_concrete(x)
+@test eltype(y) == Float16
+
 # ------------------------ Mixed Single Values and Vector
 
 dt = 4e-4
@@ -46,12 +69,25 @@ eqs = [y ~ src.output.u
 @named sys = ODESystem(eqs, t, vars, []; systems = [int, src])
 s = complete(sys)
 sys = structural_simplify(sys)
-prob = ODEProblem(sys, [], (0.0, t_end), [s.src.data => x])
+prob = ODEProblem(sys, [], (0.0, t_end), [s.src.data => x]; tofloat = false)
 @test prob.p isa Tuple{Vector{Float64}, Vector{Int}, Vector{Vector{Float64}}}
 sol = solve(prob, ImplicitEuler());
 @test sol.retcode == ReturnCode.Success
 @test sol[y][end] == x[end]
 
+#TODO: remake becomes more complicated now, how to improve?
+defs = ModelingToolkit.defaults(sys)
+defs[s.src.data] = 2x
+p′ = ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat = false)
+p′, = ModelingToolkit.split_parameters_by_type(p′) #NOTE: we need to ensure this is called now before calling remake()
+prob′ = remake(prob; p = p′)
+sol = solve(prob′, ImplicitEuler());
+@test sol.retcode == ReturnCode.Success
+@test sol[y][end] == 2x[end]
+
+prob′′ = remake(prob; p = [s.src.data => x])
+@test_broken prob′′.p isa Tuple
+
 # ------------------------ Mixed Type Converted to float (default behavior)
 
 vars = @variables y(t)=1 dy(t)=0 ddy(t)=0
@@ -77,3 +113,74 @@ prob = ODEProblem(sys, [], tspan, []; tofloat = false)
 @test prob.p isa Tuple{Vector{Float64}, Vector{Int64}}
 sol = solve(prob, ImplicitEuler());
 @test sol.retcode == ReturnCode.Success
+
+# ------------------------- Bug
+using ModelingToolkit, LinearAlgebra
+using ModelingToolkitStandardLibrary.Mechanical.Rotational
+using ModelingToolkitStandardLibrary.Blocks
+using ModelingToolkitStandardLibrary.Blocks: t
+using ModelingToolkit: connect
+
+"A wrapper function to make symbolic indexing easier"
+function wr(sys)
+    ODESystem(Equation[], ModelingToolkit.get_iv(sys), systems = [sys], name = :a_wrapper)
+end
+indexof(sym, syms) = findfirst(isequal(sym), syms)
+
+# Parameters
+m1 = 1.0
+m2 = 1.0
+k = 10.0 # Spring stiffness
+c = 3.0  # Damping coefficient
+
+@named inertia1 = Inertia(; J = m1)
+@named inertia2 = Inertia(; J = m2)
+@named spring = Spring(; c = k)
+@named damper = Damper(; d = c)
+@named torque = Torque(use_support = false)
+
+function SystemModel(u = nothing; name = :model)
+    eqs = [connect(torque.flange, inertia1.flange_a)
+        connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
+        connect(inertia2.flange_a, spring.flange_b, damper.flange_b)]
+    if u !== nothing
+        push!(eqs, connect(torque.tau, u.output))
+        return @named model = ODESystem(eqs,
+            t;
+            systems = [torque, inertia1, inertia2, spring, damper, u])
+    end
+    ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper], name)
+end
+
+model = SystemModel() # Model with load disturbance
+@named d = Step(start_time = 1.0, duration = 10.0, offset = 0.0, height = 1.0) # Disturbance
+model_outputs = [model.inertia1.w, model.inertia2.w, model.inertia1.phi, model.inertia2.phi] # This is the state realization we want to control
+inputs = [model.torque.tau.u]
+matrices, ssys = ModelingToolkit.linearize(wr(model), inputs, model_outputs)
+
+# Design state-feedback gain using LQR
+# Define cost matrices
+x_costs = [model.inertia1.w => 1.0
+    model.inertia2.w => 1.0
+    model.inertia1.phi => 1.0
+    model.inertia2.phi => 1.0]
+L = randn(1, 4) # Post-multiply by `C` to get the correct input to the controller
+
+# This old definition of MatrixGain will work because the parameter space does not include K (an Array term)
+# @component function MatrixGainAlt(K::AbstractArray; name)
+#     nout, nin = size(K, 1), size(K, 2)
+#     @named input = RealInput(; nin = nin)
+#     @named output = RealOutput(; nout = nout)
+#     eqs = [output.u[i] ~ sum(K[i, j] * input.u[j] for j in 1:nin) for i in 1:nout]
+#     compose(ODESystem(eqs, t, [], []; name = name), [input, output])
+# end
+
+@named state_feedback = MatrixGain(K = -L) # Build negative feedback into the feedback matrix
+@named add = Add(; k1 = 1.0, k2 = 1.0) # To add the control signal and the disturbance
+
+connections = [[state_feedback.input.u[i] ~ model_outputs[i] for i in 1:4]
+    connect(d.output, :d, add.input1)
+    connect(add.input2, state_feedback.output)
+    connect(add.output, :u, model.torque.tau)]
+@named closed_loop = ODESystem(connections, t, systems = [model, state_feedback, add, d])
+S = get_sensitivity(closed_loop, :u)
diff --git a/test/stream_connectors.jl b/test/stream_connectors.jl
index 9502b55be9..e4365876cf 100644
--- a/test/stream_connectors.jl
+++ b/test/stream_connectors.jl
@@ -136,29 +136,30 @@ eqns = [domain_connect(fluid, n1m1.port_a)
 
 @test_nowarn structural_simplify(n1m1Test)
 @unpack source, port_a = n1m1
-@test sort(equations(expand_connections(n1m1)), by = string) == [0 ~ port_a.m_flow
+ssort(eqs) = sort(eqs, by = string)
+@test ssort(equations(expand_connections(n1m1))) == ssort([0 ~ port_a.m_flow
     0 ~ source.port1.m_flow - port_a.m_flow
     source.port1.P ~ port_a.P
     source.port1.P ~ source.P
     source.port1.h_outflow ~ port_a.h_outflow
-    source.port1.h_outflow ~ source.h]
+    source.port1.h_outflow ~ source.h])
 @unpack port_a, port_b = pipe
-@test sort(equations(expand_connections(pipe)), by = string) ==
-      [0 ~ -port_a.m_flow - port_b.m_flow
+@test ssort(equations(expand_connections(pipe))) ==
+      ssort([0 ~ -port_a.m_flow - port_b.m_flow
     0 ~ port_a.m_flow
     0 ~ port_b.m_flow
     port_a.P ~ port_b.P
     port_a.h_outflow ~ instream(port_b.h_outflow)
-    port_b.h_outflow ~ instream(port_a.h_outflow)]
-@test sort(equations(expand_connections(sys)), by = string) ==
-      [0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow
+    port_b.h_outflow ~ instream(port_a.h_outflow)])
+@test ssort(equations(expand_connections(sys))) ==
+      ssort([0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow
     0 ~ pipe.port_b.m_flow + sink.port.m_flow
     n1m1.port_a.P ~ pipe.port_a.P
-    pipe.port_b.P ~ sink.port.P]
-@test sort(equations(expand_connections(n1m1Test)), by = string) ==
-      [0 ~ -pipe.port_a.m_flow - pipe.port_b.m_flow
-    0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow
+    pipe.port_b.P ~ sink.port.P])
+@test ssort(equations(expand_connections(n1m1Test))) ==
+      ssort([0 ~ -pipe.port_a.m_flow - pipe.port_b.m_flow
     0 ~ n1m1.source.port1.m_flow - n1m1.port_a.m_flow
+    0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow
     0 ~ pipe.port_b.m_flow + sink.port.m_flow
     fluid.m_flow ~ 0
     n1m1.port_a.P ~ pipe.port_a.P
@@ -172,7 +173,7 @@ eqns = [domain_connect(fluid, n1m1.port_a)
     pipe.port_b.h_outflow ~ n1m1.port_a.h_outflow
     sink.port.P ~ sink.P
     sink.port.h_outflow ~ sink.h_in
-    sink.port.m_flow ~ -sink.m_flow_in]
+    sink.port.m_flow ~ -sink.m_flow_in])
 
 # N1M2 model and test code.
 function N1M2(; name,
@@ -255,12 +256,12 @@ eqns = [connect(source.port, n2m2.port_a)
 @named sp2 = TwoPhaseFluidPort()
 @named sys = ODESystem([connect(sp1, sp2)], t)
 sys_exp = expand_connections(compose(sys, [sp1, sp2]))
-@test sort(equations(sys_exp), by = string) == [0 ~ -sp1.m_flow - sp2.m_flow
+@test ssort(equations(sys_exp)) == ssort([0 ~ -sp1.m_flow - sp2.m_flow
     0 ~ sp1.m_flow
     0 ~ sp2.m_flow
     sp1.P ~ sp2.P
     sp1.h_outflow ~ ModelingToolkit.instream(sp2.h_outflow)
-    sp2.h_outflow ~ ModelingToolkit.instream(sp1.h_outflow)]
+    sp2.h_outflow ~ ModelingToolkit.instream(sp1.h_outflow)])
 
 # array var
 @connector function VecPin(; name)
@@ -274,15 +275,15 @@ end
 
 @named simple = ODESystem([connect(vp1, vp2, vp3)], t)
 sys = expand_connections(compose(simple, [vp1, vp2, vp3]))
-@test sort(equations(sys), by = string) == sort([0 .~ collect(vp1.i)
-        0 .~ collect(vp2.i)
-        0 .~ collect(vp3.i)
-        vp1.v[1] ~ vp2.v[1]
-        vp1.v[2] ~ vp2.v[2]
-        vp1.v[1] ~ vp3.v[1]
-        vp1.v[2] ~ vp3.v[2]
-        0 ~ -vp1.i[1] - vp2.i[1] - vp3.i[1]
-        0 ~ -vp1.i[2] - vp2.i[2] - vp3.i[2]], by = string)
+@test ssort(equations(sys)) == ssort([0 .~ collect(vp1.i)
+    0 .~ collect(vp2.i)
+    0 .~ collect(vp3.i)
+    vp1.v[1] ~ vp2.v[1]
+    vp1.v[2] ~ vp2.v[2]
+    vp1.v[1] ~ vp3.v[1]
+    vp1.v[2] ~ vp3.v[2]
+    0 ~ -vp1.i[1] - vp2.i[1] - vp3.i[1]
+    0 ~ -vp1.i[2] - vp2.i[2] - vp3.i[2]])
 
 @connector function VectorHeatPort(; name, N = 100, T0 = 0.0, Q0 = 0.0)
     @variables (T(t))[1:N]=T0 (Q(t))[1:N]=Q0 [connect = Flow]