SciML · AayushSabharwal · Oct 6, 2025 · Oct 8, 2025 · Oct 8, 2025 · Oct 9, 2025
diff --git a/Project.toml b/Project.toml
@@ -8,7 +8,6 @@ ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
@@ -33,7 +32,6 @@ FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 ImplicitDiscreteSolve = "3263718b-31ed-49cf-8a0f-35a466e8af96"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
@@ -68,18 +66,22 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [weakdeps]
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 CasADi = "c49709b8-5c63-11e9-2fb2-69db5844192f"
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
 FMI = "14a09403-18e3-468f-ad8a-74f8dda2d9ac"
 InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57"
+JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
 Pyomo = "0e8e1daf-01b5-4eba-a626-3897743a3816"
 
 [extensions]
 MTKBifurcationKitExt = "BifurcationKit"
 MTKCasADiDynamicOptExt = "CasADi"
+MTKChainRulesCoreExt = "ChainRulesCore"
 MTKDeepDiffsExt = "DeepDiffs"
 MTKFMIExt = "FMI"
 MTKInfiniteOptExt = "InfiniteOpt"
+MTKJuliaFormatterExt = "JuliaFormatter"
 MTKLabelledArraysExt = "LabelledArrays"
 MTKPyomoDynamicOptExt = "Pyomo"
 

diff --git a/docs/Project.toml b/docs/Project.toml
@@ -32,7 +32,6 @@ StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
-Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 Attractors = "1.24"
@@ -65,4 +64,3 @@ StochasticDiffEq = "6"
 SymbolicIndexingInterface = "0.3.1"
 SymbolicUtils = "3, 4"
 Symbolics = "6"
-Unitful = "1.12"
diff --git a/docs/src/API/codegen.md b/docs/src/API/codegen.md
@@ -50,5 +50,5 @@ ModelingToolkit.calculate_A_b
 All code generation eventually calls `build_function_wrapper`.
 
 ```@docs
-build_function_wrapper
+ModelingToolkit.build_function_wrapper
 ```
diff --git a/docs/src/API/model_building.md b/docs/src/API/model_building.md
@@ -13,8 +13,6 @@ ModelingToolkit.t_nounits
 ModelingToolkit.D_nounits
 ModelingToolkit.t
 ModelingToolkit.D
-ModelingToolkit.t_unitful
-ModelingToolkit.D_unitful
 ```
 
 Users are recommended to use the appropriate common definition in their models. The required

diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md
@@ -192,9 +192,7 @@ p, replace, alias = SciMLStructures.canonicalize(Tunable(), prob.p)
 # changes to the array will be reflected in parameter values
 ```
 
