diff --git a/src/inputoutput.jl b/src/inputoutput.jl
index e07d2ed976..6ef36cddcc 100644
--- a/src/inputoutput.jl
+++ b/src/inputoutput.jl
@@ -218,6 +218,7 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
     inputs = map(x -> time_varying_as_func(value(x), sys), inputs)
 
     eqs = [eq for eq in full_equations(sys)]
+    eqs = map(subs_constants, eqs)
     if disturbance_inputs !== nothing
         # Set all disturbance *inputs* to zero (we just want to keep the disturbance state)
         subs = Dict(disturbance_inputs .=> 0)
diff --git a/test/input_output_handling.jl b/test/input_output_handling.jl
index 778e02db0a..19078bc98c 100644
--- a/test/input_output_handling.jl
+++ b/test/input_output_handling.jl
@@ -391,3 +391,14 @@ matrices, ssys = linearize(augmented_sys,
         io_sys, [x + u * t]; inputs = [u])
     @test obsfn([1.0], [2.0], nothing, 3.0) == [7.0]
 end
+
+# https://github.com/SciML/ModelingToolkit.jl/issues/2896
+@testset "Constants substitution" begin
+    @constants c = 2.0
+    @variables x(t)
+    eqs = [D(x) ~ c * x]
+    @named sys = ODESystem(eqs, t, [x], [])
+
+    f, dvs, ps = ModelingToolkit.generate_control_function(sys, simplify = true)
+    @test f[1]([0.5], nothing, nothing, 0.0) == [1.0]
+end