diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl
index 4d01bd31c9..a62580024e 100644
--- a/src/systems/optimization/optimizationsystem.jl
+++ b/src/systems/optimization/optimizationsystem.jl
@@ -114,8 +114,8 @@ function OptimizationSystem(op, unknowns, ps;
             unknowns′[i] = irrvar
         end
     end
-    op′ = substitute(op′, irreducible_subs)
-    constraints = substitute.(constraints, (irreducible_subs,))
+    op′ = fast_substitute(op′, irreducible_subs)
+    constraints = fast_substitute.(constraints, (irreducible_subs,))
 
     if !(isempty(default_u0) && isempty(default_p))
         Base.depwarn(
@@ -127,7 +127,7 @@ function OptimizationSystem(op, unknowns, ps;
         throw(ArgumentError("System names must be unique."))
     end
     defaults = todict(defaults)
-    defaults = Dict(substitute(value(k), irreducible_subs) => substitute(
+    defaults = Dict(fast_substitute(value(k), irreducible_subs) => fast_substitute(
                         value(v), irreducible_subs)
     for (k, v) in pairs(defaults) if value(v) !== nothing)
 
diff --git a/test/optimizationsystem.jl b/test/optimizationsystem.jl
index c20613441a..2ec9516721 100644
--- a/test/optimizationsystem.jl
+++ b/test/optimizationsystem.jl
@@ -1,5 +1,6 @@
 using ModelingToolkit, SparseArrays, Test, Optimization, OptimizationOptimJL,
-      OptimizationMOI, Ipopt, AmplNLWriter, Ipopt_jll, SymbolicIndexingInterface
+      OptimizationMOI, Ipopt, AmplNLWriter, Ipopt_jll, SymbolicIndexingInterface,
+      LinearAlgebra
 using ModelingToolkit: get_metadata
 
 @testset "basic" begin
@@ -388,3 +389,23 @@ end
     @test all(y -> any(x -> isequal(x, y), unknowns(sys2)), [x2, sys1.x1])
     @test all(y -> any(x -> isequal(x, y), parameters(sys2)), [p2, sys1.p1])
 end
+
+function myeigvals_1(A::AbstractMatrix)
+    eigvals(A)[1]
+end
+
+@register_symbolic myeigvals_1(A::AbstractMatrix)
+
+@testset "Issue#3473: Registered array function in objective, no irreducible variables" begin
+    p_free = @variables begin
+        p1, [bounds = (0, 1)]
+        p2, [bounds = (0, 1)]
+        p3, [bounds = (0, 1)]
+        p4, [bounds = (0, 1)]
+    end
+
+    m = diagm(p_free)
+
+    obj = myeigvals_1(m)
+    @test_nowarn OptimizationSystem(obj, p_free, []; name = :osys)
+end