fix: promote the output array of initialization to the tunable eltype

This adresses the cases where the
we generate output arrays from
fully constant RHS, which would
end up being incompatible
with the elype of parameters
in ForwardDiff and other similar
contexts.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

Author: SebastianM-C

lib/ModelingToolkitBase/src/systems/problem_utils.jl

@@ -1159,6 +1159,41 @@ end
 safe_float(x) = x
 safe_float(x::AbstractArray) = isempty(x) ? x : float(x)
 
+"""
+    PromoteToTunableEltype(observed)
+
+Wraps an `initializeprob` observed function so its output array is promoted to an
+eltype compatible with the current tunable parameters. Addresses the case where
+the observed function is generated from fully constant RHS (e.g. `initialization_eqs
+= [s ~ 0]`): the resulting `create_array(Array, nothing, …, 0, 0)` would otherwise
+produce `Vector{Int64}`, which — when downstream `remake` reinstalls it as `u0` —
+silently defeats ForwardDiff/Tracker/Measurements promotion of `u0`.
+
+Replaces the previous `safe_float` layer by subsuming it: `promote_type(Int, Float64)
+== Float64`, so plain problems still get `Vector{Float64}`; `promote_type(Int,
+ForwardDiff.Dual)` yields the Dual type.
+"""
+struct PromoteToTunableEltype{F}
+    observed::F
+end
+
+function (p::PromoteToTunableEltype)(nlsol)
+    raw = p.observed(nlsol)
+    raw isa AbstractArray || return raw
+    isempty(raw) && return raw
+    T = promote_type(eltype(raw), _tunable_eltype(parameter_values(nlsol)), Float64)
+    T === eltype(raw) ? raw : convert(AbstractArray{T}, raw)
+end
+
+_tunable_eltype(p::MTKParameters) = isempty(p.tunable) ? Bool : eltype(p.tunable)
+function _tunable_eltype(p)
+    if SciMLStructures.isscimlstructure(p)
+        tun = SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]
+        return isempty(tun) ? Bool : eltype(tun)
+    end
+    return Bool
+end
+
 """
     $(TYPEDSIGNATURES)
 
@@ -1258,8 +1293,8 @@ function maybe_build_initialization_problem(
         if isempty(solved_unknowns)
             initializeprobmap = Returns(nothing)
         else
-            initializeprobmap = u0_constructor ∘ safe_float ∘
-                getu(initializeprob, solved_unknowns)
+            initializeprobmap = u0_constructor ∘ PromoteToTunableEltype(
+                getu(initializeprob, solved_unknowns))
         end
     else
         initializeprobmap = nothing

lib/ModelingToolkitBase/test/initializationsystem.jl

@@ -911,7 +911,7 @@ end
 
     @testset "No initialization for variables" begin
         @variables x = 1.0
-        @parameters p = 10.0 
+        @parameters p = 10.0
 
         eqs = [
             0 ~ x^2 + 2p * x + 3p
@@ -2013,3 +2013,33 @@ end
         @test !(prob.f.initialization_data.initializeprob isa SCCNonlinearProblem)
     end
 end
+
+@testset "Issue #4457" begin
+    @parameters m=1.5 d=9.0
+    @variables s(t) v(t)
+
+    eqs = [
+        D(s) ~ v
+        m * D(v) ~ 1 - d * v
+    ]
+
+    sys = mtkcompile(System(eqs, t;
+        name = :model,
+        initialization_eqs = [s ~ 0, v ~ 0],
+    ))
+
+    prob = ODEProblem{true, FullSpecialize}(sys, [], (0.0, 200.0))
+    sol = solve(prob, Tsit5(); saveat = 0.1)
+    @test SciMLBase.successful_retcode(sol)
+
+    setter = setp_oop(prob, [sys.m, sys.d])
+
+    function loss(x)
+        p = setter(prob, x)
+        newprob = remake(prob; p)
+        newsol = solve(newprob, Tsit5(); saveat = 0.1)
+        sum(abs2, newsol[sys.s])
+    end
+
+    ForwardDiff.gradient(loss, [3.0, 20.0])
+end