fix: fix linearization and downstream CI

lib/ModelingToolkitBase/src/systems/analysis_points.jl

@@ -880,20 +880,20 @@ function generate_control_function(
         kwargs...
     )
     input_ap_name = canonicalize_ap(sys, input_ap_name)
-    u = []
+    u = SymbolicT[]
     for input_ap in input_ap_name
         sys, (du, _) = open_loop(sys, input_ap)
-        push!(u, du)
+        append!(u, du)
     end
     if dist_ap_name === nothing
         return ModelingToolkitBase.generate_control_function(system_modifier(sys), u; kwargs...)
     end
 
     dist_ap_name = canonicalize_ap(sys, dist_ap_name)
-    d = []
+    d = SymbolicT[]
     for dist_ap in dist_ap_name
         sys, (du, _) = open_loop(sys, dist_ap)
-        push!(d, du)
+        append!(d, du)
     end
 
     return ModelingToolkitBase.generate_control_function(system_modifier(sys), u, d; kwargs...)

lib/ModelingToolkitBase/src/systems/codegen.jl

@@ -1293,7 +1293,7 @@ function build_explicit_observed_function(
                 "For `disturbance_argument=false`, use `disturbance_inputs` as before.",
             :build_explicit_observed_function
         )
-        if known_disturbance_inputs !== nothing
+        if known_disturbance_inputs !== ()
             error("Cannot specify both `disturbance_argument=true` and `known_disturbance_inputs`")
         end
         known_disturbance_inputs = disturbance_inputs

src/linearization.jl

