diff --git a/.gitignore b/.gitignore index 6d305ad348..9193389ab8 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ Manifest.toml .vscode/* docs/src/assets/Project.toml docs/src/assets/Manifest.toml +.CondaPkg diff --git a/Project.toml b/Project.toml index b661dc402d..9538bb0894 100644 --- a/Project.toml +++ b/Project.toml @@ -64,6 +64,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [weakdeps] BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665" +CasADi = "c49709b8-5c63-11e9-2fb2-69db5844192f" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" FMI = "14a09403-18e3-468f-ad8a-74f8dda2d9ac" @@ -72,6 +73,7 @@ LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" [extensions] MTKBifurcationKitExt = "BifurcationKit" +MTKCasADiDynamicOptExt = "CasADi" MTKChainRulesCoreExt = "ChainRulesCore" MTKDeepDiffsExt = "DeepDiffs" MTKFMIExt = "FMI" @@ -86,6 +88,7 @@ BifurcationKit = "0.4" BlockArrays = "1.1" BoundaryValueDiffEqAscher = "1.1.0" BoundaryValueDiffEqMIRK = "1.4.0" +CasADi = "1.0.6" ChainRulesCore = "1" Combinatorics = "1" CommonSolve = "0.2.4" diff --git a/ext/MTKCasADiDynamicOptExt.jl b/ext/MTKCasADiDynamicOptExt.jl new file mode 100644 index 0000000000..9d16a3ea5c --- /dev/null +++ b/ext/MTKCasADiDynamicOptExt.jl @@ -0,0 +1,376 @@ +module MTKCasADiDynamicOptExt +using ModelingToolkit +using CasADi +using DiffEqBase +using UnPack +using NaNMath +const MTK = ModelingToolkit + +# NaNMath +for ff in [acos, log1p, acosh, log2, asin, tan, atanh, cos, log, sin, log10, sqrt] + f = nameof(ff) + # These need to be defined so that JuMP can trace through functions built by Symbolics + @eval NaNMath.$f(x::CasadiSymbolicObject) = Base.$f(x) +end + +# Default linear interpolation for MX objects, likely to change down the line when we support interpolation with the collocation polynomial. +struct MXLinearInterpolation + u::MX + t::Vector{Float64} + dt::Float64 +end + +struct CasADiModel + opti::Opti + U::MXLinearInterpolation + V::MXLinearInterpolation + tₛ::MX +end + +struct CasADiDynamicOptProblem{uType, tType, isinplace, P, F, K} <: + AbstractDynamicOptProblem{uType, tType, isinplace} + f::F + u0::uType + tspan::tType + p::P + model::CasADiModel + kwargs::K + + function CasADiDynamicOptProblem(f, u0, tspan, p, model, kwargs...) + new{typeof(u0), typeof(tspan), SciMLBase.isinplace(f, 5), + typeof(p), typeof(f), typeof(kwargs)}(f, u0, tspan, p, model, kwargs) + end +end + +function (M::MXLinearInterpolation)(τ) + nt = (τ - M.t[1]) / M.dt + i = 1 + floor(Int, nt) + Δ = nt - i + 1 + + (i > length(M.t) || i < 1) && error("Cannot extrapolate past the tspan.") + if i < length(M.t) + M.u[:, i] + Δ*(M.u[:, i + 1] - M.u[:, i]) + else + M.u[:, i] + end +end + +""" + CasADiDynamicOptProblem(sys::ODESystem, u0, tspan, p; dt, steps) + +Convert an ODESystem representing an optimal control system into a CasADi model +for solving using optimization. Must provide either `dt`, the timestep between collocation +points (which, along with the timespan, determines the number of points), or directly +provide the number of points as `steps`. + +The optimization variables: +- a vector-of-vectors U representing the unknowns as an interpolation array +- a vector-of-vectors V representing the controls as an interpolation array + +The constraints are: +- The set of user constraints passed to the ODESystem via `constraints` +- The solver constraints that encode the time-stepping used by the solver +""" +function MTK.CasADiDynamicOptProblem(sys::ODESystem, u0map, tspan, pmap; + dt = nothing, + steps = nothing, + guesses = Dict(), kwargs...) + MTK.warn_overdetermined(sys, u0map) + _u0map = has_alg_eqs(sys) ? u0map : merge(Dict(u0map), Dict(guesses)) + f, u0, p = MTK.process_SciMLProblem(ODEInputFunction, sys, _u0map, pmap; + t = tspan !== nothing ? tspan[1] : tspan, output_type = MX, kwargs...) + + pmap = Dict{Any, Any}(pmap) + steps, is_free_t = MTK.process_tspan(tspan, dt, steps) + model = init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t) + + CasADiDynamicOptProblem(f, u0, tspan, p, model, kwargs...) +end + +function init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t = false) + ctrls = MTK.unbound_inputs(sys) + states = unknowns(sys) + opti = CasADi.Opti() + + if is_free_t + (ts_sym, te_sym) = tspan + MTK.symbolic_type(ts_sym) !== MTK.NotSymbolic() && + error("Free initial time problems are not currently supported in CasADiDynamicOptProblem.") + tₛ = variable!(opti) + set_initial!(opti, tₛ, pmap[te_sym]) + subject_to!(opti, tₛ >= ts_sym) + hasbounds(te_sym) && begin + lo, hi = getbounds(te_sym) + subject_to!(opti, tₛ >= lo) + subject_to!(opti, tₛ >= hi) + end + pmap[te_sym] = tₛ + tsteps = LinRange(0, 1, steps) + else + tₛ = MX(1) + tsteps = LinRange(tspan[1], tspan[2], steps) + end + + U = CasADi.variable!(opti, length(states), steps) + V = CasADi.variable!(opti, length(ctrls), steps) + set_initial!(opti, U, DM(repeat(u0, 1, steps))) + c0 = MTK.value.([pmap[c] for c in ctrls]) + !isempty(c0) && set_initial!(opti, V, DM(repeat(c0, 1, steps))) + + U_interp = MXLinearInterpolation(U, tsteps, tsteps[2]-tsteps[1]) + V_interp = MXLinearInterpolation(V, tsteps, tsteps[2]-tsteps[1]) + for (i, ct) in enumerate(ctrls) + pmap[ct] = V[i, :] + end + + model = CasADiModel(opti, U_interp, V_interp, tₛ) + + set_casadi_bounds!(model, sys, pmap) + add_cost_function!(model, sys, tspan, pmap; is_free_t) + add_user_constraints!(model, sys, tspan, pmap; is_free_t) + + stidxmap = Dict([v => i for (i, v) in enumerate(states)]) + u0map = Dict([MTK.default_toterm(MTK.value(k)) => v for (k, v) in u0map]) + u0_idxs = has_alg_eqs(sys) ? collect(1:length(states)) : + [stidxmap[MTK.default_toterm(k)] for (k, v) in u0map] + add_initial_constraints!(model, u0, u0_idxs) + + model +end + +function set_casadi_bounds!(model, sys, pmap) + @unpack opti, U, V = model + for (i, u) in enumerate(unknowns(sys)) + if MTK.hasbounds(u) + lo, hi = MTK.getbounds(u) + subject_to!(opti, Symbolics.fixpoint_sub(lo, pmap) <= U.u[i, :]) + subject_to!(opti, U.u[i, :] <= Symbolics.fixpoint_sub(hi, pmap)) + end + end + for (i, v) in enumerate(MTK.unbound_inputs(sys)) + if MTK.hasbounds(v) + lo, hi = MTK.getbounds(v) + subject_to!(opti, Symbolics.fixpoint_sub(lo, pmap) <= V.u[i, :]) + subject_to!(opti, V.u[i, :] <= Symbolics.fixpoint_sub(hi, pmap)) + end + end +end + +function add_initial_constraints!(model::CasADiModel, u0, u0_idxs) + @unpack opti, U = model + for i in u0_idxs + subject_to!(opti, U.u[i, 1] == u0[i]) + end +end + +function add_user_constraints!(model::CasADiModel, sys, tspan, pmap; is_free_t) + @unpack opti, U, V, tₛ = model + + iv = MTK.get_iv(sys) + conssys = MTK.get_constraintsystem(sys) + jconstraints = isnothing(conssys) ? nothing : MTK.get_constraints(conssys) + (isnothing(jconstraints) || isempty(jconstraints)) && return nothing + + stidxmap = Dict([v => i for (i, v) in enumerate(unknowns(sys))]) + ctidxmap = Dict([v => i for (i, v) in enumerate(MTK.unbound_inputs(sys))]) + cons_unknowns = map(MTK.default_toterm, unknowns(conssys)) + + auxmap = Dict([u => MTK.default_toterm(MTK.value(u)) for u in unknowns(conssys)]) + jconstraints = substitute_casadi_vars(model, sys, pmap, jconstraints; is_free_t, auxmap) + # Manually substitute fixed-t variables + for (i, cons) in enumerate(jconstraints) + consvars = MTK.vars(cons) + for st in consvars + MTK.iscall(st) || continue + x = MTK.operation(st) + t = only(MTK.arguments(st)) + MTK.symbolic_type(t) === MTK.NotSymbolic() || continue + if haskey(stidxmap, x(iv)) + idx = stidxmap[x(iv)] + cv = U + else + idx = ctidxmap[x(iv)] + cv = V + end + cons = Symbolics.substitute(cons, Dict(x(t) => cv(t)[idx])) + end + + if cons isa Equation + subject_to!(opti, cons.lhs - cons.rhs==0) + elseif cons.relational_op === Symbolics.geq + subject_to!(opti, cons.lhs - cons.rhs≥0) + else + subject_to!(opti, cons.lhs - cons.rhs≤0) + end + end +end + +function add_cost_function!(model::CasADiModel, sys, tspan, pmap; is_free_t) + @unpack opti, U, V, tₛ = model + jcosts = copy(MTK.get_costs(sys)) + consolidate = MTK.get_consolidate(sys) + if isnothing(jcosts) || isempty(jcosts) + minimize!(opti, MX(0)) + return + end + + iv = MTK.get_iv(sys) + stidxmap = Dict([v => i for (i, v) in enumerate(unknowns(sys))]) + ctidxmap = Dict([v => i for (i, v) in enumerate(MTK.unbound_inputs(sys))]) + + jcosts = substitute_casadi_vars(model, sys, pmap, jcosts; is_free_t) + # Substitute fixed-time variables. + for i in 1:length(jcosts) + costvars = MTK.vars(jcosts[i]) + for st in costvars + MTK.iscall(st) || continue + x = operation(st) + t = only(arguments(st)) + MTK.symbolic_type(t) === MTK.NotSymbolic() || continue + if haskey(stidxmap, x(iv)) + idx = stidxmap[x(iv)] + cv = U + else + idx = ctidxmap[x(iv)] + cv = V + end + jcosts[i] = Symbolics.substitute(jcosts[i], Dict(x(t) => cv(t)[idx])) + end + end + + dt = U.t[2] - U.t[1] + intmap = Dict() + for int in MTK.collect_applied_operators(jcosts, Symbolics.Integral) + op = MTK.operation(int) + arg = only(arguments(MTK.value(int))) + lo, hi = (op.domain.domain.left, op.domain.domain.right) + !isequal((lo, hi), tspan) && error("Non-whole interval bounds for integrals are not currently supported for CasADiDynamicOptProblem.") + # Approximate integral as sum. + intmap[int] = dt * tₛ * sum(arg) + end + jcosts = map(c -> Symbolics.substitute(c, intmap), jcosts) + jcosts = MTK.value.(jcosts) + minimize!(opti, MX(MTK.value(consolidate(jcosts)))) +end + +function substitute_casadi_vars(model::CasADiModel, sys, pmap, exprs; auxmap::Dict = Dict(), is_free_t) + @unpack opti, U, V, tₛ = model + iv = MTK.get_iv(sys) + sts = unknowns(sys) + cts = MTK.unbound_inputs(sys) + + x_ops = [MTK.operation(MTK.unwrap(st)) for st in sts] + c_ops = [MTK.operation(MTK.unwrap(ct)) for ct in cts] + + exprs = map(c -> Symbolics.fixpoint_sub(c, auxmap), exprs) + exprs = map(c -> Symbolics.fixpoint_sub(c, Dict(pmap)), exprs) + # tf means different things in different contexts; a [tf] in a cost function + # should be tₛ, while a x(tf) should translate to x[1] + if is_free_t + free_t_map = Dict([[x(tₛ) => U.u[i, end] for (i, x) in enumerate(x_ops)]; + [c(tₛ) => V.u[i, end] for (i, c) in enumerate(c_ops)]]) + exprs = map(c -> Symbolics.fixpoint_sub(c, free_t_map), exprs) + end + + # for variables like x(t) + whole_interval_map = Dict([[v => U.u[i, :] for (i, v) in enumerate(sts)]; + [v => V.u[i, :] for (i, v) in enumerate(cts)]]) + exprs = map(c -> Symbolics.fixpoint_sub(c, whole_interval_map), exprs) + exprs +end + +function add_solve_constraints(prob, tableau) + @unpack A, α, c = tableau + @unpack model, f, p = prob + @unpack opti, U, V, tₛ = model + solver_opti = copy(opti) + + tsteps = U.t + dt = tsteps[2] - tsteps[1] + + nᵤ = size(U.u, 1) + nᵥ = size(V.u, 1) + + if MTK.is_explicit(tableau) + K = MX[] + for k in 1:length(tsteps)-1 + τ = tsteps[k] + for (i, h) in enumerate(c) + ΔU = sum([A[i, j] * K[j] for j in 1:(i - 1)], init = MX(zeros(nᵤ))) + Uₙ = U.u[:, k] + ΔU*dt + Vₙ = V.u[:, k] + Kₙ = tₛ * f(Uₙ, Vₙ, p, τ + h * dt) # scale the time + push!(K, Kₙ) + end + ΔU = dt * sum([α[i] * K[i] for i in 1:length(α)]) + subject_to!(solver_opti, U.u[:, k] + ΔU == U.u[:, k+1]) + empty!(K) + end + else + for k in 1:length(tsteps)-1 + τ = tsteps[k] + Kᵢ = variable!(solver_opti, nᵤ, length(α)) + ΔUs = A * Kᵢ' # the stepsize at each stage of the implicit method + for (i, h) in enumerate(c) + ΔU = ΔUs[i,:]' + Uₙ = U.u[:,k] + ΔU*dt + Vₙ = V.u[:,k] + subject_to!(solver_opti, Kᵢ[:,i] == tₛ * f(Uₙ, Vₙ, p, τ + h*dt)) + end + ΔU_tot = dt*(Kᵢ*α) + subject_to!(solver_opti, U.u[:, k] + ΔU_tot == U.u[:,k+1]) + end + end + solver_opti +end + +""" + solve(prob::CasADiDynamicOptProblem, casadi_solver, ode_solver; plugin_options, solver_options, silent) + +`plugin_options` and `solver_options` get propagated to the Opti object in CasADi. + +NOTE: the solver should be passed in as a string to CasADi. "ipopt" +""" +function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, Symbol} = "ipopt", tableau_getter = MTK.constructDefault; plugin_options::Dict = Dict(), solver_options::Dict = Dict(), silent = false) + @unpack model, u0, p, tspan, f = prob + tableau = tableau_getter() + @unpack opti, U, V, tₛ = model + + opti = add_solve_constraints(prob, tableau) + silent && (solver_options["print_level"] = 0) + solver!(opti, "$solver", plugin_options, solver_options) + + failed = false + value_getter = nothing + sol = nothing + try + sol = CasADi.solve!(opti) + value_getter = x -> CasADi.value(sol, x) + catch ErrorException + value_getter = x -> CasADi.debug_value(opti, x) + failed = true + end + + ts = value_getter(tₛ) * U.t + U_vals = value_getter(U.u) + size(U_vals, 2) == 1 && (U_vals = U_vals') + U_vals = [[U_vals[i, j] for i in 1:size(U_vals, 1)] for j in 1:length(ts)] + ode_sol = DiffEqBase.build_solution(prob, tableau_getter, ts, U_vals) + + input_sol = nothing + if prod(size(V.u)) != 0 + V_vals = value_getter(V.u) + size(V_vals, 2) == 1 && (V_vals = V_vals') + V_vals = [[V_vals[i, j] for i in 1:size(V_vals, 1)] for j in 1:length(ts)] + input_sol = DiffEqBase.build_solution(prob, tableau_getter, ts, V_vals) + end + + if failed + ode_sol = SciMLBase.solution_new_retcode(ode_sol, SciMLBase.ReturnCode.ConvergenceFailure) + !isnothing(input_sol) && (input_sol = SciMLBase.solution_new_retcode( + input_sol, SciMLBase.ReturnCode.ConvergenceFailure)) + end + + DynamicOptSolution(model, ode_sol, input_sol) +end +end diff --git a/ext/MTKInfiniteOptExt.jl b/ext/MTKInfiniteOptExt.jl index ba3a772582..e36b9c2a3e 100644 --- a/ext/MTKInfiniteOptExt.jl +++ b/ext/MTKInfiniteOptExt.jl @@ -257,8 +257,6 @@ function substitute_jump_vars(model, sys, pmap, exprs; auxmap = Dict(), is_free_ exprs end -is_explicit(tableau) = tableau isa DiffEqBase.ExplicitRKTableau - function add_infopt_solve_constraints!(model::InfiniteModel, sys, pmap; is_free_t = false) # Differential equations U = model[:U] @@ -295,7 +293,7 @@ function add_jump_solve_constraints!(prob, tableau; is_free_t = false) V = model[:V] nᵤ = length(U) nᵥ = length(V) - if is_explicit(tableau) + if MTK.is_explicit(tableau) K = Any[] for τ in tsteps for (i, h) in enumerate(c) @@ -327,23 +325,6 @@ function add_jump_solve_constraints!(prob, tableau; is_free_t = false) end end -""" -Default ODE Tableau: RadauIIA5 -""" -function constructDefault(T::Type = Float64) - sq6 = sqrt(6) - A = [11 // 45-7sq6 / 360 37 // 225-169sq6 / 1800 -2 // 225+sq6 / 75 - 37 // 225+169sq6 / 1800 11 // 45+7sq6 / 360 -2 // 225-sq6 / 75 - 4 // 9-sq6 / 36 4 // 9+sq6 / 36 1//9] - c = [2 // 5 - sq6 / 10; 2 / 5 + sq6 / 10; 1] - α = [4 // 9 - sq6 / 36; 4 // 9 + sq6 / 36; 1 // 9] - A = map(T, A) - α = map(T, α) - c = map(T, c) - - DiffEqBase.ImplicitRKTableau(A, c, α, 5) -end - """ Solve JuMPDynamicOptProblem. Arguments: - prob: a JumpDynamicOptProblem @@ -354,7 +335,7 @@ Solve JuMPDynamicOptProblem. Arguments: Returns a DynamicOptSolution, which contains both the model and the ODE solution. """ function DiffEqBase.solve( - prob::JuMPDynamicOptProblem, jump_solver, tableau_getter = constructDefault; silent = false) + prob::JuMPDynamicOptProblem, jump_solver, tableau_getter = MTK.constructDefault; silent = false) model = prob.model tableau = tableau_getter() silent && set_silent(model) diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index d2d20f902e..2f589c80e6 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -349,7 +349,7 @@ export AnalysisPoint, get_sensitivity_function, get_comp_sensitivity_function, function FMIComponent end include("systems/optimal_control_interface.jl") -export AbstractDynamicOptProblem, JuMPDynamicOptProblem, InfiniteOptDynamicOptProblem +export AbstractDynamicOptProblem, JuMPDynamicOptProblem, InfiniteOptDynamicOptProblem, CasADiDynamicOptProblem export DynamicOptSolution end # module diff --git a/src/systems/optimal_control_interface.jl b/src/systems/optimal_control_interface.jl index eb573da810..d979e7b64e 100644 --- a/src/systems/optimal_control_interface.jl +++ b/src/systems/optimal_control_interface.jl @@ -18,6 +18,7 @@ end function JuMPDynamicOptProblem end function InfiniteOptDynamicOptProblem end +function CasADiDynamicOptProblem end function warn_overdetermined(sys, u0map) constraintsys = get_constraintsystem(sys) @@ -27,6 +28,25 @@ function warn_overdetermined(sys, u0map) end end +""" +Default ODE Tableau: RadauIIA5 +""" +function constructDefault(T::Type = Float64) + sq6 = sqrt(6) + A = [11 // 45-7sq6 / 360 37 // 225-169sq6 / 1800 -2 // 225+sq6 / 75 + 37 // 225+169sq6 / 1800 11 // 45+7sq6 / 360 -2 // 225-sq6 / 75 + 4 // 9-sq6 / 36 4 // 9+sq6 / 36 1//9] + c = [2 // 5 - sq6 / 10; 2 / 5 + sq6 / 10; 1] + α = [4 // 9 - sq6 / 36; 4 // 9 + sq6 / 36; 1 // 9] + A = map(T, A) + α = map(T, α) + c = map(T, c) + + DiffEqBase.ImplicitRKTableau(A, c, α, 5) +end + +is_explicit(tableau) = tableau isa DiffEqBase.ExplicitRKTableau + """ Generate the control function f(x, u, p, t) from the ODESystem. Input variables are automatically inferred but can be manually specified. diff --git a/test/extensions/Project.toml b/test/extensions/Project.toml index 0e018b4a22..9f43e6f4a4 100644 --- a/test/extensions/Project.toml +++ b/test/extensions/Project.toml @@ -1,5 +1,6 @@ [deps] BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665" +CasADi = "c49709b8-5c63-11e9-2fb2-69db5844192f" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a" DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" @@ -25,3 +26,6 @@ SimpleDiffEq = "05bca326-078c-5bf0-a5bf-ce7c7982d7fd" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" + +[compat] +CasADi = "1.0.6" diff --git a/test/extensions/jump_control.jl b/test/extensions/dynamic_optimization.jl similarity index 83% rename from test/extensions/jump_control.jl rename to test/extensions/dynamic_optimization.jl index a93146be10..4b3b4edeb6 100644 --- a/test/extensions/jump_control.jl +++ b/test/extensions/dynamic_optimization.jl @@ -5,6 +5,9 @@ using SimpleDiffEq using OrdinaryDiffEqSDIRK, OrdinaryDiffEqVerner, OrdinaryDiffEqTsit5, OrdinaryDiffEqFIRK using Ipopt using DataInterpolations +using CasADi + +import DiffEqBase: solve const M = ModelingToolkit @testset "ODE Solution, no cost" begin @@ -28,17 +31,23 @@ const M = ModelingToolkit jsol = solve(jprob, Ipopt.Optimizer, constructRK4, silent = true) oprob = ODEProblem(sys, u0map, tspan, parammap) osol = solve(oprob, SimpleRK4(), dt = 0.01) + cprob = CasADiDynamicOptProblem(sys, u0map, tspan, parammap, dt = 0.01) + csol = solve(cprob, "ipopt", constructRK4) @test jsol.sol.u ≈ osol.u + @test csol.sol.u ≈ osol.u # Implicit method. jsol2 = solve(jprob, Ipopt.Optimizer, constructImplicitEuler, silent = true) # 63.031 ms, 26.49 MiB osol2 = solve(oprob, ImplicitEuler(), dt = 0.01, adaptive = false) # 129.375 μs, 61.91 KiB + jsol2 = solve(jprob, Ipopt.Optimizer, constructImplicitEuler, silent = true) # 63.031 ms, 26.49 MiB @test ≈(jsol2.sol.u, osol2.u, rtol = 0.001) iprob = InfiniteOptDynamicOptProblem(sys, u0map, tspan, parammap, dt = 0.01) isol = solve(iprob, Ipopt.Optimizer, derivative_method = InfiniteOpt.FiniteDifference(InfiniteOpt.Backward()), silent = true) # 11.540 ms, 4.00 MiB - @test ≈(jsol2.sol.u, osol2.u, rtol = 0.001) + @test ≈(isol.sol.u, osol2.u, rtol = 0.001) + csol2 = solve(cprob, "ipopt", constructImplicitEuler, silent = true) + @test ≈(csol2.sol.u, osol2.u, rtol = 0.001) # With a constraint u0map = Pair[] @@ -52,6 +61,11 @@ const M = ModelingToolkit @test jsol.sol(0.6)[1] ≈ 3.5 @test jsol.sol(0.3)[1] ≈ 7.0 + cprob = CasADiDynamicOptProblem(lksys, u0map, tspan, parammap; guesses = guess, dt = 0.01) + csol = solve(cprob, "ipopt", constructTsitouras5, silent = true) + @test csol.sol(0.6)[1] ≈ 3.5 + @test csol.sol(0.3)[1] ≈ 7.0 + iprob = InfiniteOptDynamicOptProblem( lksys, u0map, tspan, parammap; guesses = guess, dt = 0.01) isol = solve(iprob, Ipopt.Optimizer, @@ -66,12 +80,16 @@ const M = ModelingToolkit iprob = InfiniteOptDynamicOptProblem( lksys, u0map, tspan, parammap; guesses = guess, dt = 0.01) isol = solve(iprob, Ipopt.Optimizer, - derivative_method = InfiniteOpt.OrthogonalCollocation(3), silent = true) # 48.564 ms, 9.58 MiB + derivative_method = InfiniteOpt.OrthogonalCollocation(3), silent = true) @test all(u -> u > [1, 1], isol.sol.u) jprob = JuMPDynamicOptProblem(lksys, u0map, tspan, parammap; guesses = guess, dt = 0.01) jsol = solve(jprob, Ipopt.Optimizer, constructRadauIA3, silent = true) # 12.190 s, 9.68 GiB @test all(u -> u > [1, 1], jsol.sol.u) + + cprob = CasADiDynamicOptProblem(lksys, u0map, tspan, parammap; guesses = guess, dt = 0.01) + csol = solve(cprob, "ipopt", constructRadauIA3, silent = true) + @test all(u -> u > [1, 1], csol.sol.u) end function is_bangbang(input_sol, lbounds, ubounds, rtol = 1e-4) @@ -109,6 +127,14 @@ end @test is_bangbang(jsol.input_sol, [-1.0], [1.0]) # Test reached final position. @test ≈(jsol.sol.u[end][2], 0.25, rtol = 1e-5) + + cprob = CasADiDynamicOptProblem(block, u0map, tspan, parammap; dt = 0.01) + csol = solve(cprob, "ipopt", constructVerner8, silent = true) + # Linear systems have bang-bang controls + @test is_bangbang(csol.input_sol, [-1.0], [1.0]) + # Test reached final position. + @test ≈(csol.sol.u[end][2], 0.25, rtol = 1e-5) + # Test dynamics @parameters (u_interp::ConstantInterpolation)(..) @mtkbuild block_ode = ODESystem([D(x(t)) ~ v(t), D(v(t)) ~ u_interp(t)], t) @@ -116,6 +142,7 @@ end oprob = ODEProblem(block_ode, u0map, tspan, [u_interp => spline]) osol = solve(oprob, Vern8(), dt = 0.01, adaptive = false) @test ≈(jsol.sol.u, osol.u, rtol = 0.05) + @test ≈(csol.sol.u, osol.u, rtol = 0.05) iprob = InfiniteOptDynamicOptProblem(block, u0map, tspan, parammap; dt = 0.01) isol = solve(iprob, Ipopt.Optimizer; silent = true) @@ -147,6 +174,9 @@ end iprob = InfiniteOptDynamicOptProblem(beesys, u0map, tspan, pmap, dt = 0.01) isol = solve(iprob, Ipopt.Optimizer; silent = true) @test is_bangbang(isol.input_sol, [0.0], [1.0]) + cprob = CasADiDynamicOptProblem(beesys, u0map, tspan, pmap; dt = 0.01) + csol = solve(cprob, "ipopt", constructTsitouras5, silent = true) + @test is_bangbang(csol.input_sol, [0.0], [1.0]) @parameters (α_interp::LinearInterpolation)(..) eqs = [D(w(t)) ~ -μ * w(t) + b * s * α_interp(t) * w(t), @@ -159,6 +189,7 @@ end Dict(α_interp => ctrl_to_spline(jsol.input_sol, LinearInterpolation)))) osol = solve(oprob, Tsit5(); dt = 0.01, adaptive = false) @test ≈(osol.u, jsol.sol.u, rtol = 0.01) + @test ≈(osol.u, csol.sol.u, rtol = 0.01) osol2 = solve(oprob, ImplicitEuler(); dt = 0.01, adaptive = false) @test ≈(osol2.u, isol.sol.u, rtol = 0.01) end @@ -189,6 +220,10 @@ end jprob = JuMPDynamicOptProblem(rocket, u0map, (ts, te), pmap; dt = 0.001, cse = false) jsol = solve(jprob, Ipopt.Optimizer, constructRadauIIA5, silent = true) @test jsol.sol.u[end][1] > 1.012 + + cprob = CasADiDynamicOptProblem(rocket, u0map, (ts, te), pmap; dt = 0.001, cse = false) + csol = solve(cprob, "ipopt"; silent = true) + @test csol.sol.u[end][1] > 1.012 iprob = InfiniteOptDynamicOptProblem( rocket, u0map, (ts, te), pmap; dt = 0.001) @@ -205,6 +240,7 @@ end oprob = ODEProblem(rocket_ode, u0map, (ts, te), merge(Dict(pmap), interpmap)) osol = solve(oprob, RadauIIA5(); adaptive = false, dt = 0.001) @test ≈(jsol.sol.u, osol.u, rtol = 0.02) + @test ≈(csol.sol.u, osol.u, rtol = 0.02) interpmap1 = Dict(T_interp => ctrl_to_spline(isol.input_sol, CubicSpline)) oprob1 = ODEProblem(rocket_ode, u0map, (ts, te), merge(Dict(pmap), interpmap1)) @@ -231,6 +267,10 @@ end jsol = solve(jprob, Ipopt.Optimizer, constructTsitouras5, silent = true) @test isapprox(jsol.sol.t[end], 10.0, rtol = 1e-3) + cprob = CasADiDynamicOptProblem(rocket, u0map, (0, tf), pmap; steps = 201) + csol = solve(cprob, "ipopt", constructTsitouras5, silent = true) + @test isapprox(csol.sol.t[end], 10.0, rtol = 1e-3) + iprob = InfiniteOptDynamicOptProblem(rocket, u0map, (0, tf), pmap; steps = 200) isol = solve(iprob, Ipopt.Optimizer, silent = true) @test isapprox(isol.sol.t[end], 10.0, rtol = 1e-3) @@ -252,6 +292,10 @@ end jsol = solve(jprob, Ipopt.Optimizer, constructVerner8, silent = true) @test isapprox(jsol.sol.t[end], 2.0, atol = 1e-5) + cprob = CasADiDynamicOptProblem(block, u0map, (0, tf), parammap; steps = 51) + csol = solve(cprob, "ipopt", constructVerner8, silent = true) + @test isapprox(csol.sol.t[end], 2.0, atol = 1e-5) + iprob = InfiniteOptDynamicOptProblem(block, u0map, tspan, parammap; steps = 51) isol = solve(iprob, Ipopt.Optimizer, silent = true) @test isapprox(isol.sol.t[end], 2.0, atol = 1e-5) @@ -290,6 +334,10 @@ end jsol = solve(jprob, Ipopt.Optimizer, constructRK4, silent = true) @test jsol.sol.u[end] ≈ [π, 0, 0, 0] + cprob = CasADiDynamicOptProblem(cartpole, u0map, tspan, pmap; dt = 0.04) + csol = solve(cprob, "ipopt", constructRK4, silent = true) + @test csol.sol.u[end] ≈ [π, 0, 0, 0] + iprob = InfiniteOptDynamicOptProblem(cartpole, u0map, tspan, pmap; dt = 0.04) isol = solve(iprob, Ipopt.Optimizer, silent = true) @test isol.sol.u[end] ≈ [π, 0, 0, 0] diff --git a/test/runtests.jl b/test/runtests.jl index 80339c6606..1b8b5e03db 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -136,7 +136,7 @@ end if GROUP == "All" || GROUP == "Extensions" activate_extensions_env() - @safetestset "JuMP Collocation Solvers" include("extensions/jump_control.jl") + @safetestset "Dynamic Optimization Collocation Solvers" include("extensions/dynamic_optimization.jl") @safetestset "HomotopyContinuation Extension Test" include("extensions/homotopy_continuation.jl") @safetestset "Auto Differentiation Test" include("extensions/ad.jl") @safetestset "LabelledArrays Test" include("labelledarrays.jl")