correct error messages

Author: vyudu

src/systems/diffeqs/abstractodesystem.jl

@@ -697,6 +697,7 @@ end
 ```julia
 DiffEqBase.ODEProblem{iip}(sys::AbstractODESystem, u0map, tspan,
                            parammap = DiffEqBase.NullParameters();
+                           allow_cost = false,
                            version = nothing, tgrad = false,
                            jac = false,
                            checkbounds = false, sparse = false,
@@ -730,6 +731,7 @@ end
 function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map = [],
         tspan = get_tspan(sys),
         parammap = DiffEqBase.NullParameters();
+        allow_cost = false,
         callback = nothing,
         check_length = true,
         warn_initialize_determined = true,
@@ -745,9 +747,10 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map =
               Consider a BVProblem instead.")
     end
 
-    if !isempty(get_costs(sys))
-        error("An ODESystem with costs cannot be solved using a regular ODEProblem.
-              Solvers for optimal control problems are forthcoming.")
+    if !isempty(get_costs(sys)) && !allow_cost
+        error("ODEProblem will not optimize solutions of ODESystems that have associated cost functions. 
+              Solvers for optimal control problems are forthcoming. In order to bypass this error (e.g. 
+              to check the cost of a regular solution), pass `allow_cost` = true into the constructor.")
     end
 
     f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap;
@@ -801,21 +804,19 @@ If an ODESystem without `constraints` is specified, it will be treated as an ini
 
 ```julia
     @parameters g t_c = 0.5
-    @variables x(..) y(t) [state_priority = 10] λ(t)
+    @variables x(..) y(t) λ(t)
     eqs = [D(D(x(t))) ~ λ * x(t)
            D(D(y)) ~ λ * y - g
            x(t)^2 + y^2 ~ 1]
     cstr = [x(0.5) ~ 1]
-    @named cstrs = ConstraintsSystem(cstr, t)
-    @mtkbuild pend = ODESystem(eqs, t)
+    @mtkbuild pend = ODESystem(eqs, t; constraints = cstrs)
 
     tspan = (0.0, 1.5)
     u0map = [x(t) => 0.6, y => 0.8]
     parammap = [g => 1]
     guesses = [λ => 1]
-    constraints = [x(0.5) ~ 1]
 
-    bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; constraints, guesses, check_length = false)
+    bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, check_length = false)
 ```
 
 If the `ODESystem` has algebraic equations, like `x(t)^2 + y(t)^2`, the resulting 
@@ -844,6 +845,7 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = []
         tspan = get_tspan(sys),
         parammap = DiffEqBase.NullParameters();
         guesses = Dict(),
+        allow_cost = false,
         version = nothing, tgrad = false,
         callback = nothing,
         check_length = true,
@@ -857,9 +859,10 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = []
     end
     !isnothing(callback) && error("BVP solvers do not support callbacks.")
 
-    if !isempty(get_costs(sys))
-        error("An ODESystem with costs cannot be solved using a regular DAEProblem.
-              Solvers for optimal control problems are forthcoming.")
+    if !isempty(get_costs(sys)) && !allow_cost
+        error("BVProblem will not optimize solutions of ODESystems that have associated cost functions.
+              Solvers for optimal control problems are forthcoming. In order to bypass this error (e.g.
+              to check the cost of a regular solution), pass `allow_cost` = true into the constructor.")
     end
 
     has_alg_eqs(sys) &&
@@ -958,16 +961,19 @@ end
 
 function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan,
         parammap = DiffEqBase.NullParameters();
+        allow_cost = false,
         warn_initialize_determined = true,
         check_length = true, eval_expression = false, eval_module = @__MODULE__, kwargs...) where {iip}
     if !iscomplete(sys)
         error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `DAEProblem`")
     end
 
-    if !isempty(get_costs(sys))
-        error("An ODESystem with costs cannot be solved using a regular DAEProblem.
-              Solvers for optimal control problems are forthcoming.")
+    if !isempty(get_costs(sys)) && !allow_cost
+        error("DAEProblem will not optimize solutions of ODESystems that have associated cost functions. 
+              Solvers for optimal control problems are forthcoming. In order to bypass this error (e.g. 
+              to check the cost of a regular solution), pass `allow_cost` = true into the constructor.")
     end
+
     f, du0, u0, p = process_SciMLProblem(DAEFunction{iip}, sys, u0map, parammap;
         implicit_dae = true, du0map = du0map, check_length,
         t = tspan !== nothing ? tspan[1] : tspan,

test/optimal_control.jl

@@ -290,7 +290,8 @@ end
     @mtkbuild lksys = ODESystem(eqs, t; costs, consolidate)
     @test_throws ErrorException @mtkbuild lksys2 = ODESystem(eqs, t; costs)
 
-    prob = ODEProblem(lksys, u0map, tspan, parammap)
+    @test_throws ErrorException ODEProblem(lksys, u0map, tspan, parammap)
+    prob = ODEProblem(lksys, u0map, tspan, parammap; allow_cost = true)
     sol = solve(prob, Tsit5())
     costfn = ModelingToolkit.generate_cost_function(lksys)
     _t = tspan[2]
@@ -303,7 +304,7 @@ end
     @mtkbuild lksys = ODESystem(eqs, t; costs, consolidate)
     @test t_c ∈ Set(parameters(lksys))
     push!(parammap, t_c => 0.56)
-    prob = ODEProblem(lksys, u0map, tspan, parammap)
+    prob = ODEProblem(lksys, u0map, tspan, parammap; allow_cost = true)
     sol = solve(prob, Tsit5())
     costfn = ModelingToolkit.generate_cost_function(lksys)
     @test costfn(sol, prob.p, _t) ≈ log(sol(0.56)[2] + sol(0.)[1]) - sol(0.4)[1]^2