From fdca4e9c1168be229d14aa68680dbfe2a5326f7c Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Mon, 29 Jan 2024 15:53:10 +0530
Subject: [PATCH 1/9] refactor!: require systems to be `complete`d before
 creating an `XProblem`

---
 ext/MTKBifurcationKitExt.jl                    |  6 ++++++
 src/systems/diffeqs/abstractodesystem.jl       | 15 +++++++++++++++
 src/systems/diffeqs/sdesystem.jl               |  3 +++
 src/systems/jumps/jumpsystem.jl                |  6 ++++++
 src/systems/nonlinear/nonlinearsystem.jl       |  3 +++
 src/systems/optimization/optimizationsystem.jl |  3 +++
 6 files changed, 36 insertions(+)

diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl
index 7174f56c8b..3ef483a3db 100644
--- a/ext/MTKBifurcationKitExt.jl
+++ b/ext/MTKBifurcationKitExt.jl
@@ -88,6 +88,9 @@ function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem,
         record_from_solution = BifurcationKit.record_sol_default,
         jac = true,
         kwargs...)
+    if !ModelingToolkit.iscomplete(nsys)
+        error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `BifurcationProblem`")
+    end
     # Creates F and J functions.
     ofun = NonlinearFunction(nsys; jac = jac)
     F = ofun.f
@@ -133,6 +136,9 @@ end
 
 # When input is a ODESystem.
 function BifurcationKit.BifurcationProblem(osys::ODESystem, args...; kwargs...)
+    if !ModelingToolkit.iscomplete(osys)
+        error("A completed `ODESystem` is required. Call `complete` or `structural_simplify` on the system before creating a `BifurcationProblem`")
+    end
     nsys = NonlinearSystem([0 ~ eq.rhs for eq in equations(osys)],
         unknowns(osys),
         parameters(osys);
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 7022d5c7c1..00b132e108 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -891,6 +891,9 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map =
         callback = nothing,
         check_length = true,
         kwargs...) where {iip, specialize}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`")
+    end
     f, u0, p = process_DEProblem(ODEFunction{iip, specialize}, sys, u0map, parammap;
         t = tspan !== nothing ? tspan[1] : tspan,
         check_length, kwargs...)
@@ -959,6 +962,9 @@ end
 function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan,
         parammap = DiffEqBase.NullParameters();
         check_length = true, kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `DAEProblem`")
+    end
     f, du0, u0, p = process_DEProblem(DAEFunction{iip}, sys, u0map, parammap;
         implicit_dae = true, du0map = du0map, check_length,
         kwargs...)
@@ -984,6 +990,9 @@ function DiffEqBase.DDEProblem{iip}(sys::AbstractODESystem, u0map = [],
         callback = nothing,
         check_length = true,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `DDEProblem`")
+    end
     f, u0, p = process_DEProblem(DDEFunction{iip}, sys, u0map, parammap;
         t = tspan !== nothing ? tspan[1] : tspan,
         symbolic_u0 = true,
@@ -1042,6 +1051,9 @@ function DiffEqBase.SDDEProblem{iip}(sys::AbstractODESystem, u0map = [],
         check_length = true,
         sparsenoise = nothing,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `SDDEProblem`")
+    end
     f, u0, p = process_DEProblem(SDDEFunction{iip}, sys, u0map, parammap;
         t = tspan !== nothing ? tspan[1] : tspan,
         symbolic_u0 = true,
@@ -1216,6 +1228,9 @@ end
 function DiffEqBase.SteadyStateProblem{iip}(sys::AbstractODESystem, u0map,
         parammap = SciMLBase.NullParameters();
         check_length = true, kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `SteadyStateProblem`")
+    end
     f, u0, p = process_DEProblem(ODEFunction{iip}, sys, u0map, parammap;
         steady_state = true,
         check_length, kwargs...)
diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl
index 3527477c5a..c8bfa7a6f5 100644
--- a/src/systems/diffeqs/sdesystem.jl
+++ b/src/systems/diffeqs/sdesystem.jl
@@ -573,6 +573,9 @@ function DiffEqBase.SDEProblem{iip}(sys::SDESystem, u0map = [], tspan = get_tspa
         parammap = DiffEqBase.NullParameters();
         sparsenoise = nothing, check_length = true,
         callback = nothing, kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `SDESystem` is required. Call `complete` or `structural_simplify` on the system before creating an `SDEProblem`")
+    end
     f, u0, p = process_DEProblem(SDEFunction{iip}, sys, u0map, parammap; check_length,
         kwargs...)
     cbs = process_events(sys; callback)
diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index 7f0e43a285..c09a454436 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -297,6 +297,9 @@ function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple,
         checkbounds = false,
         use_union = true,
         kwargs...)
+    if !iscomplete(sys)
+        error("A completed `JumpSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `DiscreteProblem`")
+    end
     dvs = unknowns(sys)
     ps = parameters(sys)
 
