Finished pushing changes & tests to remaining system types.

Author: lamorton

src/systems/jumps/jumpsystem.jl

@@ -158,7 +158,7 @@ end
 
 function generate_rate_function(js::JumpSystem, rate)
     consts = collect_constants(rate)
-    if !isempty(consts) # The SymbolicUtils._build_function method of this case doesn't support preprocessing
+    if !isempty(consts) # The SymbolicUtils._build_function method of this case doesn't support postprocess_fbody
         csubs = Dict(c => getdefault(c) for c in consts)
         rate = substitute(rate, csubs)
     end
@@ -170,7 +170,7 @@ end
 
 function generate_affect_function(js::JumpSystem, affect, outputidxs)
     consts = collect_constants(affect)
-    if !isempty(consts) # The SymbolicUtils._build_function method of this case doesn't support preprocessing
+    if !isempty(consts) # The SymbolicUtils._build_function method of this case doesn't support postprocess_fbody
         csubs = Dict(c => getdefault(c) for c in consts)
         affect = substitute(affect, csubs)
     end

src/systems/nonlinear/nonlinearsystem.jl

@@ -148,7 +148,8 @@ end
 function generate_jacobian(sys::NonlinearSystem, vs = states(sys), ps = parameters(sys);
                            sparse = false, simplify = false, kwargs...)
     jac = calculate_jacobian(sys, sparse = sparse, simplify = simplify)
-    return build_function(jac, vs, ps; kwargs...)
+    pre = get_preprocess_constants(jac)
+    return build_function(jac, vs, ps; postprocess_fbody = pre, kwargs...)
 end
 
 function calculate_hessian(sys::NonlinearSystem; sparse = false, simplify = false)
@@ -165,7 +166,8 @@ end
 function generate_hessian(sys::NonlinearSystem, vs = states(sys), ps = parameters(sys);
                           sparse = false, simplify = false, kwargs...)
     hess = calculate_hessian(sys, sparse = sparse, simplify = simplify)