-See the [basic example on optimizing](https://docs.sciml.ai/ModelingToolkit/dev/examples/remake/#Optimizing-through-an-ODE-solve-and-re-creating-MTK-Problems) for combining these steps to optimizing parameters and use ForwardDiff.jl as the backend for Automatic Differentiation.
-
-# ERROR: ArgumentError: SymbolicUtils.BasicSymbolic{Real}[xˍt(t)] are missing from the variable map.
+# ERROR: ArgumentError: `[xˍt(t)]` are missing from the variable map.
 
 This error can come up after running `mtkcompile` on a system that generates dummy derivatives (i.e. variables with `ˍt`).  For example, here even though all the variables are defined with initial values, the `ODEProblem` generation will throw an error that defaults are missing from the variable map.
 

diff --git a/docs/src/basics/Validation.md b/docs/src/basics/Validation.md
@@ -155,7 +155,7 @@ future when `ModelingToolkit` is extended to support eliminating `DynamicQuantit
 
 ## Other Restrictions
 
-`Unitful` provides non-scalar units such as `dBm`, `°C`, etc. At this time, `ModelingToolkit` only supports scalar quantities. Additionally, angular degrees (`°`) are not supported because trigonometric functions will treat plain numerical values as radians, which would lead systems validated using degrees to behave erroneously when being solved.
+`DynamicQuantities` provides non-scalar units such as `°C`, etc. At this time, `ModelingToolkit` only supports scalar quantities. Additionally, angular degrees (`°`) are not supported because trigonometric functions will treat plain numerical values as radians, which would lead systems validated using degrees to behave erroneously when being solved.
 
 ## Troubleshooting & Gotchas
 
@@ -169,7 +169,7 @@ Parameter and initial condition values are supplied to problem constructors as p
 
 ```julia
 function remove_units(p::Dict)
-    Dict(k => Unitful.ustrip(ModelingToolkit.get_unit(k), v) for (k, v) in p)
+    Dict(k => DynamicQuantities.ustrip(ModelingToolkit.get_unit(k), v) for (k, v) in p)
 end
 add_units(p::Dict) = Dict(k => v * ModelingToolkit.get_unit(k) for (k, v) in p)
 ```

diff --git a/ext/MTKCasADiDynamicOptExt.jl b/ext/MTKCasADiDynamicOptExt.jl
@@ -122,7 +122,7 @@ end
 
 function MTK.lowered_var(m::CasADiModel, uv, i, t)
     X = getfield(m, uv)
-    t isa Union{Num, Symbolics.Symbolic} ? X.u[i, :] : X(t)[i]
+    t isa Union{Num, SymbolicT} ? X.u[i, :] : X(t)[i]
 end
 
 function MTK.lowered_integral(model::CasADiModel, expr, lo, hi)

diff --git a/src/adjoints.jl → ext/MTKChainRulesCoreExt.jl b/src/adjoints.jl → ext/MTKChainRulesCoreExt.jl
@@ -1,3 +1,14 @@
+module MTKChainRulesCoreExt
+
+import ChainRulesCore
+import ChainRulesCore: Tangent, ZeroTangent, NoTangent, zero_tangent, unthunk
+using ModelingToolkit: MTKParameters, NONNUMERIC_PORTION, AbstractSystem
+import ModelingToolkit
+import ModelingToolkit as MTK
+import SciMLStructures
+import SymbolicIndexingInterface: remake_buffer
+import SciMLBase: AbstractNonlinearProblem, remake
+
 function ChainRulesCore.rrule(::Type{MTKParameters}, tunables, args...)
     function mtp_pullback(dt)
         dt = unthunk(dt)
@@ -104,3 +115,11 @@ function ChainRulesCore.rrule(
 end
 
 ChainRulesCore.@non_differentiable Base.getproperty(sys::AbstractSystem, x::Symbol)
+
+function ModelingToolkit.update_initializeprob!(initprob::AbstractNonlinearProblem, prob)
+    pgetter = ChainRulesCore.@ignore_derivatives MTK.get_scimlfn(prob).initialization_data.metadata.oop_reconstruct_u0_p.pgetter
+    p = pgetter(prob, initprob)
+    return remake(initprob; p)
+end
+
+end
diff --git a/ext/MTKInfiniteOptExt.jl b/ext/MTKInfiniteOptExt.jl
@@ -122,7 +122,7 @@ end
 
 function MTK.lowered_var(m::InfiniteOptModel, uv, i, t)
     X = getfield(m, uv)
-    t isa Union{Num, Symbolics.Symbolic} ? X[i] : X[i](t)
+    t isa Union{Num, SymbolicT} ? X[i] : X[i](t)
 end
 
 function add_solve_constraints!(prob::JuMPDynamicOptProblem, tableau)
@@ -256,13 +256,13 @@ for ff in [acos, log1p, acosh, log2, asin, tan, atanh, cos, log, sin, log10, sqr
 end
 
 # JuMP variables and Symbolics variables never compare equal. When tracing through dynamics, a function argument can be either a JuMP variable or A Symbolics variable, it can never be both.
-function Base.isequal(::SymbolicUtils.Symbolic,
+function Base.isequal(::SymbolicT,
         ::Union{JuMP.GenericAffExpr, JuMP.GenericQuadExpr, InfiniteOpt.AbstractInfOptExpr})
     false
 end
 function Base.isequal(
         ::Union{JuMP.GenericAffExpr, JuMP.GenericQuadExpr, InfiniteOpt.AbstractInfOptExpr},
-        ::SymbolicUtils.Symbolic)
+        ::SymbolicT)
     false
 end
 end
diff --git a/ext/MTKJuliaFormatterExt.jl b/ext/MTKJuliaFormatterExt.jl
@@ -0,0 +1,12 @@
+module MTKJuliaFormatterExt
+
+import ModelingToolkit: readable_code, _readable_code, rec_remove_macro_linenums!
+import JuliaFormatter
+
+function readable_code(expr::Expr)
+    expr = Base.remove_linenums!(_readable_code(expr))
+    rec_remove_macro_linenums!(expr)
+    JuliaFormatter.format_text(string(expr), JuliaFormatter.SciMLStyle())
+end
+
+end
diff --git a/ext/MTKPyomoDynamicOptExt.jl b/ext/MTKPyomoDynamicOptExt.jl
@@ -5,6 +5,7 @@ using DiffEqBase
 using UnPack
 using NaNMath
 using Setfield
+import SymbolicUtils as SU
 const MTK = ModelingToolkit
 
 const SPECIAL_FUNCTIONS_DICT = Dict([acos => Pyomo.py_acos,
@@ -53,7 +54,7 @@ struct PyomoDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
     end
 end
 
-function pysym_getproperty(s::Union{Num, Symbolics.Symbolic}, name::Symbol)
+function pysym_getproperty(s::Union{Num, SymbolicT}, name::Symbol)
     Symbolics.wrap(SymbolicUtils.term(
         _getproperty, Symbolics.unwrap(s), Val{name}(), type = Symbolics.Struct{PyomoVar}))
 end
@@ -112,7 +113,7 @@ function MTK.add_constraint!(pmodel::PyomoDynamicOptModel, cons; n_idxs = 1)
         Symbolics.unwrap(expr), SPECIAL_FUNCTIONS_DICT, fold = false)
 
     cons_sym = Symbol("cons", hash(cons))
-    if occursin(Symbolics.unwrap(t_sym), expr)
+    if SU.query(isequal(Symbolics.unwrap(t_sym)), expr)
         f = eval(Symbolics.build_function(expr, model_sym, t_sym))
         setproperty!(model, cons_sym, pyomo.Constraint(model.t, rule = Pyomo.pyfunc(f)))
     else
@@ -124,7 +125,7 @@ end
 function MTK.set_objective!(pmodel::PyomoDynamicOptModel, expr)
     @unpack model, model_sym, t_sym, dummy_sym = pmodel
     expr = Symbolics.substitute(expr, SPECIAL_FUNCTIONS_DICT, fold = false)
-    if occursin(Symbolics.unwrap(t_sym), expr)
+    if SU.query(isequal(Symbolics.unwrap(t_sym)), expr)
         f = eval(Symbolics.build_function(expr, model_sym, t_sym))
         model.obj = pyomo.Objective(model.t, rule = Pyomo.pyfunc(f))
     else
@@ -165,7 +166,7 @@ end
 
 function MTK.lowered_var(m::PyomoDynamicOptModel, uv, i, t)
     X = Symbolics.value(pysym_getproperty(m.model_sym, uv))
-    var = t isa Union{Num, Symbolics.Symbolic} ? X[i, m.t_sym] : X[i, t]
+    var = t isa Union{Num, SymbolicT} ? X[i, m.t_sym] : X[i, t]
     Symbolics.unwrap(var)
 end