@@ -388,6 +391,9 @@ sol = solve(jprob, SSAStepper())
 """
 function JumpProcesses.JumpProblem(js::JumpSystem, prob, aggregator; callback = nothing,
         kwargs...)
+    if !iscomplete(js)
+        error("A completed `JumpSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `JumpProblem`")
+    end
     unknowntoid = Dict(value(unknown) => i for (i, unknown) in enumerate(unknowns(js)))
     eqs = equations(js)
     invttype = prob.tspan[1] === nothing ? Float64 : typeof(1 / prob.tspan[2])
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index f040f907f1..bb1006031c 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -366,6 +366,9 @@ end
 function DiffEqBase.NonlinearProblem{iip}(sys::NonlinearSystem, u0map,
         parammap = DiffEqBase.NullParameters();
         check_length = true, kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `NonlinearProblem`")
+    end
     f, u0, p = process_NonlinearProblem(NonlinearFunction{iip}, sys, u0map, parammap;
         check_length, kwargs...)
     pt = something(get_metadata(sys), StandardNonlinearProblem())
diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl
index e520c7490f..5f34197807 100644
--- a/src/systems/optimization/optimizationsystem.jl
+++ b/src/systems/optimization/optimizationsystem.jl
@@ -227,6 +227,9 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
         linenumbers = true, parallel = SerialForm(),
         use_union = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `OptimizationSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `OptimizationProblem`")
+    end
     if haskey(kwargs, :lcons) || haskey(kwargs, :ucons)
         Base.depwarn("`lcons` and `ucons` are deprecated. Specify constraints directly instead.",
             :OptimizationProblem, force = true)

From 4941f47bd861b69082e4a0b19dfdf6728b79143d Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Mon, 29 Jan 2024 16:08:35 +0530
Subject: [PATCH 2/9] refactor!: require systems to be `complete`d before
 creating an `XProblemExpr`

---
 src/systems/diffeqs/abstractodesystem.jl       | 9 +++++++++
 src/systems/diffeqs/sdesystem.jl               | 3 +++
 src/systems/jumps/jumpsystem.jl                | 3 +++
 src/systems/nonlinear/nonlinearsystem.jl       | 3 +++
 src/systems/optimization/optimizationsystem.jl | 3 +++
 5 files changed, 21 insertions(+)

diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 00b132e108..8b55e0ee86 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -1138,6 +1138,9 @@ struct ODEProblemExpr{iip} end
 function ODEProblemExpr{iip}(sys::AbstractODESystem, u0map, tspan,
         parammap = DiffEqBase.NullParameters(); check_length = true,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `ODEProblemExpr`")
+    end
     f, u0, p = process_DEProblem(ODEFunctionExpr{iip}, sys, u0map, parammap; check_length,
         kwargs...)
     linenumbers = get(kwargs, :linenumbers, true)
@@ -1180,6 +1183,9 @@ struct DAEProblemExpr{iip} end
 function DAEProblemExpr{iip}(sys::AbstractODESystem, du0map, u0map, tspan,
         parammap = DiffEqBase.NullParameters(); check_length = true,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `DAEProblemExpr`")
+    end
     f, du0, u0, p = process_DEProblem(DAEFunctionExpr{iip}, sys, u0map, parammap;
         implicit_dae = true, du0map = du0map, check_length,
         kwargs...)
@@ -1260,6 +1266,9 @@ function SteadyStateProblemExpr{iip}(sys::AbstractODESystem, u0map,
         parammap = SciMLBase.NullParameters();
         check_length = true,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `SteadyStateProblemExpr`")
+    end
     f, u0, p = process_DEProblem(ODEFunctionExpr{iip}, sys, u0map, parammap;
         steady_state = true,
         check_length, kwargs...)
diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl
index c8bfa7a6f5..255468d776 100644
--- a/src/systems/diffeqs/sdesystem.jl
+++ b/src/systems/diffeqs/sdesystem.jl
@@ -635,6 +635,9 @@ function SDEProblemExpr{iip}(sys::SDESystem, u0map, tspan,
         parammap = DiffEqBase.NullParameters();
         sparsenoise = nothing, check_length = true,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `SDESystem` is required. Call `complete` or `structural_simplify` on the system before creating an `SDEProblemExpr`")
+    end
     f, u0, p = process_DEProblem(SDEFunctionExpr{iip}, sys, u0map, parammap; check_length,
         kwargs...)
     linenumbers = get(kwargs, :linenumbers, true)
diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index c09a454436..de4397e917 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -356,6 +356,9 @@ function DiscreteProblemExpr{iip}(sys::JumpSystem, u0map, tspan::Union{Tuple, No
         parammap = DiffEqBase.NullParameters();
         use_union = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `JumpSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `DiscreteProblemExpr`")
+    end
     dvs = unknowns(sys)
     ps = parameters(sys)
     defs = defaults(sys)
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index bb1006031c..4f94e8971c 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -399,6 +399,9 @@ function NonlinearProblemExpr{iip}(sys::NonlinearSystem, u0map,
         parammap = DiffEqBase.NullParameters();
         check_length = true,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `NonlinearProblemExpr`")
+    end
     f, u0, p = process_NonlinearProblem(NonlinearFunctionExpr{iip}, sys, u0map, parammap;
         check_length, kwargs...)
     linenumbers = get(kwargs, :linenumbers, true)
diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl
index 5f34197807..0f26557f31 100644
--- a/src/systems/optimization/optimizationsystem.jl
+++ b/src/systems/optimization/optimizationsystem.jl
@@ -422,6 +422,9 @@ function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0map,
         linenumbers = false, parallel = SerialForm(),
         use_union = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `OptimizationSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `OptimizationProblemExpr`")
+    end
     if haskey(kwargs, :lcons) || haskey(kwargs, :ucons)
         Base.depwarn("`lcons` and `ucons` are deprecated. Specify constraints directly instead.",
             :OptimizationProblem, force = true)

From 4983626b488fb408b2c605c3bbd4142a7fad0e83 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Mon, 29 Jan 2024 16:58:57 +0530
Subject: [PATCH 3/9] docs: call `complete` where required in documentation
 before creating `XProblem`s

---
 docs/src/basics/Events.md                             | 4 ++--
 docs/src/examples/parsing.md                          | 2 +-
 docs/src/examples/sparse_jacobians.md                 | 1 +
 docs/src/tutorials/bifurcation_diagram_computation.md | 2 ++
 docs/src/tutorials/modelingtoolkitize.md              | 2 +-
 docs/src/tutorials/nonlinear.md                       | 1 +
 docs/src/tutorials/optimization.md                    | 4 ++--
 docs/src/tutorials/stochastic_diffeq.md               | 1 +
 src/systems/abstractsystem.jl                         | 2 ++
 src/systems/diffeqs/basic_transformations.jl          | 2 +-
 src/systems/diffeqs/sdesystem.jl                      | 2 +-
 src/systems/jumps/jumpsystem.jl                       | 6 +++---
 12 files changed, 18 insertions(+), 11 deletions(-)

diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md
index e23b1a6381..3cd0e85f70 100644
--- a/docs/src/basics/Events.md
+++ b/docs/src/basics/Events.md
@@ -73,7 +73,7 @@ function UnitMassWithFriction(k; name)
     ODESystem(eqs, t; continuous_events = [v ~ 0], name) # when v = 0 there is a discontinuity
 end
 @named m = UnitMassWithFriction(0.7)
-prob = ODEProblem(m, Pair[], (0, 10pi))
+prob = ODEProblem(complete(m), Pair[], (0, 10pi))
 sol = solve(prob, Tsit5())
 plot(sol)
 ```
@@ -244,7 +244,7 @@ u0 = [N => 0.0]
 tspan = (0.0, 20.0)
 p = [α => 100.0, tinject => 10.0, M => 50]
 @named osys = ODESystem(eqs, t, [N], [α, M, tinject]; discrete_events = injection)
-oprob = ODEProblem(osys, u0, tspan, p)
+oprob = ODEProblem(complete(osys), u0, tspan, p)
 sol = solve(oprob, Tsit5(); tstops = 10.0)
 plot(sol)
 ```
diff --git a/docs/src/examples/parsing.md b/docs/src/examples/parsing.md
index 47d5321a4e..e23612178f 100644
--- a/docs/src/examples/parsing.md
+++ b/docs/src/examples/parsing.md
@@ -27,6 +27,6 @@ using ModelingToolkit, NonlinearSolve
 vars = union(ModelingToolkit.vars.(eqs)...)
 @named ns = NonlinearSystem(eqs, vars, [])
 
-prob = NonlinearProblem(ns, [1.0, 1.0, 1.0])
+prob = NonlinearProblem(complete(ns), [1.0, 1.0, 1.0])
 sol = solve(prob, NewtonRaphson())
 ```
diff --git a/docs/src/examples/sparse_jacobians.md b/docs/src/examples/sparse_jacobians.md
index 27bf628357..b7d96ac569 100644
--- a/docs/src/examples/sparse_jacobians.md
+++ b/docs/src/examples/sparse_jacobians.md
@@ -55,6 +55,7 @@ Now let's use `modelingtoolkitize` to generate the symbolic version:
 
 ```@example sparsejac
 sys = modelingtoolkitize(prob);
+sys = complete(sys);
 nothing # hide
 ```
 
diff --git a/docs/src/tutorials/bifurcation_diagram_computation.md b/docs/src/tutorials/bifurcation_diagram_computation.md
index 572e032256..ccc66e37e0 100644
--- a/docs/src/tutorials/bifurcation_diagram_computation.md
+++ b/docs/src/tutorials/bifurcation_diagram_computation.md
@@ -13,6 +13,7 @@ using ModelingToolkit
 eqs = [0 ~ μ * x - x^3 + α * y,
     0 ~ -y]
 @named nsys = NonlinearSystem(eqs, [x, y], [μ, α])
+nsys = complete(nsys)
 ```
 
 we wish to compute a bifurcation diagram for this system as we vary the parameter `μ`. For this, we need to provide the following information:
@@ -97,6 +98,7 @@ D = Differential(t)
 eqs = [D(x) ~ μ * x - y - x * (x^2 + y^2),
     D(y) ~ x + μ * y - y * (x^2 + y^2)]
 @named osys = ODESystem(eqs, t)
+osys = complete(osys)
 
 bif_par = μ
 plot_var = x
diff --git a/docs/src/tutorials/modelingtoolkitize.md b/docs/src/tutorials/modelingtoolkitize.md
index 35c5acd9bf..e6e9cca046 100644
--- a/docs/src/tutorials/modelingtoolkitize.md
+++ b/docs/src/tutorials/modelingtoolkitize.md
@@ -58,5 +58,5 @@ sys = modelingtoolkitize(prob)
 Using this, we can symbolically build the Jacobian and then rebuild the ODEProblem:
 
 ```@example mtkize
-prob_jac = ODEProblem(sys, [], (0.0, 1e5), jac = true)
+prob_jac = ODEProblem(complete(sys), [], (0.0, 1e5), jac = true)
 ```
diff --git a/docs/src/tutorials/nonlinear.md b/docs/src/tutorials/nonlinear.md
index 9b8ea954ce..07ccdbcf38 100644
--- a/docs/src/tutorials/nonlinear.md
+++ b/docs/src/tutorials/nonlinear.md
@@ -17,6 +17,7 @@ eqs = [0 ~ σ * (y - x),
     0 ~ x * (ρ - z) - y,
     0 ~ x * y - β * z]
 @named ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])
+ns = complete(ns)
 
 guess = [x => 1.0,
     y => 0.0,
diff --git a/docs/src/tutorials/optimization.md b/docs/src/tutorials/optimization.md
index b4ee03141b..f572988a19 100644
--- a/docs/src/tutorials/optimization.md
+++ b/docs/src/tutorials/optimization.md
@@ -50,7 +50,7 @@ u0 = [x => 1.0
 p = [a => 1.0
     b => 100.0]
 
-prob = OptimizationProblem(sys, u0, p, grad = true, hess = true)
+prob = OptimizationProblem(complete(sys), u0, p, grad = true, hess = true)
 solve(prob, GradientDescent())
 ```
 
@@ -71,7 +71,7 @@ cons = [
 @named sys = OptimizationSystem(loss, [x, y], [a, b], constraints = cons)
 u0 = [x => 0.14
     y => 0.14]
-prob = OptimizationProblem(sys, u0, grad = true, hess = true, cons_j = true, cons_h = true)
+prob = OptimizationProblem(complete(sys), u0, grad = true, hess = true, cons_j = true, cons_h = true)
 solve(prob, IPNewton())
 ```
 
diff --git a/docs/src/tutorials/stochastic_diffeq.md b/docs/src/tutorials/stochastic_diffeq.md
index aa81449313..0cb45c73d6 100644
--- a/docs/src/tutorials/stochastic_diffeq.md
+++ b/docs/src/tutorials/stochastic_diffeq.md
@@ -24,6 +24,7 @@ noiseeqs = [0.1 * x,
     0.1 * z]
 
 @named de = SDESystem(eqs, noiseeqs, t, [x, y, z], [σ, ρ, β])
+de = complete(de)
 
 u0map = [
     x => 1.0,
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 3cfc7dcb16..c599154204 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -1330,6 +1330,8 @@ information. E.g.
 ```julia-repl
 julia> sys = debug_system(sys);
 
+julia> sys = complete(sys);
+
 julia> prob = ODEProblem(sys, [], (0, 1.0));
 
 julia> du = zero(prob.u0);
diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index abd3388a0b..61682c806e 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -31,7 +31,7 @@ u0 = [x => 1.0,
       y => 1.0,
       trJ => 1.0]
 
-prob = ODEProblem(sys2,u0,tspan,p)
+prob = ODEProblem(complete(sys2),u0,tspan,p)
 sol = solve(prob,Tsit5())
 ```
 
diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl
index 255468d776..8a845ab0fe 100644
--- a/src/systems/diffeqs/sdesystem.jl
+++ b/src/systems/diffeqs/sdesystem.jl
@@ -322,7 +322,7 @@ parammap = [
     β => 1.0
 ]
 
-probmod = SDEProblem(demod,u0modmap,(0.0,1.0),parammap)
+probmod = SDEProblem(complete(demod),u0modmap,(0.0,1.0),parammap)
 ensemble_probmod = EnsembleProblem(probmod;
           output_func = (sol,i) -> (g(sol[x,end])*sol[demod.weight,end],false),
           )
diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index de4397e917..f81dd8b0d2 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -289,7 +289,7 @@ using DiffEqBase, JumpProcesses
 u₀map = [S => 999, I => 1, R => 0]
 parammap = [β => 0.1 / 1000, γ => 0.01]
 tspan = (0.0, 250.0)
-dprob = DiscreteProblem(js, u₀map, tspan, parammap)
+dprob = DiscreteProblem(complete(js), u₀map, tspan, parammap)
 ```
 """
 function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple, Nothing},