-    return build_function(hess, vs, ps; kwargs...)
+    pre = get_preprocess_constants(hess)
+    return build_function(hess, vs, ps; postprocess_fbody = pre, kwargs...)
 end
 
 function generate_function(sys::NonlinearSystem, dvs = states(sys), ps = parameters(sys);

src/systems/optimization/optimizationsystem.jl

@@ -100,7 +100,8 @@ end
 function generate_gradient(sys::OptimizationSystem, vs = states(sys), ps = parameters(sys);
                            kwargs...)
     grad = calculate_gradient(sys)
-    return build_function(grad, vs, ps;
+    pre = get_preprocess_constants(grad)
+    return build_function(grad, vs, ps; postprocess_fbody = pre,
                           conv = AbstractSysToExpr(sys), kwargs...)
 end
 
@@ -115,13 +116,20 @@ function generate_hessian(sys::OptimizationSystem, vs = states(sys), ps = parame
     else
         hess = calculate_hessian(sys)
     end
-    return build_function(hess, vs, ps;
+    pre = get_preprocess_constants(hess)
+    return build_function(hess, vs, ps; postprocess_fbody = pre,
                           conv = AbstractSysToExpr(sys), kwargs...)
 end
 
 function generate_function(sys::OptimizationSystem, vs = states(sys), ps = parameters(sys);
                            kwargs...)
-    return build_function(equations(sys), vs, ps;
+    eqs = equations(sys)
+    consts = collect_constants(eqs)
+    if !isempty(consts) # The SymbolicUtils._build_function method of this case doesn't support postprocess_fbody
+        csubs = Dict(c => getdefault(c) for c in consts)
+        eqs = substitute(eqs, csubs)
+    end
+    return build_function(eqs, vs, ps;
                           conv = AbstractSysToExpr(sys), kwargs...)
 end
 
@@ -186,7 +194,6 @@ symbolically calculating numerical enhancements.
 function DiffEqBase.OptimizationProblem(sys::OptimizationSystem, args...; kwargs...)
     DiffEqBase.OptimizationProblem{true}(sys::OptimizationSystem, args...; kwargs...)
 end
-
 function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
                                              parammap = DiffEqBase.NullParameters();
                                              lb = nothing, ub = nothing,
@@ -211,8 +218,13 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
 
     f = generate_function(sys, checkbounds = checkbounds, linenumbers = linenumbers,
                           expression = Val{false})
-
-    obj_expr = toexpr(equations(sys))
+    eqs = equations(sys)
+    cs = collect_constants(eqs)
+    if !isempty(cs)
+        cmap = map(x -> x => getdefault(x), cs)
+        eqs = substitute(eqs, cmap)
+    end
+    obj_expr = toexpr(eqs)
     pairs_arr = p isa SciMLBase.NullParameters ?
                 [Symbol(_s) => Expr(:ref, :x, i) for (i, _s) in enumerate(dvs)] :
                 [
@@ -254,8 +266,13 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
                                  expression = Val{false})[2]
         cons_j = generate_jacobian(cons_sys; expression = Val{false}, sparse = sparse)[2]
         cons_h = generate_hessian(cons_sys; expression = Val{false}, sparse = sparse)[2]
-
-        cons_expr = toexpr(equations(cons_sys))
+        eqs = equations(cons_sys)
+        cs = collect_constants(eqs)
+        if !isempty(cs)
+            cmap = map(x -> x => getdefault(x), cs)
+            eqs = map(x -> x.lhs ~ substitute(x.rhs, cmap), eqs)
+        end
+        cons_expr = toexpr(eqs)
         rep_pars_vals!.(cons_expr, Ref(pairs_arr))
 
         if sparse
@@ -265,7 +282,6 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
             cons_jac_prototype = nothing
             cons_hess_prototype = nothing
         end
-
         _f = DiffEqBase.OptimizationFunction{iip}(f,
                                                   sys = sys,
                                                   syms = nameof.(states(sys)),
@@ -360,16 +376,23 @@ function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0,
     lb = varmap_to_vars(lb, dvs; check = false, tofloat = false, use_union)
     ub = varmap_to_vars(ub, dvs; check = false, tofloat = false, use_union)
 
-    obj_expr = toexpr(equations(sys))
+    eqs = equations(sys)
+    cs = collect_constants(eqs)
+    if !isempty(cs)
+        cmap = map(x -> x => getdefault(x), cs)
+        eqs = substitute(eqs, cmap)
+    end
+    obj_expr = toexpr(eqs)
     pairs_arr = p isa SciMLBase.NullParameters ?
                 [Symbol(_s) => Expr(:ref, :x, i) for (i, _s) in enumerate(dvs)] :
                 [
         [Symbol(_s) => Expr(:ref, :x, i) for (i, _s) in enumerate(dvs)]...,
         [Symbol(_p) => p[i] for (i, _p) in enumerate(ps)]...,
     ]
     rep_pars_vals!(obj_expr, pairs_arr)
-
+    @show sys.constraints
     if length(sys.constraints) > 0
+        
         @named cons_sys = NonlinearSystem(sys.constraints, dvs, ps)
         cons = generate_function(cons_sys, checkbounds = checkbounds,
                                  linenumbers = linenumbers,
@@ -378,7 +401,13 @@ function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0,
 
         cons_h = generate_hessian(cons_sys; expression = Val{false}, sparse = sparse)[2]
 
-        cons_expr = toexpr(equations(cons_sys))
+        eqs = equations(cons_sys)
+        cs = collect_constants(eqs)
+        if !isempty(cs)
+            cmap = map(x -> x => getdefault(x), cs)
+            eqs = map(x -> x.lhs ~ substitute(x.rhs, cmap), eqs)
+        end
+        cons_expr = toexpr(eqs)
         rep_pars_vals!.(cons_expr, Ref(pairs_arr))
 
         if sparse

src/utils.jl

@@ -513,14 +513,25 @@ function collect_constants(eqs::Vector{Equation}) #For get_substitutions_and_sol
     return constants
 end
 
-function collect_constants(eqs::AbstractArray{T}) where {T} # For generate_tgrad / generate_jacobian / generate_difference_cb
+function collect_constants(eqs::AbstractArray{T}) where {T <: Union{Num, Symbolic}} # For generate_tgrad / generate_jacobian / generate_difference_cb
     constants = T[]
     for eq in eqs
         collect_constants!(constants, unwrap(eq))
     end
     return constants
 end
 
+function collect_constants(eqs::Vector{Matrix{Num}}) # For nonlinear hessian
+    constants = Num[]
+    for m in eqs
+        for n in m
+            collect_constants!(constants, unwrap(n))
+        end
+    end
+    return constants
+end
+
+
 collect_constants(x::Num) = collect_constants(unwrap(x))
 function collect_constants(expr::Symbolic{T}) where {T} # For jump system affect / rate 
     constants = Symbolic[]

test/optimizationsystem.jl

@@ -2,16 +2,17 @@ using ModelingToolkit, SparseArrays, Test, Optimization, OptimizationOptimJL,
       OptimizationMOI, Ipopt, AmplNLWriter, Ipopt_jll
 
 @variables x y
+@constants h=1 
 @parameters a b
-loss = (a - x)^2 + b * (y - x^2)^2
+loss = (a - x)^2 + b * (y * h - x^2)^2
 sys1 = OptimizationSystem(loss, [x, y], [a, b], name = :sys1)
 
-cons2 = [x^2 + y^2 ~ 0, y * sin(x) - x ~ 0]
+cons2 = [x^2 + y^2 ~ 0, y * sin(x) - x * h ~ 0]
 sys2 = OptimizationSystem(loss, [x, y], [a, b], name = :sys2, constraints = cons2)
 
 @variables z
 @parameters β
-loss2 = sys1.x - sys2.y + z * β
+loss2 = sys1.x - sys2.y + z * β * h
 combinedsys = OptimizationSystem(loss2, [z], [β], systems = [sys1, sys2],
                                  name = :combinedsys)
 
@@ -64,7 +65,7 @@ sol = solve(prob, AmplNLWriter.Optimizer(Ipopt_jll.amplexe))
 @test sol.minimum < 1.0
 
 #equality constraint, lcons == ucons
-cons2 = [0.0 ~ x^2 + y^2]
+cons2 = [0.0 ~ h * x^2 + y^2]
 out = zeros(1)
 sys2 = OptimizationSystem(loss, [x, y], [a, b], name = :sys2, constraints = cons2)
 prob = OptimizationProblem(sys2, [x => 0.0, y => 0.0], [a => 1.0, b => 1.0], lcons = [1.0],
@@ -98,7 +99,7 @@ end
 
 # observed variable handling
 @variables OBS
-@named sys2 = OptimizationSystem(loss, [x, y], [a, b]; observed = [OBS ~ x + y])
+@named sys2 = OptimizationSystem(loss, [x, y], [a, b]; observed = [OBS ~ x * h + y])
 OBS2 = OBS
 @test isequal(OBS2, @nonamespace sys2.OBS)
 @unpack OBS = sys2

test/pde.jl

@@ -2,11 +2,12 @@ using ModelingToolkit, DiffEqBase, LinearAlgebra
 
 # Define some variables
 @parameters t x
+@constants h=1
 @variables u(..)
 Dt = Differential(t)
 Dxx = Differential(x)^2
-eq = Dt(u(t, x)) ~ Dxx(u(t, x))
-bcs = [u(0, x) ~ -x * (x - 1) * sin(x),
+eq = Dt(u(t, x)) ~ h * Dxx(u(t, x))
+bcs = [u(0, x) ~ -h * x * (x - 1) * sin(x),
     u(t, 0) ~ 0, u(t, 1) ~ 0]
 
 domains = [t ∈ (0.0, 1.0),