@@ -37,14 +37,14 @@ function linearization_function(
         initializealg = nothing,
         initialization_abstol = 1.0e-5,
         initialization_reltol = 1.0e-3,
-        op = Dict(),
+        op = Dict{SymbolicT, SymbolicT}(),
         p = DiffEqBase.NullParameters(),
         zero_dummy_der = false,
         initialization_solver_alg = nothing,
         autodiff = AutoForwardDiff(),
         eval_expression = false, eval_module = @__MODULE__,
         warn_initialize_determined = true,
-        guesses = Dict(),
+        guesses = Dict{SymbolicT, SymbolicT}(),
         warn_empty_op = true,
         t = 0.0,
         kwargs...
@@ -55,20 +55,34 @@ function linearization_function(
     end
     inputs isa AbstractVector || (inputs = [inputs])
     outputs isa AbstractVector || (outputs = [outputs])
-    inputs = mapreduce(vcat, inputs; init = []) do var
-        symbolic_type(var) == ArraySymbolic() ? collect(var) : [var]
-    end
-    outputs = mapreduce(vcat, outputs; init = []) do var
-        symbolic_type(var) == ArraySymbolic() ? collect(var) : [var]
-    end
     ssys = mtkcompile(sys; inputs, outputs, simplify, kwargs...)
     diff_idxs, alge_idxs = eq_idxs(ssys)
     if zero_dummy_der
         dummyder = setdiff(unknowns(ssys), unknowns(sys))
-        defs = Dict(x => 0.0 for x in dummyder)
-        @set! ssys.defaults = merge(defs, defaults(ssys))
-        op = merge(defs, op)
+        ics = initial_conditions(ssys)
+        for x in dummyder
+            ics[x] = Symbolics.COMMON_ZERO
+        end
+    end
+
+    _inputs = SymbolicT[]
+    _outputs = SymbolicT[]
+    for x in inputs
+        if SU.is_array_shape(SU.shape(x))
+            append!(_inputs, vec(collect(x)::Array{SymbolicT})::Vector{SymbolicT})
+        else
+            push!(_inputs, x)
+        end
     end
+    for x in outputs
+        if SU.is_array_shape(SU.shape(x))
+            append!(_outputs, vec(collect(x)::Array{SymbolicT})::Vector{SymbolicT})
+        else
+            push!(_outputs, x)
+        end
+    end
+    inputs = _inputs
+    outputs = _outputs
     sys = ssys
 
     if initializealg === nothing
@@ -591,10 +605,10 @@ function linearize_symbolic(
         B = f_u
         C = h_x
     else
-        gz = lu(g_z; check = false)
+        gz = lu(Num.(g_z); check = false)
         issuccess(gz) ||
             error("g_z not invertible, this indicates that the DAE is of index > 1.")
-        gzgx = -(gz \ g_x)
+        gzgx = -(gz \ Num.(g_x))
         A = [
             f_x f_z
             gzgx * f_x gzgx * f_z
@@ -605,7 +619,7 @@ function linearize_symbolic(
         ] # The cited paper has zeros in the bottom block, see derivation in https://github.com/SciML/ModelingToolkit.jl/pull/1691 for the correct formula
 
         C = [h_x h_z]
-        Bs = -(gz \ g_u) # This equation differ from the cited paper, the paper is likely wrong since their equaiton leads to a dimension mismatch.
+        Bs = -(gz \ Num.(g_u)) # This equation differ from the cited paper, the paper is likely wrong since their equaiton leads to a dimension mismatch.
         if !iszero(Bs)
             if !allow_input_derivatives
                 der_inds = findall(vec(any(!iszero, Bs, dims = 1)))

test/downstream/inversemodel.jl

@@ -17,7 +17,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
 connect = ModelingToolkit.connect;
 rc = 0.25 # Reference concentration
 
-@component function MixingTank(; name, c0 = 0.8, T0 = 308.5, a1 = 0.2674, a21 = 1.815, a22 = 0.4682, b = 1.5476, k0 = 1.05e14, ϵ = 34.2894)
+@component function MixingTank(; name, c0 = 0.8, T0 = 308.5, a1 = 0.2674, a21 = 1.815, a22 = 0.4682, b = 1.5476, k0 = 1.05e14, ϵ = 34.2894, xc = nothing, xT = nothing)
     pars = @parameters begin
         c0 = c0, [description = "Nominal concentration"]
         T0 = T0, [description = "Nominal temperature"]
@@ -34,11 +34,13 @@ rc = 0.25 # Reference concentration
         c = RealOutput()
         T = RealOutput()
     end
-
+
+    xc = @something(xc, c0)
+    xT = @something(xT, T0)
     vars = @variables begin
         gamma(t), [description = "Reaction speed"]
-        xc(t) = c0, [description = "Concentration"]
-        xT(t) = T0, [description = "Temperature"]
+        xc(t) = xc, [description = "Concentration"]
+        xT(t) = xT, [description = "Temperature"]
         xT_c(t), [description = "Cooling temperature"]
     end
 
@@ -73,7 +75,7 @@ begin
     function RefFilter(; name)
         sys = System(Fss; name)
         "Compute initial state that yields y0 as output"
-        empty!(ModelingToolkit.get_defaults(sys))
+        empty!(ModelingToolkit.get_initial_conditions(sys))
         return sys
     end
 end
@@ -96,7 +98,7 @@ end
     systems = @named begin
         ref = Constant(k = 0.25) # Concentration reference
         ff_gain = Gain(k = 1) # To allow turning ff off
-        controller = PI(gainPI.k = 10, T = 500)
+        controller = PI(k = 10, T = 500)
         tank = MixingTank(xc = c_start, xT = T_start, c0 = c0, T0 = T0)
         inverse_tank = MixingTank(xc = nothing, xT = T_start, c0 = c0, T0 = T0)
         feedback = Feedback()
@@ -132,7 +134,9 @@ cm = complete(model)
 
 op = Dict(
     cm.filter.y.u => 0.8 * (1 - 0.42),
-    D(cm.filter.y.u) => 0
+    D(cm.filter.y.u) => 0,
+    cm.noise_filter.x => nothing,
+    cm.controller.int.x => nothing,
 )
 tspan = (0.0, 1000.0)
 # https://github.com/SciML/ModelingToolkit.jl/issues/2786
@@ -178,18 +182,23 @@ nsys = get_named_comp_sensitivity(model, :y; op)
 @testset "Issue #3319" begin
     op1 = Dict(
         cm.filter.y.u => 0.8 * (1 - 0.42),
-        cm.tank.xc => 0.65
+        cm.tank.xc => 0.65,
+        cm.noise_filter.x => nothing,
     )
 
     op2 = Dict(
         cm.filter.y.u => 0.8 * (1 - 0.42),
-        cm.tank.xc => 0.45
+        cm.tank.xc => 0.45,
+        cm.noise_filter.x => nothing,
     )
 
     output = :y
     # we need to provide `op` so the initialization system knows which
     # values to hold constant
-    lin_fun, ssys = get_sensitivity_function(model, output; op = op1)
+    #
+    # Printing `lin_fun` in a tuple (so it doesn't hit the pretty-printing method)
+    # will likely segfault Julia. Somehow.
+    lin_fun, ssys = get_sensitivity_function(model, output; op = op1);
     matrices1, extras1 = linearize(ssys, lin_fun, op = op1)
     matrices2, extras2 = linearize(ssys, lin_fun, op = op2)
     @test extras1.x != extras2.x

test/downstream/linearization_dd.jl

@@ -7,16 +7,23 @@ using ModelingToolkitStandardLibrary.Blocks
 using ModelingToolkitStandardLibrary.Mechanical.MultiBody2D
 using ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition
 using Test
+import Symbolics
 import NonlinearSolve
+using Setfield: @set
 
 using ControlSystemsMTK
 using ControlSystemsMTK.ControlSystemsBase: sminreal, minreal, poles
 connect = ModelingToolkit.connect
 
+function rm_bindings(sys)
+    @set sys.bindings = empty(bindings(sys))
+end
+
 @independent_variables t
 D = Differential(t)
 
 @named link1 = Link(; m = 0.2, l = 10, I = 1, g = -9.807)
+link1 = rm_bindings(link1)
 @named cart = TranslationalPosition.Mass(; m = 1, s = 0)
 @named fixed = Fixed()
 @named force = Force(use_support = false)
@@ -64,13 +71,13 @@ lsyss,
 
 dummyder = setdiff(unknowns(sysss), unknowns(model))
 # op2 = merge(ModelingToolkit.guesses(model), op, Dict(x => 0.0 for x in dummyder))
-op2 = merge(ModelingToolkit.defaults(syss), op)
+op2 = merge(ModelingToolkit.initial_conditions(syss), op)
 op2[link1.fy1] = -op2[link1.g] * op2[link1.m]
 op2[cart.f] = 0
 
-@test substitute(lsyss.A, op2) ≈ lsys.A
+@test Symbolics.value.(substitute(lsyss.A, op2; fold = Val(true))) ≈ lsys.A
 # We cannot pivot symbolically, so the part where a linear solve is required
 # is not reliable.
-@test substitute(lsyss.B, op2)[1:6, 1] ≈ lsys.B[1:6, 1]
-@test substitute(lsyss.C, op2) ≈ lsys.C
-@test substitute(lsyss.D, op2) ≈ lsys.D
+@test Symbolics.value.(substitute(lsyss.B, op2; fold = Val(true)))[1:6, 1] ≈ lsys.B[1:6, 1]
+@test Symbolics.value.(substitute(lsyss.C, op2; fold = Val(true))) ≈ lsys.C
+@test Symbolics.value.(substitute(lsyss.D, op2; fold = Val(true))) ≈ lsys.D