@@ -347,7 +347,7 @@ using DiffEqBase, JumpProcesses
 u₀map = [S => 999, I => 1, R => 0]
 parammap = [β => 0.1 / 1000, γ => 0.01]
 tspan = (0.0, 250.0)
-dprob = DiscreteProblem(js, u₀map, tspan, parammap)
+dprob = DiscreteProblem(complete(js), u₀map, tspan, parammap)
 ```
 """
 struct DiscreteProblemExpr{iip} end
@@ -388,7 +388,7 @@ Generates a JumpProblem from a JumpSystem.
 Continuing the example from the [`DiscreteProblem`](@ref) definition:
 
 ```julia
-jprob = JumpProblem(js, dprob, Direct())
+jprob = JumpProblem(complete(js), dprob, Direct())
 sol = solve(jprob, SSAStepper())
 ```
 """

From caead09d6d5f6c4ddd1aa026cf3ae5e2047d6683 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Mon, 29 Jan 2024 21:21:16 +0530
Subject: [PATCH 4/9] TEMP: change version so CI runs

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 2d327d1d0b..2a4c33f0ce 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <mayingbo5@gmail.com>", "Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "9.0.0"
+version = "8.76.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

From a7f6c6fc882a20027154afd25f9befe0d21fa83c Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Tue, 30 Jan 2024 19:25:44 +0530
Subject: [PATCH 5/9] fixup! refactor!: require systems to be `complete`d
 before creating an `XProblem`

---
 ext/MTKBifurcationKitExt.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl
index 3ef483a3db..20a8aaa6bd 100644
--- a/ext/MTKBifurcationKitExt.jl
+++ b/ext/MTKBifurcationKitExt.jl
@@ -143,7 +143,7 @@ function BifurcationKit.BifurcationProblem(osys::ODESystem, args...; kwargs...)
         unknowns(osys),
         parameters(osys);
         name = nameof(osys))
-    return BifurcationKit.BifurcationProblem(nsys, args...; kwargs...)
+    return BifurcationKit.BifurcationProblem(complete(nsys), args...; kwargs...)
 end
 
 end # module

From 2282053c8209f239686941a11ac646fbf32a27d3 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Tue, 30 Jan 2024 19:26:01 +0530
Subject: [PATCH 6/9] fix: add constructor so `JumpSystem`s can be `complete`d

---
 src/systems/jumps/jumpsystem.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index f81dd8b0d2..f349a1dcf0 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -119,6 +119,7 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem
             connector_type, devents, metadata, gui_metadata, complete)
     end
 end
+JumpSystem(tag, ap, iv, states, ps, var_to_name, args...; kwargs...) = JumpSystem{typeof(ap)}(tag, ap, iv, states, ps, var_to_name, args...; kwargs...)
 
 function JumpSystem(eqs, iv, unknowns, ps;
         observed = Equation[],

From c555a29c4c2a1abf75f868279191fee3a23e882f Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Tue, 30 Jan 2024 19:26:20 +0530
Subject: [PATCH 7/9] fix: call `complete` when `OptimizationSystem`s are
 simplified

---
 src/systems/optimization/optimizationsystem.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl
index 0f26557f31..0ea8348d02 100644
--- a/src/systems/optimization/optimizationsystem.jl
+++ b/src/systems/optimization/optimizationsystem.jl
@@ -621,5 +621,6 @@ function structural_simplify(sys::OptimizationSystem; kwargs...)
     neweqs = fixpoint_sub.(equations(sys), (subs,))
     @set! sys.op = length(neweqs) == 1 ? first(neweqs) : neweqs
     @set! sys.unknowns = newsts
+    sys = complete(sys)
     return sys
 end

From 9b416177d66e5a2035d0c5be8bcbda9e5864f5ca Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Tue, 30 Jan 2024 19:26:49 +0530
Subject: [PATCH 8/9] test: fix tests

---
 test/constants.jl                             |  2 +-
 test/dae_jacobian.jl                          |  2 +-
 test/extensions/bifurcationkit.jl             |  4 +--
 test/funcaffect.jl                            | 10 ++++----
 test/jacobiansparsity.jl                      |  4 +--
 test/jumpsystem.jl                            |  7 ++++++
 test/lowering_solving.jl                      |  4 ++-
 test/modelingtoolkitize.jl                    | 10 ++++----
 test/nonlinearsystem.jl                       |  7 ++++--
 test/odesystem.jl                             | 11 ++++++--
 test/optimizationsystem.jl                    | 25 ++++++++++---------
 test/sdesystem.jl                             | 10 +++++---
 test/serialization.jl                         |  5 ++--
 test/steadystatesystems.jl                    |  1 +
 .../index_reduction.jl                        |  4 +--
 test/structural_transformation/tearing.jl     |  2 +-
 test/symbolic_events.jl                       | 10 +++++---
 test/symbolic_parameters.jl                   |  4 ++-
 18 files changed, 77 insertions(+), 45 deletions(-)

diff --git a/test/constants.jl b/test/constants.jl
index ee9c9fa618..373c37464a 100644
--- a/test/constants.jl
+++ b/test/constants.jl
@@ -10,7 +10,7 @@ UMT = ModelingToolkit.UnitfulUnitCheck
 D = Differential(t)
 eqs = [D(x) ~ a]
 @named sys = ODESystem(eqs)
-prob = ODEProblem(sys, [0], [0.0, 1.0], [])
+prob = ODEProblem(complete(sys), [0], [0.0, 1.0], [])
 sol = solve(prob, Tsit5())
 
 newsys = MT.eliminate_constants(sys)
diff --git a/test/dae_jacobian.jl b/test/dae_jacobian.jl
index 9047a91bc2..4ca6ce865c 100644
--- a/test/dae_jacobian.jl
+++ b/test/dae_jacobian.jl
@@ -49,7 +49,7 @@ du0 = [0.5, -2.0]
 p = [p1 => 1.5,
     p2 => 3.0]
 
-prob = DAEProblem(sys, du0, u0, tspan, p, jac = true, sparse = true)
+prob = DAEProblem(complete(sys), du0, u0, tspan, p, jac = true, sparse = true)
 sol = solve(prob, IDA(linear_solver = :KLU))
 
 @test maximum(sol - sol1) < 1e-12
diff --git a/test/extensions/bifurcationkit.jl b/test/extensions/bifurcationkit.jl
index 53c7369485..a3ed2bb813 100644
--- a/test/extensions/bifurcationkit.jl
+++ b/test/extensions/bifurcationkit.jl
@@ -9,7 +9,7 @@ let
     eqs = [0 ~ μ * x - x^3 + α * y,
         0 ~ -y]
     @named nsys = NonlinearSystem(eqs, [x, y], [μ, α])
-
+    nsys = complete(nsys)
     # Creates BifurcationProblem 
     bif_par = μ
     p_start = [μ => -1.0, α => 1.0]
@@ -62,7 +62,7 @@ let
     eqs = [D(x) ~ -x + a * y + x^2 * y,
         D(y) ~ b - a * y - x^2 * y]
     @named sys = ODESystem(eqs)
-
+    sys = complete(sys)
     # Creates BifurcationProblem
     bprob = BifurcationProblem(sys,
         [x => 1.5, y => 1.0],
diff --git a/test/funcaffect.jl b/test/funcaffect.jl
index bb953988b5..782f98a62d 100644
--- a/test/funcaffect.jl
+++ b/test/funcaffect.jl
@@ -11,7 +11,7 @@ affect1!(integ, u, p, ctx) = integ.u[u.u] += 10
 
 @named sys = ODESystem(eqs, t, [u], [],
     discrete_events = [[4.0] => (affect1!, [u], [], nothing)])
-prob = ODEProblem(sys, [u => 10.0], (0, 10.0))
+prob = ODEProblem(complete(sys), [u => 10.0], (0, 10.0))
 sol = solve(prob, Tsit5())
 i4 = findfirst(==(4.0), sol[:t])
 @test sol.u[i4 + 1][1] > 10.0
@@ -62,7 +62,7 @@ end
 ctx1 = [10.0]
 @named sys = ODESystem(eqs, t, [u], [],
     discrete_events = [[4.0, 8.0] => (affect2!, [u], [], ctx1)])
-prob = ODEProblem(sys, [u => 10.0], (0, 10.0))
+prob = ODEProblem(complete(sys), [u => 10.0], (0, 10.0))
 sol = solve(prob, Tsit5())
 i4 = findfirst(==(4.0), sol[:t])
 @test sol.u[i4 + 1][1] > 10.0
@@ -79,7 +79,7 @@ end
 @parameters a = 10.0
 @named sys = ODESystem(eqs, t, [u], [a],
     discrete_events = [[4.0, 8.0] => (affect3!, [u], [a], nothing)])
-prob = ODEProblem(sys, [u => 10.0], (0, 10.0))
+prob = ODEProblem(complete(sys), [u => 10.0], (0, 10.0))
 
 sol = solve(prob, Tsit5())
 i4 = findfirst(==(4.0), sol[:t])
@@ -97,7 +97,7 @@ end
     discrete_events = [
         [4.0, 8.0] => (affect3!, [u], [a => :b], nothing),
     ])
-prob = ODEProblem(sys, [u => 10.0], (0, 10.0))
+prob = ODEProblem(complete(sys), [u => 10.0], (0, 10.0))
 
 sol = solve(prob, Tsit5())
 i4 = findfirst(==(4.0), sol[:t])
@@ -225,7 +225,7 @@ balls = compose(balls, [ball1, ball2])
 @test ModelingToolkit.has_discrete_events(balls)
 @test length(ModelingToolkit.affects(ModelingToolkit.discrete_events(balls))) == 2
 
-prob = ODEProblem(balls, [ball1.x => 10.0, ball1.v => 0, ball2.x => 10.0, ball2.v => 0],
+prob = ODEProblem(complete(balls), [ball1.x => 10.0, ball1.v => 0, ball2.x => 10.0, ball2.v => 0],
     (0, 3.0))
 sol = solve(prob, Tsit5())
 
diff --git a/test/jacobiansparsity.jl b/test/jacobiansparsity.jl
index f427c953c0..fe71e51f2e 100644
--- a/test/jacobiansparsity.jl
+++ b/test/jacobiansparsity.jl
@@ -40,7 +40,7 @@ end
 u0 = init_brusselator_2d(xyd_brusselator)
 prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
     u0, (0.0, 11.5), p)
-sys = modelingtoolkitize(prob_ode_brusselator_2d)
+sys = complete(modelingtoolkitize(prob_ode_brusselator_2d))
 
 # test sparse jacobian pattern only.
 prob = ODEProblem(sys, u0, (0, 11.5), sparse = true, jac = false)
@@ -74,7 +74,7 @@ f = DiffEqBase.ODEFunction(sys, u0 = nothing, sparse = true, jac = false)
 u0 = similar(init_brusselator_2d(xyd_brusselator), Float32)
 prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
     u0, (0.0, 11.5), p)
-sys = modelingtoolkitize(prob_ode_brusselator_2d)
+sys = complete(modelingtoolkitize(prob_ode_brusselator_2d))
 
 prob = ODEProblem(sys, u0, (0, 11.5), sparse = true, jac = false)
 @test eltype(prob.f.jac_prototype) == Float32
diff --git a/test/jumpsystem.jl b/test/jumpsystem.jl
index fdb6356426..9911c777e8 100644
--- a/test/jumpsystem.jl
+++ b/test/jumpsystem.jl
@@ -59,6 +59,7 @@ rate₃ = γ * I * h
 affect₃ = [I ~ I * h - 1, R ~ R + 1]
 j₃ = ConstantRateJump(rate₃, affect₃)
 @named js2 = JumpSystem([j₁, j₃], t, [S, I, R], [β, γ])
+js2 = complete(js2)
 u₀ = [999, 1, 0];
 p = (0.1 / 1000, 0.01);
 tspan = (0.0, 250.0);
@@ -80,6 +81,7 @@ m = getmean(jprob, Nsims)
 @variables S2(t)
 obs = [S2 ~ 2 * S]
 @named js2b = JumpSystem([j₁, j₃], t, [S, I, R], [β, γ], observed = obs)
+js2b = complete(js2b)
 dprob = DiscreteProblem(js2b, u₀map, tspan, parammap)
 jprob = JumpProblem(js2b, dprob, Direct(), save_positions = (false, false), rng = rng)
 sol = solve(jprob, SSAStepper(), saveat = tspan[2] / 10)
@@ -132,6 +134,7 @@ m2 = getmean(jprob, Nsims)
 maj1 = MassActionJump(2 * β / 2, [S => 1, I => 1], [S => -1, I => 1])
 maj2 = MassActionJump(γ, [I => 1], [I => -1, R => 1])
 @named js3 = JumpSystem([maj1, maj2], t, [S, I, R], [β, γ])
+js3 = complete(js3)
 dprob = DiscreteProblem(js3, u₀map, tspan, parammap)
 jprob = JumpProblem(js3, dprob, Direct(), rng = rng)
 m3 = getmean(jprob, Nsims)
@@ -139,6 +142,7 @@ m3 = getmean(jprob, Nsims)
 
 # maj jump test with various dep graphs
 @named js3b = JumpSystem([maj1, maj2], t, [S, I, R], [β, γ])
+js3b = complete(js3b)
 jprobb = JumpProblem(js3b, dprob, NRM(), rng = rng)
 m4 = getmean(jprobb, Nsims)
 @test abs(m - m4) / m < 0.01
@@ -150,6 +154,7 @@ m4 = getmean(jprobc, Nsims)
 maj1 = MassActionJump(2.0, [0 => 1], [S => 1])
 maj2 = MassActionJump(γ, [S => 1], [S => -1])
 @named js4 = JumpSystem([maj1, maj2], t, [S], [β, γ])
+js4 = complete(js4)
 dprob = DiscreteProblem(js4, [S => 999], (0, 1000.0), [β => 100.0, γ => 0.01])
 jprob = JumpProblem(js4, dprob, Direct(), rng = rng)
 m4 = getmean(jprob, Nsims)
@@ -159,6 +164,7 @@ m4 = getmean(jprob, Nsims)
 maj1 = MassActionJump(2.0, [0 => 1], [S => 1])
 maj2 = MassActionJump(γ, [S => 2], [S => -1])
 @named js4 = JumpSystem([maj1, maj2], t, [S], [β, γ])
+js4 = complete(js4)
 dprob = DiscreteProblem(js4, [S => 999], (0, 1000.0), [β => 100.0, γ => 0.01])
 jprob = JumpProblem(js4, dprob, Direct(), rng = rng)
 sol = solve(jprob, SSAStepper());
@@ -177,6 +183,7 @@ end
 maj1 = MassActionJump(k1 * k3, [0 => 1], [A => -1, B => 1])
 maj2 = MassActionJump(k2, [B => 1], [A => 1, B => -1])
 @named js5 = JumpSystem([maj1, maj2], t, [A, B], [k1, k2, k3])
+js5 = complete(js5)
 p = [k1 => 2.0, k2 => 0.0, k3 => 0.5]
 u₀ = [A => 100, B => 0]
 tspan = (0.0, 2000.0)
diff --git a/test/lowering_solving.jl b/test/lowering_solving.jl
index 712a2a3595..58dca5e2b1 100644
--- a/test/lowering_solving.jl
+++ b/test/lowering_solving.jl
@@ -30,6 +30,8 @@ p = [σ => 28.0,
     β => 8 / 3]
 
 tspan = (0.0, 100.0)
+
+sys = complete(sys)
 prob = ODEProblem(sys, u0, tspan, p, jac = true)
 probexpr = ODEProblemExpr(sys, u0, tspan, p, jac = true)
 sol = solve(prob, Tsit5())
@@ -52,7 +54,7 @@ lorenz2 = ODESystem(eqs, name = :lorenz2)
 @parameters γ
 connections = [0 ~ lorenz1.x + lorenz2.y + α * γ]
 @named connected = ODESystem(connections, t, [α], [γ], systems = [lorenz1, lorenz2])
-
+connected = complete(connected)
 u0 = [lorenz1.x => 1.0,
     lorenz1.y => 0.0,
     lorenz1.z => 0.0,
diff --git a/test/modelingtoolkitize.jl b/test/modelingtoolkitize.jl
index fea90491a3..be804d5c71 100644
--- a/test/modelingtoolkitize.jl
+++ b/test/modelingtoolkitize.jl
@@ -52,7 +52,7 @@ x0 = zeros(2)
 p = [1.0, 100.0]
 
 prob = OptimizationProblem(rosenbrock, x0, p)
-sys = modelingtoolkitize(prob) # symbolicitize me captain!
+sys = complete(modelingtoolkitize(prob)) # symbolicitize me captain!
 
 prob = OptimizationProblem(sys, x0, p, grad = true, hess = true)
 sol = solve(prob, NelderMead())
@@ -141,7 +141,7 @@ problem = ODEProblem(SIRD_ac!, ℬ, 𝒯, 𝒫)
 @time solution = solve(problem, Tsit5(), saveat = 1:final_time);
 
 problem = ODEProblem(SIRD_ac!, ℬ, 𝒯, 𝒫)
-sys = modelingtoolkitize(problem)
+sys = complete(modelingtoolkitize(problem))
 fast_problem = ODEProblem(sys, ℬ, 𝒯, 𝒫)
 @time solution = solve(fast_problem, Tsit5(), saveat = 1:final_time)
 
@@ -179,7 +179,7 @@ u0 = [1.0, 0, 0, 0, 0]
 p = [9.8, 1]
 tspan = (0, 10.0)
 pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)
-pendulum_sys_org = modelingtoolkitize(pendulum_prob)
+pendulum_sys_org = complete(modelingtoolkitize(pendulum_prob))
 sts = unknowns(pendulum_sys_org)
 pendulum_sys = dae_index_lowering(pendulum_sys_org)
 prob = ODEProblem(pendulum_sys, Pair[], tspan)
@@ -250,7 +250,7 @@ u0 = @LArray [9998.0, 1.0, 1.0, 1.0] (:S, :I, :R, :C)
 
 # Initiate ODE problem
 problem = ODEProblem(SIR!, u0, tspan, p)
-sys = modelingtoolkitize(problem)
+sys = complete(modelingtoolkitize(problem))
 
 @parameters t
 @test all(isequal.(parameters(sys), getproperty.(@variables(β, η, ω, φ, σ, μ), :val)))
@@ -304,6 +304,6 @@ noiseeqs = [0.1 * x,
     0.1 * z]
 
 @named sys = SDESystem(eqs, noiseeqs, t, [x, y, z], [sig, rho, beta]; tspan = (0, 1000.0))
-prob = SDEProblem(sys)
+prob = SDEProblem(complete(sys))
 sys = modelingtoolkitize(prob)
 @test ModelingToolkit.has_tspan(sys)
diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl
index aa1cd97973..3c329c887e 100644
--- a/test/nonlinearsystem.jl
+++ b/test/nonlinearsystem.jl
@@ -66,6 +66,7 @@ eqs = [0 ~ σ * a * h,
     0 ~ x * (ρ - z) - y,
     0 ~ x * y - β * z]
 @named ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])
+ns = complete(ns)
 nlsys_func = generate_function(ns, [x, y, z], [σ, ρ, β])
 nf = NonlinearFunction(ns)
 jac = calculate_jacobian(ns)
@@ -98,7 +99,7 @@ eqs1 = [
 
 lorenz = name -> NonlinearSystem(eqs1, [x, y, z, u, F], [σ, ρ, β], name = name)
 lorenz1 = lorenz(:lorenz1)
-@test_throws ArgumentError NonlinearProblem(lorenz1, zeros(5))
+@test_throws ArgumentError NonlinearProblem(complete(lorenz1), zeros(5))
 lorenz2 = lorenz(:lorenz2)
 @named connected = NonlinearSystem([s ~ a + lorenz1.x
         lorenz2.y ~ s * h
@@ -128,7 +129,7 @@ eqs = [0 ~ σ * (y - x),
     0 ~ x * (ρ - z) - y,
     0 ~ x * y - β * z * h]
 @named ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])
-np = NonlinearProblem(ns, [0, 0, 0], [1, 2, 3], jac = true, sparse = true)
+np = NonlinearProblem(complete(ns), [0, 0, 0], [1, 2, 3], jac = true, sparse = true)
 @test calculate_jacobian(ns, sparse = true) isa SparseMatrixCSC
 
 # issue #819
@@ -199,6 +200,7 @@ eq = [v1 ~ sin(2pi * t * h)
         u[4] ~ h]
 
     sys = NonlinearSystem(eqs, collect(u[1:4]), Num[], defaults = Dict([]), name = :test)
+    sys = complete(sys)
     prob = NonlinearProblem(sys, ones(length(unknowns(sys))))
 
     sol = NonlinearSolve.solve(prob, NewtonRaphson())
@@ -223,6 +225,7 @@ testdict = Dict([:test => 1])
         0 ~ x * (b - z) - y,
         0 ~ x * y - c * z]
     @named sys = NonlinearSystem(eqs, [x, y, z], [a, b, c], defaults = Dict(x => 2.0))
+    sys = complete(sys)
     prob = NonlinearProblem(sys, ones(length(unknowns(sys))))
 
     prob_ = remake(prob, u0 = [1.0, 2.0, 3.0], p = [1.1, 1.2, 1.3])
diff --git a/test/odesystem.jl b/test/odesystem.jl
index 1320fa7634..5ccbaec536 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -232,6 +232,7 @@ eqs = [D(y₁) ~ -k₁ * y₁ + k₃ * y₂ * y₃,
     0 ~ y₁ + y₂ + y₃ - 1,
     D(y₂) ~ k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃ * κ]
 @named sys = ODESystem(eqs, defaults = [k₁ => 100, k₂ => 3e7, y₁ => 1.0])
+sys = complete(sys)
 u0 = Pair[]
 push!(u0, y₂ => 0.0)
 push!(u0, y₃ => 0.0)
@@ -279,7 +280,7 @@ sol_dpmap = solve(prob_dpmap, Rodas5())
         sys1 = makesys(:sys1)
         sys2 = makesys(:sys2)
         @parameters t b=1.0
-        ODESystem(Equation[], t, [], [b]; systems = [sys1, sys2], name = :foo)
+        complete(ODESystem(Equation[], t, [], [b]; systems = [sys1, sys2], name = :foo))
     end
 
     sys = makecombinedsys()
@@ -435,6 +436,7 @@ default_u0 = [D(x) => 0.0, x => 10.0]
 default_p = [M => 1.0, b => 1.0, k => 1.0]
 @named sys = ODESystem(eqs, t, [x], ps; defaults = [default_u0; default_p], tspan)
 sys = ode_order_lowering(sys)
+sys = complete(sys)
 prob = ODEProblem(sys)
 sol = solve(prob, Tsit5())
 @test sol.t[end] == tspan[end]
@@ -448,6 +450,7 @@ prob = ODEProblem{false}(sys; u0_constructor = x -> SVector(x...))
 D = Differential(t)
 eqs = [D(x1) ~ -x1]
 @named sys = ODESystem(eqs, t, [x1, x2], [])
+sys = complete(sys)
 @test_throws ArgumentError ODEProblem(sys, [1.0, 1.0], (0.0, 1.0))
 @test_nowarn ODEProblem(sys, [1.0, 1.0], (0.0, 1.0), check_length = false)
 
@@ -545,6 +548,7 @@ eqs = [D(x) ~ foo(x, ms); D.(ms) .~ 1]
 @named sys = ODESystem(eqs, t, [x; ms], [])
 @named emptysys = ODESystem(Equation[], t)
 @named outersys = compose(emptysys, sys)
+outersys = complete(outersys)
 prob = ODEProblem(outersys, [sys.x => 1.0; collect(sys.ms) .=> 1:3], (0, 1.0))
 @test_nowarn solve(prob, Tsit5())
 
@@ -631,6 +635,7 @@ eqs[end] = D(D(z)) ~ α * x - β * y
     end
 
     @named sys = ODESystem(eqs, t, us, ps; defaults = defs, preface = preface)
+    sys = complete(sys)
     prob = ODEProblem(sys, [], (0.0, 1.0))
     sol = solve(prob, Euler(); dt = 0.1)
 
@@ -690,6 +695,7 @@ let
     D = Differential(t)
     eqs = [D(A) ~ -k1 * k2 * A]
     @named sys = ODESystem(eqs, t)
+    sys = complete(sys)
     u0map = [A => 1.0]
     pmap = (k1 => 1.0, k2 => 1)
     tspan = (0.0, 1.0)
@@ -843,6 +849,7 @@ let
     D = Differential(t)
     eqs = [D(A) ~ -k * A]
     @named osys = ODESystem(eqs, t)
+    osys = complete(osys)
     oprob = ODEProblem(osys, [A => 1.0], (0.0, 10.0), [k => 1.0]; check_length = false)
     @test_nowarn sol = solve(oprob, Tsit5())
 end
@@ -958,7 +965,7 @@ testdict = Dict([:name => "test"])
 eqs = [∂t(Q) ~ 1 / sin(P)
     ∂t(P) ~ log(-cos(Q))]
 @named sys = ODESystem(eqs, t, [P, Q], [])
-sys = debug_system(sys);
+sys = complete(debug_system(sys));
 prob = ODEProblem(sys, [], (0, 1.0));
 du = zero(prob.u0);
 if VERSION < v"1.8"
diff --git a/test/optimizationsystem.jl b/test/optimizationsystem.jl
index 4e57e633a6..5d9b462d32 100644
--- a/test/optimizationsystem.jl
+++ b/test/optimizationsystem.jl
@@ -14,8 +14,8 @@ using ModelingToolkit: get_metadata
     @variables z
     @parameters β
     loss2 = sys1.x - sys2.y + z * β
-    combinedsys = OptimizationSystem(loss2, [z], [β], systems = [sys1, sys2],
-        name = :combinedsys)
+    combinedsys = complete(OptimizationSystem(loss2, [z], [β], systems = [sys1, sys2],
+        name = :combinedsys))
 
     equations(combinedsys)
     unknowns(combinedsys)
@@ -27,7 +27,7 @@ using ModelingToolkit: get_metadata
     generate_gradient(combinedsys)
     generate_hessian(combinedsys)
     hess_sparsity = ModelingToolkit.hessian_sparsity(sys1)
-    sparse_prob = OptimizationProblem(sys1, [x, y], [a, b], grad = true, sparse = true)
+    sparse_prob = OptimizationProblem(complete(sys1), [x, y], [a, b], grad = true, sparse = true)
     @test sparse_prob.f.hess_prototype.rowval == hess_sparsity.rowval
     @test sparse_prob.f.hess_prototype.colptr == hess_sparsity.colptr
 
@@ -57,7 +57,7 @@ end
         x^2 + y^2 ≲ 1.0,
     ]
     @named sys = OptimizationSystem(loss, [x, y], [a, b], constraints = cons)
-
+    sys = complete(sys)
     prob = OptimizationProblem(sys, [x => 0.0, y => 0.0], [a => 1.0, b => 1.0],
         grad = true, hess = true, cons_j = true, cons_h = true)
     @test prob.f.sys === sys
@@ -149,8 +149,8 @@ end
     ]
     sys1 = OptimizationSystem(o1, [x], [a], name = :sys1, constraints = c1)
     sys2 = OptimizationSystem(o2, [y], [], name = :sys2, constraints = c2)
-    sys = OptimizationSystem(0, [], []; name = :sys, systems = [sys1, sys2],
-        constraints = [sys1.x + sys2.y ~ 2], checks = false)
+    sys = complete(OptimizationSystem(0, [], []; name = :sys, systems = [sys1, sys2],
+        constraints = [sys1.x + sys2.y ~ 2], checks = false))
     prob = OptimizationProblem(sys, [0.0, 0.0])
     @test isequal(constraints(sys), vcat(sys1.x + sys2.y ~ 2, sys1.x ~ 1, sys2.y ~ 1))
     @test isequal(equations(sys), (sys1.x - sys1.a)^2 + (sys2.y - 1 / 2)^2)
@@ -190,8 +190,8 @@ end
     @variables z
     @parameters β
     loss2 = sys1.x - sys2.y + z * β
-    combinedsys = OptimizationSystem(loss2, [z], [β], systems = [sys1, sys2],
-        name = :combinedsys)
+    combinedsys = complete(OptimizationSystem(loss2, [z], [β], systems = [sys1, sys2],
+        name = :combinedsys))
 
     u0 = [sys1.x => 1.0
         sys1.y => 2.0
@@ -243,7 +243,7 @@ end
             x[1]^2 + x[2]^2 ≲ 2.0,
         ])
 
-    prob = OptimizationProblem(sys, [x[1] => 2.0, x[2] => 0.0], [], grad = true,
+    prob = OptimizationProblem(complete(sys), [x[1] => 2.0, x[2] => 0.0], [], grad = true,
         hess = true, cons_j = true, cons_h = true)
     sol = Optimization.solve(prob, Ipopt.Optimizer(); print_level = 0)
     @test sol.u ≈ [1, 0]
@@ -257,7 +257,7 @@ end
     @parameters a b
     loss = (a - x)^2 + b * (y - x^2)^2
     @named sys = OptimizationSystem(loss, [x, y], [a, b, c])
-    prob = OptimizationProblem(sys, [x => 0.0, y => 0.0], [a => 1.0, b => 100.0])
+    prob = OptimizationProblem(complete(sys), [x => 0.0, y => 0.0], [a => 1.0, b => 100.0])
     @test prob.lb == [-Inf, 0.0]
     @test prob.ub == [Inf, Inf]
 end
@@ -269,12 +269,12 @@ end
     cons = [
         x₁^2 + x₂^2 ≲ 1.0,
     ]
-    sys1 = OptimizationSystem(loss, [x₁, x₂], [α₁, α₂], name = :sys1, constraints = cons)
+    sys1 = complete(OptimizationSystem(loss, [x₁, x₂], [α₁, α₂], name = :sys1, constraints = cons))
 
     prob1 = OptimizationProblem(sys1, [x₁ => 0.0, x₂ => 0.0], [α₁ => 1.0, α₂ => 100.0],
         grad = true, hess = true, cons_j = true, cons_h = true)
 
-    sys2 = modelingtoolkitize(prob1)
+    sys2 = complete(modelingtoolkitize(prob1))
     prob2 = OptimizationProblem(sys2, [x₁ => 0.0, x₂ => 0.0], [α₁ => 1.0, α₂ => 100.0],
         grad = true, hess = true, cons_j = true, cons_h = true)
 
@@ -289,6 +289,7 @@ end
     @parameters a b
     loss = (a - x)^2 + b * (y - x^2)^2
     @named sys = OptimizationSystem(loss, [x, y], [a, b], constraints = [x^2 + y^2 ≲ 0.0])
+    sys = complete(sys)
     @test_throws ArgumentError OptimizationProblem(sys,
         [x => 0.0, y => 0.0],
         [a => 1.0, b => 100.0],
diff --git a/test/sdesystem.jl b/test/sdesystem.jl
index a630deac99..6c92fb0cb5 100644
--- a/test/sdesystem.jl
+++ b/test/sdesystem.jl
@@ -21,6 +21,7 @@ noiseeqs = [0.1 * x,
 @test SDESystem(sys, noiseeqs, name = :foo) isa SDESystem
 
 @named de = SDESystem(eqs, noiseeqs, t, [x, y, z], [σ, ρ, β], tspan = (0.0, 10.0))
+de = complete(de)
 f = eval(generate_diffusion_function(de)[1])
 @test f(ones(3), rand(3), nothing) == 0.1ones(3)
 
@@ -42,6 +43,7 @@ noiseeqs_nd = [0.01*x 0.01*x*y 0.02*x*z
     σ 0.01*y 0.02*x*z
     ρ β 0.01*z]
 @named de = SDESystem(eqs, noiseeqs_nd, t, [x, y, z], [σ, ρ, β])
+de = complete(de)
 f = eval(generate_diffusion_function(de)[1])
 @test f([1, 2, 3.0], [0.1, 0.2, 0.3], nothing) == [0.01*1 0.01*1*2 0.02*1*3
     0.1 0.01*2 0.02*1*3
@@ -504,13 +506,14 @@ noiseeqs = [0.1 * x]
     ]
 
     @named de = SDESystem(eqs, noiseeqs, t, [x], [α, β], observed = [weight ~ x * 10])
-
+    de = complete(de)
     prob = SDEProblem(de, u0map, (0.0, 1.0), parammap)
     sol = solve(prob, EM(), dt = dt)
     @test observed(de) == [weight ~ x * 10]
     @test sol[weight] == 10 * sol[x]
 
     @named ode = ODESystem(eqs, t, [x], [α, β], observed = [weight ~ x * 10])
+    ode = complete(ode)
     odeprob = ODEProblem(ode, u0map, (0.0, 1.0), parammap)
     solode = solve(odeprob, Tsit5())
     @test observed(ode) == [weight ~ x * 10]
@@ -528,7 +531,7 @@ end
     noiseeqs = [β * x]
 
     @named de = SDESystem(eqs, noiseeqs, t, [x], [α, β])
-
+    de = complete(de)
     g(x) = x[1]^2
     dt = 1 // 2^(7)
     x0 = 0.1
@@ -564,7 +567,7 @@ end
 
     ## Variance reduction method
     u = x
-    demod = ModelingToolkit.Girsanov_transform(de, u; θ0 = 0.1)
+    demod = complete(ModelingToolkit.Girsanov_transform(de, u; θ0 = 0.1))
 
     probmod = SDEProblem(demod, u0map, (0.0, 1.0), parammap)
 
@@ -609,6 +612,7 @@ diffusion_eqs = [s*x 0
     (s * x * z)-s * z 0]
 
 sys2 = SDESystem(drift_eqs, diffusion_eqs, t, sts, ps, name = :sys1)
+sys2 = complete(sys2)
 @test sys1 == sys2
 
 prob = SDEProblem(sys1, sts .=> [1.0, 0.0, 0.0],
diff --git a/test/serialization.jl b/test/serialization.jl
index 79f6f34bdf..430a8765d1 100644
--- a/test/serialization.jl
+++ b/test/serialization.jl
@@ -5,6 +5,7 @@ using ModelingToolkit, SciMLBase, Serialization, OrdinaryDiffEq
 D = Differential(t)
 
 @named sys = ODESystem([D(x) ~ -0.5 * x], defaults = Dict(x => 1.0))
+sys = complete(sys)
 for prob in [
     eval(ModelingToolkit.ODEProblem{false}(sys, nothing, nothing,
         SciMLBase.NullParameters())),
@@ -37,7 +38,7 @@ sol = solve(prob, ImplicitEuler())
 
 ## Check ODESystem with Observables ----------
 ss_exp = ModelingToolkit.toexpr(ss)
-ss_ = eval(ss_exp)
+ss_ = complete(eval(ss_exp))
 prob_ = ODEProblem(ss_, [], (0, 0.1))
 sol_ = solve(prob_, ImplicitEuler())
 @test sol[all_obs] == sol_[all_obs]
@@ -64,5 +65,5 @@ end)
 probexpr = ODEProblemExpr{true}(ss, [], (0, 0.1); observedfun_exp);
 prob_obs = eval(probexpr)
 sol_obs = solve(prob_obs, ImplicitEuler())
-
+@show all_obs
 @test sol_obs[all_obs] == sol[all_obs]
diff --git a/test/steadystatesystems.jl b/test/steadystatesystems.jl
index 9d31844b23..9b5a6c603f 100644
--- a/test/steadystatesystems.jl
+++ b/test/steadystatesystems.jl
@@ -7,6 +7,7 @@ using Test
 D = Differential(t)
 eqs = [D(x) ~ x^2 - r]
 @named de = ODESystem(eqs)
+de = complete(de)
 
 for factor in [1e-1, 1e0, 1e10],
     u0_p in [(2.34, 2.676), (22.34, 1.632), (0.3, 15.676), (0.3, 0.006)]
diff --git a/test/structural_transformation/index_reduction.jl b/test/structural_transformation/index_reduction.jl
index a56a04f4c4..40236b9d94 100644
--- a/test/structural_transformation/index_reduction.jl
+++ b/test/structural_transformation/index_reduction.jl
@@ -71,7 +71,7 @@ prob = ODEProblem(ODEFunction(first_order_idx1_pendulum),
 sol = solve(prob, Rodas5());
 #plot(sol, idxs=(1, 2))
 
-new_sys = dae_index_lowering(ModelingToolkit.ode_order_lowering(pendulum2))
+new_sys = complete(dae_index_lowering(ModelingToolkit.ode_order_lowering(pendulum2)))
 
 prob_auto = ODEProblem(new_sys,
     [D(x) => 0,
@@ -98,7 +98,7 @@ pendulum2 = ODESystem(eqs2, t, [x, y, T], [L, g], name = :pendulum)
 first_order_sys = ModelingToolkit.ode_order_lowering(pendulum2)
 
 # Perform index reduction to get an Index 1 DAE
-new_sys = dae_index_lowering(first_order_sys)
+new_sys = complete(dae_index_lowering(first_order_sys))
 
 u0 = [
     D(x) => 0.0,
diff --git a/test/structural_transformation/tearing.jl b/test/structural_transformation/tearing.jl
index 818daf6c95..b0152828ba 100644
--- a/test/structural_transformation/tearing.jl
+++ b/test/structural_transformation/tearing.jl
@@ -119,7 +119,7 @@ end
 # unknowns: u5
 # solve for
 # 0 = u5 - hypot(sin(u5), hypot(cos(sin(u5)), hypot(sin(u5), cos(sin(u5)))))
-tornsys = tearing(sys)
+tornsys = complete(tearing(sys))
 @test isequal(equations(tornsys), [0 ~ u5 - hypot(u4, u1)])
 prob = NonlinearProblem(tornsys, ones(1))
 sol = solve(prob, NewtonRaphson())
diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl
index d85f2118fd..5ac0604774 100644
--- a/test/symbolic_events.jl
+++ b/test/symbolic_events.jl
@@ -136,6 +136,8 @@ fsys = flatten(sys)
 @test isequal(ModelingToolkit.continuous_events(sys2)[1].eqs[], x ~ 2)
 @test isequal(ModelingToolkit.continuous_events(sys2)[2].eqs[], sys.x ~ 1)
 
+sys = complete(sys)
+sys2 = complete(sys2)
 # Functions should be generated for root-finding equations
 prob = ODEProblem(sys, Pair[], (0.0, 2.0))
 p0 = 0
@@ -190,6 +192,7 @@ sol = solve(prob, Tsit5())
 @test minimum(t -> abs(t - 2), sol.t) < 1e-10 # test that the solver stepped at the second root
 
 @named sys = ODESystem(eqs, continuous_events = [x ~ 1, x ~ 2]) # two root eqs using the same unknown
+sys = complete(sys)
 prob = ODEProblem(sys, Pair[], (0.0, 3.0))
 @test get_callback(prob) isa ModelingToolkit.DiffEqCallbacks.VectorContinuousCallback
 sol = solve(prob, Tsit5())
@@ -341,7 +344,7 @@ sys = structural_simplify(model)
 let
     function testsol(osys, u0, p, tspan; tstops = Float64[], skipparamtest = false,
             kwargs...)
-        oprob = ODEProblem(osys, u0, tspan, p; kwargs...)
+        oprob = ODEProblem(complete(osys), u0, tspan, p; kwargs...)
         sol = solve(oprob, Tsit5(); tstops = tstops, abstol = 1e-10, reltol = 1e-10)
         @test isapprox(sol(1.0000000001)[1] - sol(0.999999999)[1], 1.0; rtol = 1e-6)
         !skipparamtest && (@test oprob.p[1] == 1.0)
@@ -392,7 +395,7 @@ let
     end
     cb2‵‵ = [2.0] => (affect!, [], [k], nothing)
     @named osys4 = ODESystem(eqs, t, [A], [k, t1], discrete_events = [cb1, cb2‵‵])
-    oprob4 = ODEProblem(osys4, u0, tspan, p)
+    oprob4 = ODEProblem(complete(osys4), u0, tspan, p)
     testsol(osys4, u0, p, tspan; tstops = [1.0])
 
     # mixing with symbolic condition in the func affect
@@ -415,7 +418,7 @@ end
 let
     function testsol(ssys, u0, p, tspan; tstops = Float64[], skipparamtest = false,
             kwargs...)
-        sprob = SDEProblem(ssys, u0, tspan, p; kwargs...)
+        sprob = SDEProblem(complete(ssys), u0, tspan, p; kwargs...)
         sol = solve(sprob, RI5(); tstops = tstops, abstol = 1e-10, reltol = 1e-10)
         @test isapprox(sol(1.0000000001)[1] - sol(0.999999999)[1], 1.0; rtol = 1e-4)
         !skipparamtest && (@test sprob.p[1] == 1.0)
@@ -495,6 +498,7 @@ end
 let rng = rng
     function testsol(jsys, u0, p, tspan; tstops = Float64[], skipparamtest = false,
             N = 40000, kwargs...)
+        jsys = complete(jsys)
         dprob = DiscreteProblem(jsys, u0, tspan, p)
         jprob = JumpProblem(jsys, dprob, Direct(); kwargs...)
         sol = solve(jprob, SSAStepper(); tstops = tstops)
diff --git a/test/symbolic_parameters.jl b/test/symbolic_parameters.jl
index d779676b49..4e4df0e063 100644
--- a/test/symbolic_parameters.jl
+++ b/test/symbolic_parameters.jl
@@ -25,7 +25,7 @@ resolved = ModelingToolkit.varmap_to_vars(Dict(), parameters(ns),
     defaults = ModelingToolkit.defaults(ns))
 @test resolved == [1, 0.1 + 1, (0.1 + 1) * 1.1]
 
-prob = NonlinearProblem(ns, [u => 1.0], Pair[])
+prob = NonlinearProblem(complete(ns), [u => 1.0], Pair[])
 @test prob.u0 == [1.0, 1.1, 0.9]
 @show sol = solve(prob, NewtonRaphson())
 
@@ -39,6 +39,7 @@ res = ModelingToolkit.varmap_to_vars(Dict(), parameters(top),
     defaults = ModelingToolkit.defaults(top))
 @test res == [0.5, 1, 0.1 + 1, (0.1 + 1) * 1.1]
 
+top = complete(top)
 prob = NonlinearProblem(top, [unknowns(ns, u) => 1.0, a => 1.0], [])
 @test prob.u0 == [1.0, 0.5, 1.1, 0.9]
 @show sol = solve(prob, NewtonRaphson())
@@ -59,6 +60,7 @@ end)
 der = Differential(t)
 eqs = [der(x) ~ x]
 @named sys = ODESystem(eqs, t, vars, [x0])
+sys = complete(sys)
 pars = [
     x0 => 10.0,
 ]

From a779f4b6a1febc4bf65935ac26cde2eae1d5124f Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Wed, 31 Jan 2024 13:15:29 +0530
Subject: [PATCH 9/9] refactor!: require systems to be completed before
 creating an `XFunction`/`XFunctionExpr`

---
 src/systems/diffeqs/abstractodesystem.jl       | 18 ++++++++++++++++++
 src/systems/diffeqs/sdesystem.jl               |  6 ++++++
 src/systems/nonlinear/nonlinearsystem.jl       |  6 ++++++
 test/distributed.jl                            |  1 +
 test/function_registration.jl                  |  4 ++++
 test/labelledarrays.jl                         |  1 +
 test/mass_matrix.jl                            |  1 +
 test/odesystem.jl                              | 11 +++++++----
 test/precompile_test/ODEPrecompileTest.jl      |  1 +
 .../index_reduction.jl                         |  2 +-
 10 files changed, 46 insertions(+), 5 deletions(-)

diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 8b55e0ee86..da151c9b0f 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -308,6 +308,9 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, dvs = u
         analytic = nothing,
         split_idxs = nothing,
         kwargs...) where {iip, specialize}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEFunction`")
+    end
     f_gen = generate_function(sys, dvs, ps; expression = Val{eval_expression},
         expression_module = eval_module, checkbounds = checkbounds,
         kwargs...)
@@ -504,6 +507,9 @@ function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys)
         eval_module = @__MODULE__,
         checkbounds = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `DAEFunction`")
+    end
     f_gen = generate_function(sys, dvs, ps; implicit_dae = true,
         expression = Val{eval_expression},
         expression_module = eval_module, checkbounds = checkbounds,
@@ -579,6 +585,9 @@ function DiffEqBase.DDEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys)
         eval_module = @__MODULE__,
         checkbounds = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `DDEFunction`")
+    end
     f_gen = generate_function(sys, dvs, ps; isdde = true,
         expression = Val{true},
         expression_module = eval_module, checkbounds = checkbounds,
@@ -603,6 +612,9 @@ function DiffEqBase.SDDEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys
         eval_module = @__MODULE__,
         checkbounds = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `SDDEFunction`")
+    end
     f_gen = generate_function(sys, dvs, ps; isdde = true,
         expression = Val{true},
         expression_module = eval_module, checkbounds = checkbounds,
@@ -656,6 +668,9 @@ function ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = unknowns(sys),
         sparsity = false,
         observedfun_exp = nothing,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEFunctionExpr`")
+    end
     f_oop, f_iip = generate_function(sys, dvs, ps; expression = Val{true}, kwargs...)
 
     dict = Dict()
@@ -830,6 +845,9 @@ function DAEFunctionExpr{iip}(sys::AbstractODESystem, dvs = unknowns(sys),
         linenumbers = false,
         sparse = false, simplify = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `DAEFunctionExpr`")
+    end
     f_oop, f_iip = generate_function(sys, dvs, ps; expression = Val{true},
         implicit_dae = true, kwargs...)
     fsym = gensym(:f)
diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl
index 8a845ab0fe..7e11a33f52 100644
--- a/src/systems/diffeqs/sdesystem.jl
+++ b/src/systems/diffeqs/sdesystem.jl
@@ -393,6 +393,9 @@ function DiffEqBase.SDEFunction{iip}(sys::SDESystem, dvs = unknowns(sys),
         jac = false, Wfact = false, eval_expression = true,
         checkbounds = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `SDESystem` is required. Call `complete` or `structural_simplify` on the system before creating an `SDEFunction`")
+    end
     dvs = scalarize.(dvs)
     ps = scalarize.(ps)
 
@@ -515,6 +518,9 @@ function SDEFunctionExpr{iip}(sys::SDESystem, dvs = unknowns(sys),
         jac = false, Wfact = false,
         sparse = false, linenumbers = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `SDESystem` is required. Call `complete` or `structural_simplify` on the system before creating an `SDEFunctionExpr`")
+    end
     idx = iip ? 2 : 1
     f = generate_function(sys, dvs, ps; expression = Val{true}, kwargs...)[idx]
     g = generate_diffusion_function(sys, dvs, ps; expression = Val{true}, kwargs...)[idx]
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index 4f94e8971c..e25fa2ffaf 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -234,6 +234,9 @@ function SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(s
         eval_expression = true,
         sparse = false, simplify = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `NonlinearFunction`")
+    end
     f_gen = generate_function(sys, dvs, ps; expression = Val{eval_expression}, kwargs...)
     f_oop, f_iip = eval_expression ?
                    (drop_expr(@RuntimeGeneratedFunction(ex)) for ex in f_gen) : f_gen
@@ -296,6 +299,9 @@ function NonlinearFunctionExpr{iip}(sys::NonlinearSystem, dvs = unknowns(sys),
         linenumbers = false,
         sparse = false, simplify = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `NonlinearFunctionExpr`")
+    end
     idx = iip ? 2 : 1
     f = generate_function(sys, dvs, ps; expression = Val{true}, kwargs...)[idx]
 
diff --git a/test/distributed.jl b/test/distributed.jl
index fb89dfda1e..ef30a1b96f 100644
--- a/test/distributed.jl
+++ b/test/distributed.jl
@@ -14,6 +14,7 @@ addprocs(2)
     D(z) ~ x * y - β * z]
 
 @everywhere @named de = ODESystem(eqs)
+@everywhere de = complete(de)
 @everywhere ode_func = ODEFunction(de, [x, y, z], [σ, ρ, β])
 
 @everywhere u0 = [19.0, 20.0, 50.0]
diff --git a/test/function_registration.jl b/test/function_registration.jl
index 183d71343b..efa7764914 100644
--- a/test/function_registration.jl
+++ b/test/function_registration.jl
@@ -18,6 +18,7 @@ end
 
 eq = Dt(u) ~ do_something(x) + MyModule.do_something(x)
 @named sys = ODESystem([eq], t, [u], [x])
+sys = complete(sys)
 fun = ODEFunction(sys)
 
 u0 = 5.0
@@ -40,6 +41,7 @@ end
 
 eq = Dt(u) ~ do_something_2(x) + MyNestedModule.do_something_2(x)
 @named sys = ODESystem([eq], t, [u], [x])
+sys = complete(sys)
 fun = ODEFunction(sys)
 
 u0 = 3.0
@@ -61,6 +63,7 @@ end
 
 eq = Dt(u) ~ do_something_3(x) + (@__MODULE__).do_something_3(x)
 @named sys = ODESystem([eq], t, [u], [x])
+sys = complete(sys)
 fun = ODEFunction(sys)
 
 u0 = 7.0
@@ -99,6 +102,7 @@ function build_ode()
     Dt = Differential(t)
     eq = Dt(u) ~ do_something_4(x) + (@__MODULE__).do_something_4(x)
     @named sys = ODESystem([eq], t, [u], [x])
+    sys = complete(sys)
     fun = ODEFunction(sys, eval_expression = false)
 end
 function run_test()
diff --git a/test/labelledarrays.jl b/test/labelledarrays.jl
index 62890a1dfa..c47b87c1e2 100644
--- a/test/labelledarrays.jl
+++ b/test/labelledarrays.jl
@@ -13,6 +13,7 @@ eqs = [D(x) ~ σ * (y - x),
     D(z) ~ x * y - β * z]
 
 @named de = ODESystem(eqs)
+de = complete(de)
 ff = ODEFunction(de, [x, y, z], [σ, ρ, β], jac = true)
 
 a = @SVector [1.0, 2.0, 3.0]
diff --git a/test/mass_matrix.jl b/test/mass_matrix.jl
index 1e54f58ef1..0c92e6f328 100644
--- a/test/mass_matrix.jl
+++ b/test/mass_matrix.jl
@@ -9,6 +9,7 @@ eqs = [D(y[1]) ~ -k[1] * y[1] + k[3] * y[2] * y[3],
     0 ~ y[1] + y[2] + y[3] - 1]
 
 @named sys = ODESystem(eqs, t, y, k)
+sys = complete(sys)
 @test_throws ArgumentError ODESystem(eqs, y[1])
 M = calculate_massmatrix(sys)
 @test M == [1 0 0
diff --git a/test/odesystem.jl b/test/odesystem.jl
index 5ccbaec536..d7ed5a41ad 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -51,6 +51,7 @@ jac_expr = generate_jacobian(de)
 jac = calculate_jacobian(de)
 jacfun = eval(jac_expr[2])
 
+de = complete(de)
 for f in [
     ODEFunction(de, [x, y, z], [σ, ρ, β], tgrad = true, jac = true),
     eval(ODEFunctionExpr(de, [x, y, z], [σ, ρ, β], tgrad = true, jac = true)),
@@ -167,7 +168,7 @@ lowered_eqs = [D(uˍtt) ~ 2uˍtt + uˍt + xˍt + 1
 
 test_diffeq_inference("first-order transform", de1, t, [uˍtt, xˍt, uˍt, u, x], [])
 du = zeros(5)
-ODEFunction(de1, [uˍtt, xˍt, uˍt, u, x], [])(du, ones(5), nothing, 0.1)
+ODEFunction(complete(de1), [uˍtt, xˍt, uˍt, u, x], [])(du, ones(5), nothing, 0.1)
 @test du == [5.0, 3.0, 1.0, 1.0, 1.0]
 
 # Internal calculations
@@ -182,7 +183,7 @@ jac = calculate_jacobian(de)
 @test ModelingToolkit.jacobian_sparsity(de).colptr == sparse(jac).colptr
 @test ModelingToolkit.jacobian_sparsity(de).rowval == sparse(jac).rowval
 
-f = ODEFunction(de, [x, y, z], [σ, ρ, β])
+f = ODEFunction(complete(de), [x, y, z], [σ, ρ, β])
 
 D = Differential(t)
 @parameters A B C
@@ -208,7 +209,7 @@ function lotka(u, p, t)
 end
 
 prob = ODEProblem(ODEFunction{false}(lotka), [1.0, 1.0], (0.0, 1.0), [1.5, 1.0, 3.0, 1.0])
-de = modelingtoolkitize(prob)
+de = complete(modelingtoolkitize(prob))
 ODEFunction(de)(similar(prob.u0), prob.u0, prob.p, 0.1)
 
 function lotka(du, u, p, t)
@@ -220,7 +221,7 @@ end
 
 prob = ODEProblem(lotka, [1.0, 1.0], (0.0, 1.0), [1.5, 1.0, 3.0, 1.0])
 
-de = modelingtoolkitize(prob)
+de = complete(modelingtoolkitize(prob))
 ODEFunction(de)(similar(prob.u0), prob.u0, prob.p, 0.1)
 
 # automatic unknown detection for DAEs
@@ -579,11 +580,13 @@ eqs = [
 ]
 
 @named sys = ODESystem(eqs, t, [x, y, z], [α, β])
+sys = complete(sys)
 @test_throws Any ODEFunction(sys)
 
 eqs = copy(eqs)
 eqs[end] = D(D(z)) ~ α * x - β * y
 @named sys = ODESystem(eqs, t, [x, y, z], [α, β])
+sys = complete(sys)
 @test_throws Any ODEFunction(sys)
 
 @testset "Preface tests" begin
diff --git a/test/precompile_test/ODEPrecompileTest.jl b/test/precompile_test/ODEPrecompileTest.jl
index 4a8eb64d5c..caf411f88d 100644
--- a/test/precompile_test/ODEPrecompileTest.jl
+++ b/test/precompile_test/ODEPrecompileTest.jl
@@ -13,6 +13,7 @@ function system(; kwargs...)
         D(z) ~ x * y - β * z]
 
     @named de = ODESystem(eqs)
+    de = complete(de)
     return ODEFunction(de, [x, y, z], [σ, ρ, β]; kwargs...)
 end
 
diff --git a/test/structural_transformation/index_reduction.jl b/test/structural_transformation/index_reduction.jl
index 40236b9d94..f51e9d5060 100644
--- a/test/structural_transformation/index_reduction.jl
+++ b/test/structural_transformation/index_reduction.jl
@@ -59,7 +59,7 @@ idx1_pendulum = [D(x) ~ w,
     # substitute the rhs
     0 ~ 2x * (T * x) + 2 * xˍt * xˍt + 2y * (T * y - g) + 2 * yˍt * yˍt]
 @named idx1_pendulum = ODESystem(idx1_pendulum, t, [x, y, w, z, xˍt, yˍt, T], [L, g])
-first_order_idx1_pendulum = ode_order_lowering(idx1_pendulum)
+first_order_idx1_pendulum = complete(ode_order_lowering(idx1_pendulum))
 
 using OrdinaryDiffEq
 using LinearAlgebra