From 1093c1cfe3a63a73953fd6c594825c7e0ca77ca3 Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 8 Apr 2025 00:40:15 -0400 Subject: [PATCH 01/28] feat: solver tableau debugging --- Project.toml | 4 + ext/MTKJuMPControlExt.jl | 218 ++++++++++++++++++++++++ test/extensions/Project.toml | 1 + test/extensions/jump_control.jl | 291 +++----------------------------- 4 files changed, 249 insertions(+), 265 deletions(-) create mode 100644 ext/MTKJuMPControlExt.jl diff --git a/Project.toml b/Project.toml index b661dc402d..b29b2b9566 100644 --- a/Project.toml +++ b/Project.toml @@ -31,6 +31,7 @@ FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e" FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5" Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" @@ -66,8 +67,10 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" +DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" FMI = "14a09403-18e3-468f-ad8a-74f8dda2d9ac" InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" [extensions] @@ -97,6 +100,7 @@ DeepDiffs = "1" DelayDiffEq = "5.50" DiffEqBase = "6.170.1" DiffEqCallbacks = "2.16, 3, 4" +DiffEqDevTools = "2.48.0" DiffEqNoiseProcess = "5" DiffRules = "0.1, 1.0" DifferentiationInterface = "0.6.47" diff --git a/ext/MTKJuMPControlExt.jl b/ext/MTKJuMPControlExt.jl new file mode 100644 index 0000000000..1e788b05ec --- /dev/null +++ b/ext/MTKJuMPControlExt.jl @@ -0,0 +1,218 @@ +module MTKJuMPControlExt +using ModelingToolkit +using JuMP, InfiniteOpt +using DiffEqDevTools, DiffEqBase, SciMLBase +using LinearAlgebra +const MTK = ModelingToolkit + +struct JuMPControlProblem{uType, tType, P, F, K} + f::F + u0::uType + tspan::tType + p::P + model::InfiniteModel + kwargs::K + + function JuMPControlProblem(f, u0, tspan, p, model; kwargs...) + new{typeof(u0), typeof(tspan), typeof(p), typeof(f), typeof(kwargs)}(f, u0, tspan, p, model, kwargs) + end +end + +""" + JuMPControlProblem(sys::ODESystem, u0, tspan, p; dt) + +Convert an ODESystem representing an optimal control system into a JuMP model +for solving using optimization. Must provide `dt` for determining the length +of the interpolation arrays. + +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.JuMPControlProblem(sys::ODESystem, u0map, tspan, pmap; dt = error("dt must be provided for JuMPControlProblem."), kwargs...) + ts = tspan[1] + te = tspan[2] + steps = ts:dt:te + ctrls = controls(sys) + states = unknowns(sys) + constraintsys = MTK.get_constraintsystem(sys) + + if !isnothing(constraintsys) + (length(constraints(constraintsys)) + length(u0map) > length(sts)) && + @warn "The JuMPControlProblem is overdetermined. The total number of conditions (# constraints + # fixed initial values given by u0map) exceeds the total number of states. The solvers will default to doing a nonlinear least-squares optimization." + end + + f, u0, p = MTK.process_SciMLProblem(ODEFunction, sys, u0map, pmap; + t = tspan !== nothing ? tspan[1] : tspan, kwargs...) + + model = InfiniteModel() + @infinite_parameter(model, t in [ts, te], num_supports = length(steps), derivative_method = OrthogonalCollocation(2)) + @variable(model, U[1:length(states)], Infinite(t), start = ts) + @variable(model, V[1:length(ctrls)], Infinite(t), start = ts) + + add_jump_cost_function!(model, sys) + add_user_constraints!(model, sys) + + stidxmap = Dict([v => i for (i, v) in enumerate(states)]) + u0_idxs = has_alg_eqs(sys) ? collect(1:length(states)) : [stidxmap[k] for (k, v) in u0map] + add_initial_constraints!(model, u0, u0_idxs, tspan) + + JuMPControlProblem(f, u0, tspan, p, model, kwargs...) +end + +function add_jump_cost_function!(model, sys) + jcosts = MTK.get_costs(sys) + consolidate = MTK.get_consolidate(sys) + if isnothing(consolidate) + @objective(model, Min, 0) + return + end + iv = MTK.get_iv(sys) + + stidxmap = Dict([v => i for (i, v) in enumerate(unknowns(sys))]) + cidxmap = Dict([v => i for (i, v) in enumerate(controls(sys))]) + + for st in unknowns(sys) + x = operation(st) + t = only(arguments(st)) + idx = stidxmap[x(iv)] + subval = isequal(t, iv) ? model[:U][idx] : model[:U][idx](t) + jcosts = Symbolics.substitute(jcosts, Dict(x(t) => subval)) + end + + for ct in controls(sys) + p = operation(ct) + t = only(arguments(ct)) + idx = cidxmap[p(iv)] + subval = isequal(t, iv) ? model[:V][idx] : model[:V][idx](t) + jcosts = Symbolics.substitute(jcosts, Dict(x(t) => subval)) + end + + @objective(model, Min, consolidate(jcosts)) +end + +function add_user_constraints!(model, sys) + jconstraints = if !(csys = MTK.get_constraintsystem(sys) isa Nothing) + MTK.get_constraints(csys) + else + nothing + end + isnothing(jconstraints) && return nothing + + iv = MTK.get_iv(sys) + stidxmap = Dict([v => i for (i, v) in enumerate(unknowns(sys))]) + cidxmap = Dict([v => i for (i, v) in enumerate(controls(sys))]) + + for st in unknowns(sys) + x = operation(st) + t = only(arguments(st)) + idx = stidxmap[x(iv)] + subval = isequal(t, iv) ? model[:U][idx] : model[:U][idx](t) + jconstraints = Symbolics.substitute(jconstraints, Dict(x(t) => subval)) + end + + for ct in controls(sys) + p = operation(ct) + t = only(arguments(ct)) + idx = cidxmap[p(iv)] + subval = isequal(t, iv) ? model[:V][idx] : model[:V][idx](t) + jconstraints = Symbolics.substitute(jconstraints, Dict(p(t) => subval)) + end + + for (i, cons) in enumerate(jconstraints) + if cons isa Equation + @constraint(model, user[i], cons.lhs - cons.rhs == 0) + elseif cons.relational_op === Symbolics.geq + @constraint(model, user[i], cons.lhs - cons.rhs ≥ 0) + else + @constraint(model, user[i], cons.lhs - cons.rhs ≤ 0) + end + end +end + +function add_initial_constraints!(model, u0, u0_idxs, tspan) + ts = tspan[1] + @constraint(model, init_u0_idx[i in u0_idxs], model[:U][i](ts) == u0[i]) +end + +is_explicit(tableau) = tableau isa DiffEqDevTools.ExplicitRKTableau + +function add_solve_constraints!(prob, tableau) + A = tableau.A + α = tableau.α + c = tableau.c + model = prob.model + f = prob.f + p = prob.p + tsteps = supports(model[:t]) + pop!(tsteps) + dt = tsteps[2] - tsteps[1] + + U = model[:U] + nᵤ = length(U) + if is_explicit(tableau) + K = Any[] + for τ in tsteps + for (i, h) in enumerate(c) + ΔU = sum([A[i, j] * K[j] for j in 1:i-1], init = zeros(nᵤ)) + Uₙ = [U[i](τ) + ΔU[i]*dt for i in 1:nᵤ] + Kₙ = f(Uₙ, p, τ + h*dt) + push!(K, Kₙ) + end + ΔU = sum([α[i] * K[i] for i in 1:length(α)]) + @constraint(model, [n = 1:nᵤ], U[n](τ) + ΔU[n] == U[n](τ + dt)) + empty!(K) + end + else + @variable(model, K[1:length(a), 1:nᵤ], Infinite(t), start = tsteps[1]) + for τ in tsteps + ΔUs = [A * K(τ)] + for (i, h) in enumerate(c) + ΔU = ΔUs[i] + Uₙ = [U[j](τ) + ΔU[j](τ)*dt for j in 1:nᵤ] + @constraint(model, K[i](τ) == f(Uₙ, p, τ + h*dt)) + end + ΔU = sum([α[i] * K[i] for i in 1:length(α)]) + @constraint(model, U(τ) + dot(α, K(τ)) == U(τ + dt)) + end + end +end + +""" +""" +struct JuMPControlSolution + model::InfiniteModel + sol::ODESolution +end + +""" +Solve JuMPControlProblem. Arguments: +- prob: a JumpControlProblem +- jump_solver: a LP solver such as HiGHS +- ode_solver: Takes in a symbol representing the solver. Acceptable solvers may be found at https://docs.sciml.ai/DiffEqDevDocs/stable/internals/tableaus/. Note that the symbol may be different than the typical name of the solver, e.g. :Tsitouras5 rather than Tsit5. +""" +function DiffEqBase.solve(prob::JuMPControlProblem, jump_solver, ode_solver::Symbol) + model = prob.model + tableau_getter = Symbol(:construct, ode_solver) + @eval tableau = $tableau_getter() + ts = supports(model[:t]) + add_solve_constraints!(prob, tableau) + + set_optimizer(model, jump_solver) + optimize!(model) + + if is_solved_and_feasible(model) + sol = DiffEqBase.build_solution(prob, ode_solver, ts, value.(U)) + JuMPControlSolution(model, sol) + else + sol = DiffEqBase.build_solution(prob, ode_solver, ts, value.(U)) + sol = SciMLBase.solution_new_retcode(sol, SciMLBase.ReturnCode.ConvergenceFailure) + JuMPControlSolution(model, sol) + end +end + +end diff --git a/test/extensions/Project.toml b/test/extensions/Project.toml index 0e018b4a22..bd51ea5203 100644 --- a/test/extensions/Project.toml +++ b/test/extensions/Project.toml @@ -6,6 +6,7 @@ DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" diff --git a/test/extensions/jump_control.jl b/test/extensions/jump_control.jl index a93146be10..c0a071e31d 100644 --- a/test/extensions/jump_control.jl +++ b/test/extensions/jump_control.jl @@ -1,296 +1,57 @@ using ModelingToolkit -import JuMP, InfiniteOpt +using JuMP, InfiniteOpt using DiffEqDevTools, DiffEqBase using SimpleDiffEq -using OrdinaryDiffEqSDIRK, OrdinaryDiffEqVerner, OrdinaryDiffEqTsit5, OrdinaryDiffEqFIRK -using Ipopt -using DataInterpolations +using HiGHS const M = ModelingToolkit @testset "ODE Solution, no cost" begin # Test solving without anything attached. @parameters α=1.5 β=1.0 γ=3.0 δ=1.0 - @variables x(..) y(..) - t = M.t_nounits - D = M.D_nounits + M.@variables x(..) y(..) + t = M.t_nounits; D = M.D_nounits eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y(t), D(y(t)) ~ -γ * y(t) + δ * x(t) * y(t)] - @mtkbuild sys = ODESystem(eqs, t) tspan = (0.0, 1.0) u0map = [x(t) => 4.0, y(t) => 2.0] - parammap = [α => 1.5, β => 1.0, γ => 3.0, δ => 1.0] + parammap = [α => 7.5, β => 4, γ => 8.0, δ => 5.0] + @mtkbuild sys = ODESystem(eqs, t) # Test explicit method. - jprob = JuMPDynamicOptProblem(sys, u0map, tspan, parammap, dt = 0.01) - @test JuMP.num_constraints(jprob.model) == 2 # initials - jsol = solve(jprob, Ipopt.Optimizer, constructRK4, silent = true) + jprob = JuMPControlProblem(sys, u0map, tspan, parammap, dt = 0.01) + @test num_constraints(jprob.model) == 2 # initials + jsol = solve(jprob, Ipopt.Optimizer, :Tsitouras5) oprob = ODEProblem(sys, u0map, tspan, parammap) - osol = solve(oprob, SimpleRK4(), dt = 0.01) + osol = solve(oprob, SimpleTsit5(), adaptive = false) @test jsol.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 - @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) + jsol2 = solve(prob, Ipopt.Optimizer, :RK4) + osol2 = solve(oprob, SimpleRK4(), adaptive = false) + @test jsol2.sol.u ≈ osol2.u # With a constraint - u0map = Pair[] + @parameters α=1.5 β=1.0 γ=3.0 δ=1.0 + @variables x(..) y(..) + + eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y(t), + D(y(t)) ~ -γ * y(t) + δ * x(t) * y(t)] + + u0map = [] + tspan = (0.0, 1.0) guess = [x(t) => 4.0, y(t) => 2.0] constr = [x(0.6) ~ 3.5, x(0.3) ~ 7.0] @mtkbuild lksys = ODESystem(eqs, t; constraints = constr) - jprob = JuMPDynamicOptProblem(lksys, u0map, tspan, parammap; guesses = guess, dt = 0.01) - @test JuMP.num_constraints(jprob.model) == 2 - jsol = solve(jprob, Ipopt.Optimizer, constructTsitouras5, silent = true) # 12.190 s, 9.68 GiB - @test jsol.sol(0.6)[1] ≈ 3.5 - @test jsol.sol(0.3)[1] ≈ 7.0 - - 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 - sol = isol.sol + jprob = JuMPControlProblem(sys, u0map, tspan, parammap; guesses, dt = 0.01) + @test num_constraints(jprob.model) == 2 == num_variables(jprob.model) == 2 + jsol = solve(prob, HiGHS.Optimizer, :Tsitouras5) + sol = jsol.sol @test sol(0.6)[1] ≈ 3.5 @test sol(0.3)[1] ≈ 7.0 - - # Test whole-interval constraints - constr = [x(t) ≳ 1, y(t) ≳ 1] - @mtkbuild lksys = ODESystem(eqs, t; constraints = constr) - 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 - @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) -end - -function is_bangbang(input_sol, lbounds, ubounds, rtol = 1e-4) - for v in 1:(length(input_sol.u[1]) - 1) - all(i -> ≈(i[v], bounds[v]; rtol) || ≈(i[v], bounds[u]; rtol), input_sol.u) || - return false - end - true end -function ctrl_to_spline(inputsol, splineType) - us = reduce(vcat, inputsol.u) - ts = reduce(vcat, inputsol.t) - splineType(us, ts) -end - -@testset "Linear systems" begin - # Double integrator - t = M.t_nounits - D = M.D_nounits - @variables x(..) [bounds = (0.0, 0.25)] v(..) - @variables u(..) [bounds = (-1.0, 1.0), input = true] - constr = [v(1.0) ~ 0.0] - cost = [-x(1.0)] # Maximize the final distance. - @named block = ODESystem( - [D(x(t)) ~ v(t), D(v(t)) ~ u(t)], t; costs = cost, constraints = constr) - block, input_idxs = structural_simplify(block, ([u(t)], [])) - - u0map = [x(t) => 0.0, v(t) => 0.0] - tspan = (0.0, 1.0) - parammap = [u(t) => 0.0] - jprob = JuMPDynamicOptProblem(block, u0map, tspan, parammap; dt = 0.01) - jsol = solve(jprob, Ipopt.Optimizer, constructVerner8, silent = true) - # Linear systems have bang-bang controls - @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) - # Test dynamics - @parameters (u_interp::ConstantInterpolation)(..) - @mtkbuild block_ode = ODESystem([D(x(t)) ~ v(t), D(v(t)) ~ u_interp(t)], t) - spline = ctrl_to_spline(jsol.input_sol, ConstantInterpolation) - 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) - - iprob = InfiniteOptDynamicOptProblem(block, u0map, tspan, parammap; dt = 0.01) - isol = solve(iprob, Ipopt.Optimizer; silent = true) - @test is_bangbang(isol.input_sol, [-1.0], [1.0]) - @test ≈(isol.sol.u[end][2], 0.25, rtol = 1e-5) - osol = solve(oprob, ImplicitEuler(); dt = 0.01, adaptive = false) - @test ≈(isol.sol.u, osol.u, rtol = 0.05) - - ################### - ### Bee example ### - ################### - # From Lawrence Evans' notes - @variables w(..) q(..) α(t) [input = true, bounds = (0, 1)] - @parameters b c μ s ν - - tspan = (0, 4) - eqs = [D(w(t)) ~ -μ * w(t) + b * s * α * w(t), - D(q(t)) ~ -ν * q(t) + c * (1 - α) * s * w(t)] - costs = [-q(tspan[2])] - - @named beesys = ODESystem(eqs, t; costs) - beesys, input_idxs = structural_simplify(beesys, ([α], [])) - u0map = [w(t) => 40, q(t) => 2] - pmap = [b => 1, c => 1, μ => 1, s => 1, ν => 1, α => 1] - - jprob = JuMPDynamicOptProblem(beesys, u0map, tspan, pmap, dt = 0.01) - jsol = solve(jprob, Ipopt.Optimizer, constructTsitouras5, silent = true) - @test is_bangbang(jsol.input_sol, [0.0], [1.0]) - 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]) - - @parameters (α_interp::LinearInterpolation)(..) - eqs = [D(w(t)) ~ -μ * w(t) + b * s * α_interp(t) * w(t), - D(q(t)) ~ -ν * q(t) + c * (1 - α_interp(t)) * s * w(t)] - @mtkbuild beesys_ode = ODESystem(eqs, t) - oprob = ODEProblem(beesys_ode, - u0map, - tspan, - merge(Dict(pmap), - 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) - osol2 = solve(oprob, ImplicitEuler(); dt = 0.01, adaptive = false) - @test ≈(osol2.u, isol.sol.u, rtol = 0.01) -end - -@testset "Rocket launch" begin - t = M.t_nounits - D = M.D_nounits - - @parameters h_c m₀ h₀ g₀ D_c c Tₘ m_c - @variables h(..) v(..) m(..) [bounds = (m_c, 1)] T(..) [input = true, bounds = (0, Tₘ)] - drag(h, v) = D_c * v^2 * exp(-h_c * (h - h₀) / h₀) - gravity(h) = g₀ * (h₀ / h) - - eqs = [D(h(t)) ~ v(t), - D(v(t)) ~ (T(t) - drag(h(t), v(t))) / m(t) - gravity(h(t)), - D(m(t)) ~ -T(t) / c] - - (ts, te) = (0.0, 0.2) - costs = [-h(te)] - cons = [T(te) ~ 0, m(te) ~ m_c] - @named rocket = ODESystem(eqs, t; costs, constraints = cons) - rocket, input_idxs = structural_simplify(rocket, ([T(t)], [])) - - u0map = [h(t) => h₀, m(t) => m₀, v(t) => 0] - pmap = [ - g₀ => 1, m₀ => 1.0, h_c => 500, c => 0.5 * √(g₀ * h₀), D_c => 0.5 * 620 * m₀ / g₀, - Tₘ => 3.5 * g₀ * m₀, T(t) => 0.0, h₀ => 1, m_c => 0.6] - 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 - - iprob = InfiniteOptDynamicOptProblem( - rocket, u0map, (ts, te), pmap; dt = 0.001) - isol = solve(iprob, Ipopt.Optimizer, silent = true) - @test isol.sol.u[end][1] > 1.012 - - # Test solution - @parameters (T_interp::CubicSpline)(..) - eqs = [D(h(t)) ~ v(t), - D(v(t)) ~ (T_interp(t) - drag(h(t), v(t))) / m(t) - gravity(h(t)), - D(m(t)) ~ -T_interp(t) / c] - @mtkbuild rocket_ode = ODESystem(eqs, t) - interpmap = Dict(T_interp => ctrl_to_spline(jsol.input_sol, CubicSpline)) - 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) - - interpmap1 = Dict(T_interp => ctrl_to_spline(isol.input_sol, CubicSpline)) - oprob1 = ODEProblem(rocket_ode, u0map, (ts, te), merge(Dict(pmap), interpmap1)) - osol1 = solve(oprob1, ImplicitEuler(); adaptive = false, dt = 0.001) - @test ≈(isol.sol.u, osol1.u, rtol = 0.01) -end - -@testset "Free final time problems" begin - t = M.t_nounits - D = M.D_nounits - - @variables x(..) u(..) [input = true, bounds = (0, 1)] - @parameters tf - eqs = [D(x(t)) ~ -2 + 0.5 * u(t)] - # Integral cost function - costs = [-Symbolics.Integral(t in (0, tf))(x(t) - u(t)), -x(tf)] - consolidate(u) = u[1] + u[2] - @named rocket = ODESystem(eqs, t; costs, consolidate) - rocket, input_idxs = structural_simplify(rocket, ([u(t)], [])) - - u0map = [x(t) => 17.5] - pmap = [u(t) => 0.0, tf => 8] - jprob = JuMPDynamicOptProblem(rocket, u0map, (0, tf), pmap; steps = 201) - jsol = solve(jprob, Ipopt.Optimizer, constructTsitouras5, silent = true) - @test isapprox(jsol.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) - - @variables x(..) v(..) - @variables u(..) [bounds = (-1.0, 1.0), input = true] - @parameters tf - constr = [v(tf) ~ 0, x(tf) ~ 0] - cost = [tf] # Minimize time - - @named block = ODESystem( - [D(x(t)) ~ v(t), D(v(t)) ~ u(t)], t; costs = cost, constraints = constr) - block, input_idxs = structural_simplify(block, ([u(t)], [])) - - u0map = [x(t) => 1.0, v(t) => 0.0] - tspan = (0.0, tf) - parammap = [u(t) => 0.0, tf => 1.0] - jprob = JuMPDynamicOptProblem(block, u0map, tspan, parammap; steps = 51) - jsol = solve(jprob, Ipopt.Optimizer, constructVerner8, silent = true) - @test isapprox(jsol.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) -end - -@testset "Cart-pole problem" begin - t = M.t_nounits - D = M.D_nounits - # gravity, length, moment of Inertia, drag coeff - @parameters g l mₚ mₖ - @variables x(..) θ(..) u(t) [input = true, bounds = (-10, 10)] - - s = sin(θ(t)) - c = cos(θ(t)) - H = [mₖ+mₚ mₚ*l*c - mₚ*l*c mₚ*l^2] - C = [0 -mₚ*D(θ(t))*l*s - 0 0] - qd = [D(x(t)), D(θ(t))] - G = [0, mₚ * g * l * s] - B = [1, 0] - - tf = 5 - rhss = -H \ Vector(C * qd + G - B * u) - eqs = [D(D(x(t))) ~ rhss[1], D(D(θ(t))) ~ rhss[2]] - cons = [θ(tf) ~ π, x(tf) ~ 0, D(θ(tf)) ~ 0, D(x(tf)) ~ 0] - costs = [Symbolics.Integral(t in (0, tf))(u^2)] - tspan = (0, tf) - - @named cartpole = ODESystem(eqs, t; costs, constraints = cons) - cartpole, input_idxs = structural_simplify(cartpole, ([u], [])) - - u0map = [D(x(t)) => 0.0, D(θ(t)) => 0.0, θ(t) => 0.0, x(t) => 0.0] - pmap = [mₖ => 1.0, mₚ => 0.2, l => 0.5, g => 9.81, u => 0] - jprob = JuMPDynamicOptProblem(cartpole, u0map, tspan, pmap; dt = 0.04) - jsol = solve(jprob, Ipopt.Optimizer, constructRK4, silent = true) - @test jsol.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] +@testset "Optimal control problem" begin end From d075449f86139c64d6988f75ee75f9a667c679ce Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 23 Apr 2025 19:09:18 -0400 Subject: [PATCH 02/28] test: more test fixes --- ext/MTKJuMPControlExt.jl | 153 ++++++++++++++++++++++++++----- src/systems/diffeqs/odesystem.jl | 1 + src/systems/systems.jl | 5 + test/runtests.jl | 11 +++ 4 files changed, 146 insertions(+), 24 deletions(-) diff --git a/ext/MTKJuMPControlExt.jl b/ext/MTKJuMPControlExt.jl index 1e788b05ec..631ec71ff0 100644 --- a/ext/MTKJuMPControlExt.jl +++ b/ext/MTKJuMPControlExt.jl @@ -18,6 +18,21 @@ struct JuMPControlProblem{uType, tType, P, F, K} end end +struct InfiniteOptControlProblem{uType, tType, isinplace, P, F, K} <: + AbstractOptimalControlProblem{uType, tType, isinplace} + f::F + u0::uType + tspan::tType + p::P + model::InfiniteModel + kwargs::K + + function InfiniteOptControlProblem(f, u0, tspan, p, model, kwargs...) + new{typeof(u0), typeof(tspan), SciMLBase.isinplace(f), + typeof(p), typeof(f), typeof(kwargs)}(f, u0, tspan, p, model, kwargs) + end +end + """ JuMPControlProblem(sys::ODESystem, u0, tspan, p; dt) @@ -33,20 +48,40 @@ 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.JuMPControlProblem(sys::ODESystem, u0map, tspan, pmap; dt = error("dt must be provided for JuMPControlProblem."), kwargs...) - ts = tspan[1] - te = tspan[2] - steps = ts:dt:te - ctrls = controls(sys) - states = unknowns(sys) - constraintsys = MTK.get_constraintsystem(sys) +function MTK.JuMPControlProblem(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, kwargs...) if !isnothing(constraintsys) (length(constraints(constraintsys)) + length(u0map) > length(sts)) && @warn "The JuMPControlProblem is overdetermined. The total number of conditions (# constraints + # fixed initial values given by u0map) exceeds the total number of states. The solvers will default to doing a nonlinear least-squares optimization." end - f, u0, p = MTK.process_SciMLProblem(ODEFunction, sys, u0map, pmap; + JuMPControlProblem(f, u0, tspan, p, model, kwargs...) +end + +""" + InfiniteOptControlProblem(sys::ODESystem, u0map, tspan, pmap; dt) + +Convert an ODESystem representing an optimal control system into a InfiniteOpt model +for solving using optimization. Must provide `dt` for determining the length +of the interpolation arrays. + +Related to `JuMPControlProblem`, but directly adds the differential equations +of the system as derivative constraints, rather than using a solver tableau. +""" +function MTK.InfiniteOptControlProblem(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, kwargs...) model = InfiniteModel() @@ -57,9 +92,22 @@ function MTK.JuMPControlProblem(sys::ODESystem, u0map, tspan, pmap; dt = error(" add_jump_cost_function!(model, sys) add_user_constraints!(model, sys) + @infinite_parameter(model, t in [tspan[1], tspan[2]], num_supports=steps) + @variable(model, U[i = 1:length(states)], Infinite(t), start=u0[i]) + c0 = MTK.value.([pmap[c] for c in ctrls]) + @variable(model, V[i = 1:length(ctrls)], Infinite(t), start=c0[i]) + + set_jump_bounds!(model, sys, pmap) + add_jump_cost_function!(model, sys, (tspan[1], tspan[2]), pmap; is_free_t) + add_user_constraints!(model, sys, pmap; is_free_t) + stidxmap = Dict([v => i for (i, v) in enumerate(states)]) - u0_idxs = has_alg_eqs(sys) ? collect(1:length(states)) : [stidxmap[k] for (k, v) in u0map] - add_initial_constraints!(model, u0, u0_idxs, tspan) + 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, tspan[1]) + return model +end JuMPControlProblem(f, u0, tspan, p, model, kwargs...) end @@ -115,14 +163,10 @@ function add_user_constraints!(model, sys) jconstraints = Symbolics.substitute(jconstraints, Dict(x(t) => subval)) end - for ct in controls(sys) - p = operation(ct) - t = only(arguments(ct)) - idx = cidxmap[p(iv)] - subval = isequal(t, iv) ? model[:V][idx] : model[:V][idx](t) - jconstraints = Symbolics.substitute(jconstraints, Dict(p(t) => subval)) - end - + auxmap = Dict([u => MTK.default_toterm(MTK.value(u)) for u in unknowns(conssys)]) + jconstraints = substitute_jump_vars(model, sys, pmap, jconstraints; auxmap) + + # Substitute to-term'd variables for (i, cons) in enumerate(jconstraints) if cons isa Equation @constraint(model, user[i], cons.lhs - cons.rhs == 0) @@ -139,6 +183,31 @@ function add_initial_constraints!(model, u0, u0_idxs, tspan) @constraint(model, init_u0_idx[i in u0_idxs], model[:U][i](ts) == u0[i]) end +function substitute_jump_vars(model, sys, pmap, exprs; auxmap = Dict()) + iv = MTK.get_iv(sys) + sts = unknowns(sys) + cts = MTK.unbound_inputs(sys) + U = model[:U] + V = model[:V] + exprs = map(c -> Symbolics.fixpoint_sub(c, auxmap), exprs) + + # for variables like x(t) + whole_interval_map = Dict([[v => U[i] for (i, v) in enumerate(sts)]; + [v => V[i] for (i, v) in enumerate(cts)]]) + exprs = map(c -> Symbolics.fixpoint_sub(c, whole_interval_map), exprs) + + # for variables like x(1.0) + x_ops = [MTK.operation(MTK.unwrap(st)) for st in sts] + c_ops = [MTK.operation(MTK.unwrap(ct)) for ct in cts] + fixed_t_map = Dict([[x_ops[i] => U[i] for i in 1:length(U)]; + [c_ops[i] => V[i] for i in 1:length(V)]]) + + exprs = map(c -> Symbolics.fixpoint_sub(c, fixed_t_map), exprs) + + exprs = map(c -> Symbolics.fixpoint_sub(c, Dict(pmap)), exprs) + exprs +end + is_explicit(tableau) = tableau isa DiffEqDevTools.ExplicitRKTableau function add_solve_constraints!(prob, tableau) @@ -148,7 +217,10 @@ function add_solve_constraints!(prob, tableau) model = prob.model f = prob.f p = prob.p - tsteps = supports(model[:t]) + + t = model[:t] + tsteps = supports(t) + tmax = tsteps[end] pop!(tsteps) dt = tsteps[2] - tsteps[1] @@ -167,17 +239,19 @@ function add_solve_constraints!(prob, tableau) @constraint(model, [n = 1:nᵤ], U[n](τ) + ΔU[n] == U[n](τ + dt)) empty!(K) end + @show num_variables(model) else @variable(model, K[1:length(a), 1:nᵤ], Infinite(t), start = tsteps[1]) for τ in tsteps ΔUs = [A * K(τ)] for (i, h) in enumerate(c) - ΔU = ΔUs[i] - Uₙ = [U[j](τ) + ΔU[j](τ)*dt for j in 1:nᵤ] - @constraint(model, K[i](τ) == f(Uₙ, p, τ + h*dt)) + ΔU = @view ΔUs[i, :] + Uₙ = U + ΔU * dt + @constraint(model, [j = 1:nᵤ], K[i, j](τ)==(tₛ * f(Uₙ, V, p, τ + h * dt)[j]), + DomainRestrictions(t => min(τ + h * dt, tmax)), base_name="solve_K($τ)") end - ΔU = sum([α[i] * K[i] for i in 1:length(α)]) - @constraint(model, U(τ) + dot(α, K(τ)) == U(τ + dt)) + @constraint(model, [n = 1:nᵤ], U[n](τ) + ΔU_tot[n]==U[n](min(τ + dt, tmax)), + DomainRestrictions(t => τ), base_name="solve_U($τ)") end end end @@ -202,6 +276,37 @@ function DiffEqBase.solve(prob::JuMPControlProblem, jump_solver, ode_solver::Sym ts = supports(model[:t]) add_solve_constraints!(prob, tableau) + # Unregister current solver constraints + for con in all_constraints(model) + if occursin("solve", JuMP.name(con)) + unregister(model, Symbol(JuMP.name(con))) + delete(model, con) + end + end + unregister(model, :K) + for var in all_variables(model) + if occursin("K", JuMP.name(var)) + unregister(model, Symbol(JuMP.name(var))) + delete(model, var) + end + end + add_jump_solve_constraints!(prob, tableau; is_free_t = haskey(model, :tf)) + _solve(prob, jump_solver, ode_solver) +end + +""" +`derivative_method` kwarg refers to the method used by InfiniteOpt to compute derivatives. The list of possible options can be found at https://infiniteopt.github.io/InfiniteOpt.jl/stable/guide/derivative/. Defaults to FiniteDifference(Backward()). +""" +function DiffEqBase.solve(prob::InfiniteOptControlProblem, jump_solver; + derivative_method = InfiniteOpt.FiniteDifference(Backward()), silent = false) + model = prob.model + silent && set_silent(model) + set_derivative_method(model[:t], derivative_method) + _solve(prob, jump_solver, derivative_method) +end + +function _solve(prob::AbstractOptimalControlProblem, jump_solver, solver) + model = prob.model set_optimizer(model, jump_solver) optimize!(model) diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index 4a13d7ccf1..c13670b5c7 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -767,6 +767,7 @@ function process_constraint_system( collect_vars!(constraintsts, constraintps, cons, iv) union!(constraintsts, collect_applied_operators(cons, Differential)) end + @show constraintsts # Validate the states. validate_vars_and_find_ps!(constraintsts, constraintps, sts, iv) diff --git a/src/systems/systems.jl b/src/systems/systems.jl index 52f93afb9b..9da7249300 100644 --- a/src/systems/systems.jl +++ b/src/systems/systems.jl @@ -163,6 +163,11 @@ function __structural_simplify( end end +function toterm_auxsystems(system::ODESystem) + constraints = system.constraintsystem.constraints + +end + """ $(TYPEDSIGNATURES) diff --git a/test/runtests.jl b/test/runtests.jl index 80339c6606..620394cc9e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -22,6 +22,12 @@ function activate_downstream_env() Pkg.instantiate() end +function activate_dynamic_optimization_env() + Pkg.activate("dynamic_optimization") + Pkg.develop(PackageSpec(path = dirname(@__DIR__))) + Pkg.instantiate() +end + @time begin if GROUP == "All" || GROUP == "InterfaceI" @testset "InterfaceI" begin @@ -143,4 +149,9 @@ end @safetestset "BifurcationKit Extension Test" include("extensions/bifurcationkit.jl") @safetestset "InfiniteOpt Extension Test" include("extensions/test_infiniteopt.jl") end + + if GROUP == "All" || GROUP == "Dynamic Optimization" + activate_dynamic_optimization_env() + @safetestset "JuMP Collocation Solvers" include("dynamic_optimization/jump_control") + end end From e983d4f9c2ac7a2eda5e1c75d44624921f0b31df Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 25 Apr 2025 14:56:55 -0400 Subject: [PATCH 03/28] more test fixes --- ext/MTKJuMPControlExt.jl | 11 ++++++----- src/systems/diffeqs/odesystem.jl | 1 - 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/ext/MTKJuMPControlExt.jl b/ext/MTKJuMPControlExt.jl index 631ec71ff0..de3aa96269 100644 --- a/ext/MTKJuMPControlExt.jl +++ b/ext/MTKJuMPControlExt.jl @@ -239,16 +239,17 @@ function add_solve_constraints!(prob, tableau) @constraint(model, [n = 1:nᵤ], U[n](τ) + ΔU[n] == U[n](τ + dt)) empty!(K) end - @show num_variables(model) else - @variable(model, K[1:length(a), 1:nᵤ], Infinite(t), start = tsteps[1]) + @variable(model, K[1:length(α), 1:nᵤ], Infinite(t)) + ΔUs = A * K + ΔU_tot = dt * (K' * α) for τ in tsteps ΔUs = [A * K(τ)] for (i, h) in enumerate(c) ΔU = @view ΔUs[i, :] - Uₙ = U + ΔU * dt - @constraint(model, [j = 1:nᵤ], K[i, j](τ)==(tₛ * f(Uₙ, V, p, τ + h * dt)[j]), - DomainRestrictions(t => min(τ + h * dt, tmax)), base_name="solve_K($τ)") + Uₙ = U + ΔU * h * dt + @constraint(model, [j = 1:nᵤ], K[i, j]==(tₛ * f(Uₙ, V, p, τ + h * dt)[j]), + DomainRestrictions(t => τ), base_name="solve_K$i($τ)") end @constraint(model, [n = 1:nᵤ], U[n](τ) + ΔU_tot[n]==U[n](min(τ + dt, tmax)), DomainRestrictions(t => τ), base_name="solve_U($τ)") diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index c13670b5c7..4a13d7ccf1 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -767,7 +767,6 @@ function process_constraint_system( collect_vars!(constraintsts, constraintps, cons, iv) union!(constraintsts, collect_applied_operators(cons, Differential)) end - @show constraintsts # Validate the states. validate_vars_and_find_ps!(constraintsts, constraintps, sts, iv) From d558fabb1383176545e95d89be3a6362ce7a5680 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 25 Apr 2025 23:06:54 -0400 Subject: [PATCH 04/28] test fixes --- src/structural_transformation/utils.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl index fa02f0b3cd..0d8a12a06d 100644 --- a/src/structural_transformation/utils.jl +++ b/src/structural_transformation/utils.jl @@ -96,6 +96,7 @@ function check_consistency(state::TransformationState, orig_inputs; nothrow = fa fullvars = get_fullvars(state) neqs = n_concrete_eqs(state) @unpack graph, var_to_diff = state.structure + @show equations(state.sys) highest_vars = computed_highest_diff_variables(complete!(state.structure)) n_highest_vars = 0 for (v, h) in enumerate(highest_vars) From 4ac62b4905662aa6e9c367708bec0881a14ebacb Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 28 Apr 2025 15:20:43 -0400 Subject: [PATCH 05/28] fix more tests --- src/structural_transformation/utils.jl | 1 - test/runtests.jl | 11 ----------- 2 files changed, 12 deletions(-) diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl index 0d8a12a06d..fa02f0b3cd 100644 --- a/src/structural_transformation/utils.jl +++ b/src/structural_transformation/utils.jl @@ -96,7 +96,6 @@ function check_consistency(state::TransformationState, orig_inputs; nothrow = fa fullvars = get_fullvars(state) neqs = n_concrete_eqs(state) @unpack graph, var_to_diff = state.structure - @show equations(state.sys) highest_vars = computed_highest_diff_variables(complete!(state.structure)) n_highest_vars = 0 for (v, h) in enumerate(highest_vars) diff --git a/test/runtests.jl b/test/runtests.jl index 620394cc9e..80339c6606 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -22,12 +22,6 @@ function activate_downstream_env() Pkg.instantiate() end -function activate_dynamic_optimization_env() - Pkg.activate("dynamic_optimization") - Pkg.develop(PackageSpec(path = dirname(@__DIR__))) - Pkg.instantiate() -end - @time begin if GROUP == "All" || GROUP == "InterfaceI" @testset "InterfaceI" begin @@ -149,9 +143,4 @@ end @safetestset "BifurcationKit Extension Test" include("extensions/bifurcationkit.jl") @safetestset "InfiniteOpt Extension Test" include("extensions/test_infiniteopt.jl") end - - if GROUP == "All" || GROUP == "Dynamic Optimization" - activate_dynamic_optimization_env() - @safetestset "JuMP Collocation Solvers" include("dynamic_optimization/jump_control") - end end From b1bc16094f49103786ef265724a325f6b621b8a0 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 2 May 2025 12:47:39 -0400 Subject: [PATCH 06/28] refactor: start systems test file --- ext/MTKCasADiDynamicOptExt.jl | 115 +++++++++++++++++++++++++ test/downstream/dynamic_opt_systems.jl | 19 ++++ 2 files changed, 134 insertions(+) create mode 100644 ext/MTKCasADiDynamicOptExt.jl create mode 100644 test/downstream/dynamic_opt_systems.jl diff --git a/ext/MTKCasADiDynamicOptExt.jl b/ext/MTKCasADiDynamicOptExt.jl new file mode 100644 index 0000000000..af2067896a --- /dev/null +++ b/ext/MTKCasADiDynamicOptExt.jl @@ -0,0 +1,115 @@ +module MTKCasADiDynamicOptExt +using ModelingToolkit +using CasADi +using DiffEqDevTools, DiffEqBase +using DataInterpolations +const MTK = MOdelingToolkit + +struct CasADiDynamicOptProblem{uType, tType, isinplace, P, F, K} <: + AbstractDynamicOptProblem{uType, tType, isinplace} + f::F + u0::uType + tspan::tType + p::P + model::Opti + 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 + +struct CasADiModel + opti::Opti + U::MX + V::MX +end + +struct TimedMX +end + +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, kwargs...) + + pmap = Dict{Any, Any}(pmap) + steps, is_free_t = MTK.process_tspan(tspan, dt, steps) + model = init_model() +end + +function init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t) + ctrls = MTK.unbound_inputs(sys) + states = unknowns(sys) + model = CasADi.Opti() + + U = CasADi.variable!(model, length(states), steps) + V = CasADi.variable!(model, length(ctrls), steps) +end + +function add_initial_constraints!() + +end + +function add_user_constraints!(model::CasADiModel, sys, pmap; is_free_t = false) + +end + +function add_cost_function!(model) + +end + +function add_solve_constraints!(prob, tableau; is_free_t) + A = tableau.A + α = tableau.α + c = tableau.c + model = prob.model + f = prob.f + p = prob.p + + opti = model.opti + t = model[:t] + tsteps = supports(t) + tmax = tsteps[end] + pop!(tsteps) + tₛ = is_free_t ? model[:tf] : 1 + dt = tsteps[2] - tsteps[1] + + U = model.U + V = model.V + nᵤ = length(U) + nᵥ = length(V) + + if is_explicit(tableau) + K = Any[] + for τ in tsteps + for (i, h) in enumerate(c) + ΔU = sum([A[i, j] * K[j] for j in 1:(i - 1)], init = zeros(nᵤ)) + Uₙ = [U[i](τ) + ΔU[i] * dt for i in 1:nᵤ] + Vₙ = [V[i](τ) for i in 1:nᵥ] + 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!(model, U[n](τ) + ΔU[n]==U[n](τ + dt)) + empty!(K) + end + else + end +end + +function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, Symbol}, ode_solver::Union{String, Symbol}; silent = false) + model = prob.model + opti = model.opti + + solver!(opti, solver) + sol = solve(opti) + DynamicOptSolution(model, sol, input_sol) +end + +end diff --git a/test/downstream/dynamic_opt_systems.jl b/test/downstream/dynamic_opt_systems.jl new file mode 100644 index 0000000000..30d8feb1cd --- /dev/null +++ b/test/downstream/dynamic_opt_systems.jl @@ -0,0 +1,19 @@ +function lotkavolterra() + +end + +function () + +end + +function rocket_fft() + +end + +function rocket() + +end + +function cartpole() + +end From 85e3ed666eab6943caade2465f9309d83afefb59 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 2 May 2025 20:59:04 -0400 Subject: [PATCH 07/28] feat: implementation for CasADiProblem --- Project.toml | 4 + ext/MTKCasADiDynamicOptExt.jl | 229 +++++++++++++++++++++---- test/downstream/dynamic_opt_systems.jl | 28 ++- 3 files changed, 219 insertions(+), 42 deletions(-) diff --git a/Project.toml b/Project.toml index b29b2b9566..2ea71a8ed3 100644 --- a/Project.toml +++ b/Project.toml @@ -65,7 +65,9 @@ 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" +DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" FMI = "14a09403-18e3-468f-ad8a-74f8dda2d9ac" @@ -75,6 +77,7 @@ LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" [extensions] MTKBifurcationKitExt = "BifurcationKit" +MTKCasADiDynamicOptExt = ["CasADi", "DataInterpolations", "DiffEqDevTools"] MTKChainRulesCoreExt = "ChainRulesCore" MTKDeepDiffsExt = "DeepDiffs" MTKFMIExt = "FMI" @@ -89,6 +92,7 @@ BifurcationKit = "0.4" BlockArrays = "1.1" BoundaryValueDiffEqAscher = "1.1.0" BoundaryValueDiffEqMIRK = "1.4.0" +CasADi = "1.0.0" ChainRulesCore = "1" Combinatorics = "1" CommonSolve = "0.2.4" diff --git a/ext/MTKCasADiDynamicOptExt.jl b/ext/MTKCasADiDynamicOptExt.jl index af2067896a..e16ce1b6b5 100644 --- a/ext/MTKCasADiDynamicOptExt.jl +++ b/ext/MTKCasADiDynamicOptExt.jl @@ -3,6 +3,7 @@ using ModelingToolkit using CasADi using DiffEqDevTools, DiffEqBase using DataInterpolations +using UnPack const MTK = MOdelingToolkit struct CasADiDynamicOptProblem{uType, tType, isinplace, P, F, K} <: @@ -22,16 +23,15 @@ end struct CasADiModel opti::Opti - U::MX - V::MX -end - -struct TimedMX + U::AbstractInterpolation + V::AbstractInterpolation + tₛ::Union{Number, MX} end function MTK.CasADiDynamicOptProblem(sys::ODESystem, u0map, tspan, pmap; dt = nothing, - steps = nothing, + steps = nothing, + interpolation_method::AbstractInterpolation = LinearInterpolation, guesses = Dict(), kwargs...) MTK.warn_overdetermined(sys, u0map) _u0map = has_alg_eqs(sys) ? u0map : merge(Dict(u0map), Dict(guesses)) @@ -43,73 +43,228 @@ function MTK.CasADiDynamicOptProblem(sys::ODESystem, u0map, tspan, pmap; model = init_model() end -function init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t) +function init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t = false, interp_type::AbstractInterpolation) ctrls = MTK.unbound_inputs(sys) states = unknowns(sys) - model = CasADi.Opti() + 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.") + tₛ = variable!(opti) + tsteps = LinRange(0, 1, steps) + else + tₛ = 1 + tsteps = LinRange(tspan[1], tspan[2], steps) + end - U = CasADi.variable!(model, length(states), steps) - V = CasADi.variable!(model, length(ctrls), steps) + U = CasADi.variable!(opti, length(states), steps) + V = CasADi.variable!(opti, length(ctrls), steps) + + U_interp = interp_type(Matrix(U), tsteps) + V_interp = interp_type(Matrix(V), tsteps) + + CasADiModel(opti, U_interp, V_interp, tₛ) end -function add_initial_constraints!() - +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, lo <= U[i, :] <= hi) + end + end + for (i, v) in enumerate(MTK.unbound_inputs(sys)) + if MTK.hasbounds(v) + lo, hi = MTK.getbounds(v) + subject_to!(opti, lo <= V[i, :] <= hi) + end + end +end + +function add_initial_constraints!(model::CasADiModel, u0, u0_idxs, ts) + @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, pmap; is_free_t = false) - + @unpack opti, U, V, tₛ = model + + iv = 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(sts)]) + pidxmap = Dict([v => i for (i, v) in enumerate(ps)]) + cons_unknowns = map(MTK.default_toterm, unknowns(conssys)) + for st in cons_unknowns + x = operation(st) + t = only(argments(st)) + idx = stidxmap[x(iv)] + + jconstraints = map(c -> Symbolics.substitute(c, Dict(x(t) => U(t)[idx])), jconstraints) + end + jconstraints = substitute_casadi_vars(model, sys, pmap, jconstraints) + + for (i, cons) in enumerate(jconstraints) + if cons isa Equation + subject_to!(opti, cons.lhs - cons.rhs==0) + elseif cons.relational_op === Symbolics.geq + subject_to!(model, cons.lhs - cons.rhs≥0) + else + subject_to!(model, cons.lhs - cons.rhs≤0) + end + end end -function add_cost_function!(model) +function add_cost_function!(model::CasADiModel, sys, tspan, pmap) + @unpack opti, U, V, tₛ = model + jcosts = MTK.get_costs(sys) + consolidate = MTK.get_consolidate(sys) + + if isnothing(jcosts) || isempty(jcosts) + minimize!(opti, 0) + return + end + stidxmap = Dict([v => i for (i, v) in enumerate(sts)]) + pidxmap = Dict([v => i for (i, v) in enumerate(ps)]) + for i in 1:length(jcosts) + vars = vars(jcosts[i]) + for st in vars + x = operation(st) + t = only(arguments(st)) + t isa Union{Num, MTK.Symbolic} && continue + idx = stidxmap[x(iv)] + jcosts[i] = Symbolics.substitute(jcosts[i], Dict(x(t) => U(t)[idx])) + end + end + jcosts = substitute_casadi_vars(model::CasADiModel, sys, pmap, jcosts; auxmap) + jcosts = map( + c -> Symbolics.substitute(c, MTK.∫() => Symbolics.Integral(iv in tspan)), jcosts) + + 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) + (lo, hi) !== tspan && error("Non-whole interval bounds for integrals are not currently supported.") + intmap[int] = dt * tₛ * sum(arg) + end + jcosts = map(c -> Symbolics.substitute(c, intmap), jcosts) + minimize!(opti, consolidate(jcosts)) +end + +function substitute_casadi_vars(model::CasADiModel, sys, pmap, exprs; auxmap = Dict()) + @unpack opti, U, V = 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) + + # 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; is_free_t) - A = tableau.A - α = tableau.α - c = tableau.c - model = prob.model - f = prob.f - p = prob.p + @unpack A, α, c = tableau + @unpack model, f, p = prob + @unpack opti, U, V, tₛ = model - opti = model.opti - t = model[:t] - tsteps = supports(t) - tmax = tsteps[end] - pop!(tsteps) - tₛ = is_free_t ? model[:tf] : 1 + tsteps = U.t dt = tsteps[2] - tsteps[1] - U = model.U - V = model.V nᵤ = length(U) nᵥ = length(V) if is_explicit(tableau) K = Any[] - for τ in tsteps + for k in 1:length(tsteps)-1 for (i, h) in enumerate(c) ΔU = sum([A[i, j] * K[j] for j in 1:(i - 1)], init = zeros(nᵤ)) - Uₙ = [U[i](τ) + ΔU[i] * dt for i in 1:nᵤ] - Vₙ = [V[i](τ) for i in 1: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!(model, U[n](τ) + ΔU[n]==U[n](τ + dt)) - empty!(K) + subject_to!(opti, U.u[:, k] + ΔU == U.u[:, k+1]) end else + ΔU_tot = dt * (K' * α) + for k in 1:length(tsteps)-1 + Kᵢ = variable!(opti, length(α), nᵤ) + ΔUs = A * Kᵢ # the stepsize at each stage of the implicit method + for (i, h) in enumerate(c) + ΔU = @view ΔUs[i, :] + Uₙ = U.u[:,k] + ΔU + Vₙ = V.u[:,k] + subject_to!(opti, K[i,:] == tₛ * f(Uₙ, Vₙ, p, τ + h*dt)) + end + ΔU_tot = dt*(Kᵢ'*α) + subject_to!(opti, U.u[:, k] + ΔU_tot == U.u[:,k+1]) + end end end -function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, Symbol}, ode_solver::Union{String, Symbol}; silent = false) +is_explicit(tableau) = tableau isa DiffEqDevTools.ExplicitRKTableau + +""" + solve(prob::CasADiDynamicOptProblem, casadi_solver, ode_solver; plugin_options, solver_options) + +`plugin_options` and `solver_options` get propagated to the Opti object in CasADi. +""" +function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, Symbol}, ode_solver::Union{String, Symbol}; plugin_options::Dict = Dict(), solver_options::Dict = Dict(), silent = false) model = prob.model opti = model.opti - solver!(opti, solver) - sol = solve(opti) - DynamicOptSolution(model, sol, input_sol) + solver!(opti, solver, plugin_options, solver_options) + add_casadi_solve_constraints!(prob, tableau) + solver!(cmodel, "$solver", plugin_options, solver_options) + + failed = false + try + sol = solve(opti) + value_getter = x -> CasADi.value(sol, x) + catch ErrorException + value_getter = x -> CasADi.debug_value(opti, x) + failed = true + continue + end + + ts = value_getter(tₛ) * U.t + U_vals = value_getter(U) + U_vals = [[U_vals[i][j] for i in 1:length(U_vals)] for j in 1:length(ts)] + sol = DiffEqBase.build_solution(prob, ode_solver, ts, U_vals) + + input_sol = nothing + if !isempty(V) + V_vals = value_getter(V) + V_vals = [[V_vals[i][j] for i in 1:length(V_vals)] for j in 1:length(ts)] + input_sol = DiffEqBase.build_solution(prob, ode_solver, ts, V_vals) + end + + if failed + sol = SciMLBase.solution_new_retcode(sol, SciMLBase.ReturnCode.ConvergenceFailure) + !isnothing(input_sol) && (input_sol = SciMLBase.solution_new_retcode( + input_sol, SciMLBase.ReturnCode.ConvergenceFailure)) + end + + DynamicOptSolution(cmodel, sol, input_sol) end end diff --git a/test/downstream/dynamic_opt_systems.jl b/test/downstream/dynamic_opt_systems.jl index 30d8feb1cd..f13b6446ef 100644 --- a/test/downstream/dynamic_opt_systems.jl +++ b/test/downstream/dynamic_opt_systems.jl @@ -1,9 +1,27 @@ -function lotkavolterra() - -end +function build_lotkavolterra(; with_constraint = false) + @parameters α=1.5 β=1.0 γ=3.0 δ=1.0 + @variables x(..) y(..) + t = M.t_nounits + D = M.D_nounits -function () - + eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y(t), + D(y(t)) ~ -γ * y(t) + δ * x(t) * y(t)] + + tspan = (0.0, 1.0) + parammap = [α => 1.5, β => 1.0, γ => 3.0, δ => 1.0] + + if with_constraint + constr = [x(0.6) ~ 3.5, x(0.3) ~ 7.0] + guess = [x(t) => 4.0, y(t) => 2.0] + u0map = Pair[] + else + constr = nothing + guess = Pair[] + u0map = [x(t) => 4.0, y(t) => 2.0] + end + + @mtkbuild sys = ODESystem(eqs, t; constraints = constr) + sys, u0map, tspan, parammap, guess end function rocket_fft() From 1d5223634c9bdfb99dbcb04822c76b8359661d85 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 3 May 2025 00:26:08 -0400 Subject: [PATCH 08/28] refactor: move tableau generation to interface --- ext/MTKInfiniteOptExt.jl | 17 -- ext/MTKJuMPControlExt.jl | 324 --------------------------------------- 2 files changed, 341 deletions(-) delete mode 100644 ext/MTKJuMPControlExt.jl diff --git a/ext/MTKInfiniteOptExt.jl b/ext/MTKInfiniteOptExt.jl index ba3a772582..5a48f9b96f 100644 --- a/ext/MTKInfiniteOptExt.jl +++ b/ext/MTKInfiniteOptExt.jl @@ -327,23 +327,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 diff --git a/ext/MTKJuMPControlExt.jl b/ext/MTKJuMPControlExt.jl deleted file mode 100644 index de3aa96269..0000000000 --- a/ext/MTKJuMPControlExt.jl +++ /dev/null @@ -1,324 +0,0 @@ -module MTKJuMPControlExt -using ModelingToolkit -using JuMP, InfiniteOpt -using DiffEqDevTools, DiffEqBase, SciMLBase -using LinearAlgebra -const MTK = ModelingToolkit - -struct JuMPControlProblem{uType, tType, P, F, K} - f::F - u0::uType - tspan::tType - p::P - model::InfiniteModel - kwargs::K - - function JuMPControlProblem(f, u0, tspan, p, model; kwargs...) - new{typeof(u0), typeof(tspan), typeof(p), typeof(f), typeof(kwargs)}(f, u0, tspan, p, model, kwargs) - end -end - -struct InfiniteOptControlProblem{uType, tType, isinplace, P, F, K} <: - AbstractOptimalControlProblem{uType, tType, isinplace} - f::F - u0::uType - tspan::tType - p::P - model::InfiniteModel - kwargs::K - - function InfiniteOptControlProblem(f, u0, tspan, p, model, kwargs...) - new{typeof(u0), typeof(tspan), SciMLBase.isinplace(f), - typeof(p), typeof(f), typeof(kwargs)}(f, u0, tspan, p, model, kwargs) - end -end - -""" - JuMPControlProblem(sys::ODESystem, u0, tspan, p; dt) - -Convert an ODESystem representing an optimal control system into a JuMP model -for solving using optimization. Must provide `dt` for determining the length -of the interpolation arrays. - -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.JuMPControlProblem(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, kwargs...) - - if !isnothing(constraintsys) - (length(constraints(constraintsys)) + length(u0map) > length(sts)) && - @warn "The JuMPControlProblem is overdetermined. The total number of conditions (# constraints + # fixed initial values given by u0map) exceeds the total number of states. The solvers will default to doing a nonlinear least-squares optimization." - end - - JuMPControlProblem(f, u0, tspan, p, model, kwargs...) -end - -""" - InfiniteOptControlProblem(sys::ODESystem, u0map, tspan, pmap; dt) - -Convert an ODESystem representing an optimal control system into a InfiniteOpt model -for solving using optimization. Must provide `dt` for determining the length -of the interpolation arrays. - -Related to `JuMPControlProblem`, but directly adds the differential equations -of the system as derivative constraints, rather than using a solver tableau. -""" -function MTK.InfiniteOptControlProblem(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, kwargs...) - - model = InfiniteModel() - @infinite_parameter(model, t in [ts, te], num_supports = length(steps), derivative_method = OrthogonalCollocation(2)) - @variable(model, U[1:length(states)], Infinite(t), start = ts) - @variable(model, V[1:length(ctrls)], Infinite(t), start = ts) - - add_jump_cost_function!(model, sys) - add_user_constraints!(model, sys) - - @infinite_parameter(model, t in [tspan[1], tspan[2]], num_supports=steps) - @variable(model, U[i = 1:length(states)], Infinite(t), start=u0[i]) - c0 = MTK.value.([pmap[c] for c in ctrls]) - @variable(model, V[i = 1:length(ctrls)], Infinite(t), start=c0[i]) - - set_jump_bounds!(model, sys, pmap) - add_jump_cost_function!(model, sys, (tspan[1], tspan[2]), pmap; is_free_t) - add_user_constraints!(model, sys, 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, tspan[1]) - return model -end - - JuMPControlProblem(f, u0, tspan, p, model, kwargs...) -end - -function add_jump_cost_function!(model, sys) - jcosts = MTK.get_costs(sys) - consolidate = MTK.get_consolidate(sys) - if isnothing(consolidate) - @objective(model, Min, 0) - return - end - iv = MTK.get_iv(sys) - - stidxmap = Dict([v => i for (i, v) in enumerate(unknowns(sys))]) - cidxmap = Dict([v => i for (i, v) in enumerate(controls(sys))]) - - for st in unknowns(sys) - x = operation(st) - t = only(arguments(st)) - idx = stidxmap[x(iv)] - subval = isequal(t, iv) ? model[:U][idx] : model[:U][idx](t) - jcosts = Symbolics.substitute(jcosts, Dict(x(t) => subval)) - end - - for ct in controls(sys) - p = operation(ct) - t = only(arguments(ct)) - idx = cidxmap[p(iv)] - subval = isequal(t, iv) ? model[:V][idx] : model[:V][idx](t) - jcosts = Symbolics.substitute(jcosts, Dict(x(t) => subval)) - end - - @objective(model, Min, consolidate(jcosts)) -end - -function add_user_constraints!(model, sys) - jconstraints = if !(csys = MTK.get_constraintsystem(sys) isa Nothing) - MTK.get_constraints(csys) - else - nothing - end - isnothing(jconstraints) && return nothing - - iv = MTK.get_iv(sys) - stidxmap = Dict([v => i for (i, v) in enumerate(unknowns(sys))]) - cidxmap = Dict([v => i for (i, v) in enumerate(controls(sys))]) - - for st in unknowns(sys) - x = operation(st) - t = only(arguments(st)) - idx = stidxmap[x(iv)] - subval = isequal(t, iv) ? model[:U][idx] : model[:U][idx](t) - jconstraints = Symbolics.substitute(jconstraints, Dict(x(t) => subval)) - end - - auxmap = Dict([u => MTK.default_toterm(MTK.value(u)) for u in unknowns(conssys)]) - jconstraints = substitute_jump_vars(model, sys, pmap, jconstraints; auxmap) - - # Substitute to-term'd variables - for (i, cons) in enumerate(jconstraints) - if cons isa Equation - @constraint(model, user[i], cons.lhs - cons.rhs == 0) - elseif cons.relational_op === Symbolics.geq - @constraint(model, user[i], cons.lhs - cons.rhs ≥ 0) - else - @constraint(model, user[i], cons.lhs - cons.rhs ≤ 0) - end - end -end - -function add_initial_constraints!(model, u0, u0_idxs, tspan) - ts = tspan[1] - @constraint(model, init_u0_idx[i in u0_idxs], model[:U][i](ts) == u0[i]) -end - -function substitute_jump_vars(model, sys, pmap, exprs; auxmap = Dict()) - iv = MTK.get_iv(sys) - sts = unknowns(sys) - cts = MTK.unbound_inputs(sys) - U = model[:U] - V = model[:V] - exprs = map(c -> Symbolics.fixpoint_sub(c, auxmap), exprs) - - # for variables like x(t) - whole_interval_map = Dict([[v => U[i] for (i, v) in enumerate(sts)]; - [v => V[i] for (i, v) in enumerate(cts)]]) - exprs = map(c -> Symbolics.fixpoint_sub(c, whole_interval_map), exprs) - - # for variables like x(1.0) - x_ops = [MTK.operation(MTK.unwrap(st)) for st in sts] - c_ops = [MTK.operation(MTK.unwrap(ct)) for ct in cts] - fixed_t_map = Dict([[x_ops[i] => U[i] for i in 1:length(U)]; - [c_ops[i] => V[i] for i in 1:length(V)]]) - - exprs = map(c -> Symbolics.fixpoint_sub(c, fixed_t_map), exprs) - - exprs = map(c -> Symbolics.fixpoint_sub(c, Dict(pmap)), exprs) - exprs -end - -is_explicit(tableau) = tableau isa DiffEqDevTools.ExplicitRKTableau - -function add_solve_constraints!(prob, tableau) - A = tableau.A - α = tableau.α - c = tableau.c - model = prob.model - f = prob.f - p = prob.p - - t = model[:t] - tsteps = supports(t) - tmax = tsteps[end] - pop!(tsteps) - dt = tsteps[2] - tsteps[1] - - U = model[:U] - nᵤ = length(U) - if is_explicit(tableau) - K = Any[] - for τ in tsteps - for (i, h) in enumerate(c) - ΔU = sum([A[i, j] * K[j] for j in 1:i-1], init = zeros(nᵤ)) - Uₙ = [U[i](τ) + ΔU[i]*dt for i in 1:nᵤ] - Kₙ = f(Uₙ, p, τ + h*dt) - push!(K, Kₙ) - end - ΔU = sum([α[i] * K[i] for i in 1:length(α)]) - @constraint(model, [n = 1:nᵤ], U[n](τ) + ΔU[n] == U[n](τ + dt)) - empty!(K) - end - else - @variable(model, K[1:length(α), 1:nᵤ], Infinite(t)) - ΔUs = A * K - ΔU_tot = dt * (K' * α) - for τ in tsteps - ΔUs = [A * K(τ)] - for (i, h) in enumerate(c) - ΔU = @view ΔUs[i, :] - Uₙ = U + ΔU * h * dt - @constraint(model, [j = 1:nᵤ], K[i, j]==(tₛ * f(Uₙ, V, p, τ + h * dt)[j]), - DomainRestrictions(t => τ), base_name="solve_K$i($τ)") - end - @constraint(model, [n = 1:nᵤ], U[n](τ) + ΔU_tot[n]==U[n](min(τ + dt, tmax)), - DomainRestrictions(t => τ), base_name="solve_U($τ)") - end - end -end - -""" -""" -struct JuMPControlSolution - model::InfiniteModel - sol::ODESolution -end - -""" -Solve JuMPControlProblem. Arguments: -- prob: a JumpControlProblem -- jump_solver: a LP solver such as HiGHS -- ode_solver: Takes in a symbol representing the solver. Acceptable solvers may be found at https://docs.sciml.ai/DiffEqDevDocs/stable/internals/tableaus/. Note that the symbol may be different than the typical name of the solver, e.g. :Tsitouras5 rather than Tsit5. -""" -function DiffEqBase.solve(prob::JuMPControlProblem, jump_solver, ode_solver::Symbol) - model = prob.model - tableau_getter = Symbol(:construct, ode_solver) - @eval tableau = $tableau_getter() - ts = supports(model[:t]) - add_solve_constraints!(prob, tableau) - - # Unregister current solver constraints - for con in all_constraints(model) - if occursin("solve", JuMP.name(con)) - unregister(model, Symbol(JuMP.name(con))) - delete(model, con) - end - end - unregister(model, :K) - for var in all_variables(model) - if occursin("K", JuMP.name(var)) - unregister(model, Symbol(JuMP.name(var))) - delete(model, var) - end - end - add_jump_solve_constraints!(prob, tableau; is_free_t = haskey(model, :tf)) - _solve(prob, jump_solver, ode_solver) -end - -""" -`derivative_method` kwarg refers to the method used by InfiniteOpt to compute derivatives. The list of possible options can be found at https://infiniteopt.github.io/InfiniteOpt.jl/stable/guide/derivative/. Defaults to FiniteDifference(Backward()). -""" -function DiffEqBase.solve(prob::InfiniteOptControlProblem, jump_solver; - derivative_method = InfiniteOpt.FiniteDifference(Backward()), silent = false) - model = prob.model - silent && set_silent(model) - set_derivative_method(model[:t], derivative_method) - _solve(prob, jump_solver, derivative_method) -end - -function _solve(prob::AbstractOptimalControlProblem, jump_solver, solver) - model = prob.model - set_optimizer(model, jump_solver) - optimize!(model) - - if is_solved_and_feasible(model) - sol = DiffEqBase.build_solution(prob, ode_solver, ts, value.(U)) - JuMPControlSolution(model, sol) - else - sol = DiffEqBase.build_solution(prob, ode_solver, ts, value.(U)) - sol = SciMLBase.solution_new_retcode(sol, SciMLBase.ReturnCode.ConvergenceFailure) - JuMPControlSolution(model, sol) - end -end - -end From b863d5dc28f48ca22f1d9ba5a8f5479fefc0f1de Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 3 May 2025 00:26:45 -0400 Subject: [PATCH 09/28] rebase --- Project.toml | 2 +- ext/MTKCasADiDynamicOptExt.jl | 13 +++++++++---- src/systems/optimal_control_interface.jl | 20 ++++++++++++++++++++ 3 files changed, 30 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 2ea71a8ed3..060a54c1f8 100644 --- a/Project.toml +++ b/Project.toml @@ -77,7 +77,7 @@ LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" [extensions] MTKBifurcationKitExt = "BifurcationKit" -MTKCasADiDynamicOptExt = ["CasADi", "DataInterpolations", "DiffEqDevTools"] +MTKCasADiDynamicOptExt = ["CasADi", "DataInterpolations"] MTKChainRulesCoreExt = "ChainRulesCore" MTKDeepDiffsExt = "DeepDiffs" MTKFMIExt = "FMI" diff --git a/ext/MTKCasADiDynamicOptExt.jl b/ext/MTKCasADiDynamicOptExt.jl index e16ce1b6b5..bc04111714 100644 --- a/ext/MTKCasADiDynamicOptExt.jl +++ b/ext/MTKCasADiDynamicOptExt.jl @@ -1,10 +1,10 @@ module MTKCasADiDynamicOptExt using ModelingToolkit using CasADi -using DiffEqDevTools, DiffEqBase +using DiffEqBase using DataInterpolations using UnPack -const MTK = MOdelingToolkit +const MTK = ModelingToolkit struct CasADiDynamicOptProblem{uType, tType, isinplace, P, F, K} <: AbstractDynamicOptProblem{uType, tType, isinplace} @@ -221,17 +221,22 @@ function add_solve_constraints!(prob, tableau; is_free_t) end end -is_explicit(tableau) = tableau isa DiffEqDevTools.ExplicitRKTableau """ solve(prob::CasADiDynamicOptProblem, casadi_solver, ode_solver; plugin_options, solver_options) `plugin_options` and `solver_options` get propagated to the Opti object in CasADi. """ -function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, Symbol}, ode_solver::Union{String, Symbol}; plugin_options::Dict = Dict(), solver_options::Dict = Dict(), silent = false) +function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, Symbol}, ode_solver::Symbol = :Default; plugin_options::Dict = Dict(), solver_options::Dict = Dict(), silent = false) model = prob.model opti = model.opti + if ode_solver == :Default + tableau = MTK.constructDefault() + else + tableau_getter = Symbol(:construct, ode_solver) + tableau = @eval Main.tableau_getter() + end solver!(opti, solver, plugin_options, solver_options) add_casadi_solve_constraints!(prob, tableau) solver!(cmodel, "$solver", plugin_options, solver_options) 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. From 8db9a75d5adf00bbe09f02de354376dbbc484036 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 3 May 2025 01:42:20 -0400 Subject: [PATCH 10/28] apply stash --- ext/MTKCasADiDynamicOptExt.jl | 1 - src/ModelingToolkit.jl | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/ext/MTKCasADiDynamicOptExt.jl b/ext/MTKCasADiDynamicOptExt.jl index bc04111714..44f6d254ea 100644 --- a/ext/MTKCasADiDynamicOptExt.jl +++ b/ext/MTKCasADiDynamicOptExt.jl @@ -2,7 +2,6 @@ module MTKCasADiDynamicOptExt using ModelingToolkit using CasADi using DiffEqBase -using DataInterpolations using UnPack const MTK = ModelingToolkit 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 From e81e9c4beee98e770c3627e6e89da450eb1cac86 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 3 May 2025 01:42:42 -0400 Subject: [PATCH 11/28] drop DataInteprolations --- Project.toml | 7 +------ ext/MTKCasADiDynamicOptExt.jl | 34 +++++++++++++++++++++++----------- 2 files changed, 24 insertions(+), 17 deletions(-) diff --git a/Project.toml b/Project.toml index 060a54c1f8..33280a5064 100644 --- a/Project.toml +++ b/Project.toml @@ -31,7 +31,6 @@ FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e" FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5" Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" @@ -67,17 +66,14 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665" CasADi = "c49709b8-5c63-11e9-2fb2-69db5844192f" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" -DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" FMI = "14a09403-18e3-468f-ad8a-74f8dda2d9ac" InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" -JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" [extensions] MTKBifurcationKitExt = "BifurcationKit" -MTKCasADiDynamicOptExt = ["CasADi", "DataInterpolations"] +MTKCasADiDynamicOptExt = "CasADi" MTKChainRulesCoreExt = "ChainRulesCore" MTKDeepDiffsExt = "DeepDiffs" MTKFMIExt = "FMI" @@ -104,7 +100,6 @@ DeepDiffs = "1" DelayDiffEq = "5.50" DiffEqBase = "6.170.1" DiffEqCallbacks = "2.16, 3, 4" -DiffEqDevTools = "2.48.0" DiffEqNoiseProcess = "5" DiffRules = "0.1, 1.0" DifferentiationInterface = "0.6.47" diff --git a/ext/MTKCasADiDynamicOptExt.jl b/ext/MTKCasADiDynamicOptExt.jl index 44f6d254ea..a4cfcbd779 100644 --- a/ext/MTKCasADiDynamicOptExt.jl +++ b/ext/MTKCasADiDynamicOptExt.jl @@ -20,17 +20,32 @@ struct CasADiDynamicOptProblem{uType, tType, isinplace, P, F, K} <: end 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::AbstractInterpolation - V::AbstractInterpolation + U::MXLinearInterpolation + V::MXLinearInterpolation tₛ::Union{Number, MX} end +function (M::MXLinearInterpolation)(τ) + nt = (τ - M.t) / M.dt + i = 1 + floor(Int, nt) + Δ = nt - i + 1 + + (i > length(M.t) || i < 1) && error("Cannot extrapolate past the tspan.") + M.u[i] + Δ*(M.u[i + 1] - M.u[i]) +end + function MTK.CasADiDynamicOptProblem(sys::ODESystem, u0map, tspan, pmap; dt = nothing, steps = nothing, - interpolation_method::AbstractInterpolation = LinearInterpolation, guesses = Dict(), kwargs...) MTK.warn_overdetermined(sys, u0map) _u0map = has_alg_eqs(sys) ? u0map : merge(Dict(u0map), Dict(guesses)) @@ -39,10 +54,10 @@ function MTK.CasADiDynamicOptProblem(sys::ODESystem, u0map, tspan, pmap; pmap = Dict{Any, Any}(pmap) steps, is_free_t = MTK.process_tspan(tspan, dt, steps) - model = init_model() + model = init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t) end -function init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t = false, interp_type::AbstractInterpolation) +function init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t = false) ctrls = MTK.unbound_inputs(sys) states = unknowns(sys) opti = CasADi.Opti() @@ -61,8 +76,8 @@ function init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t = false, inter U = CasADi.variable!(opti, length(states), steps) V = CasADi.variable!(opti, length(ctrls), steps) - U_interp = interp_type(Matrix(U), tsteps) - V_interp = interp_type(Matrix(V), tsteps) + U_interp = MXLinearInterpolation(U, tsteps, tsteps[2]-tsteps[1]) + V_interp = MXLinearInterpolation(V, tsteps, tsteps[2]-tsteps[1]) CasADiModel(opti, U_interp, V_interp, tₛ) end @@ -220,9 +235,8 @@ function add_solve_constraints!(prob, tableau; is_free_t) end end - """ - solve(prob::CasADiDynamicOptProblem, casadi_solver, ode_solver; plugin_options, solver_options) + 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. """ @@ -247,7 +261,6 @@ function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, S catch ErrorException value_getter = x -> CasADi.debug_value(opti, x) failed = true - continue end ts = value_getter(tₛ) * U.t @@ -270,5 +283,4 @@ function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, S DynamicOptSolution(cmodel, sol, input_sol) end - end From 280224403591da99e9327a60b925cbd80feac9d9 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 5 May 2025 10:33:39 -0400 Subject: [PATCH 12/28] fix: don't use eval --- ext/MTKCasADiDynamicOptExt.jl | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/ext/MTKCasADiDynamicOptExt.jl b/ext/MTKCasADiDynamicOptExt.jl index a4cfcbd779..4fea4a94e0 100644 --- a/ext/MTKCasADiDynamicOptExt.jl +++ b/ext/MTKCasADiDynamicOptExt.jl @@ -43,6 +43,22 @@ function (M::MXLinearInterpolation)(τ) M.u[i] + Δ*(M.u[i + 1] - M.u[i]) 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, @@ -240,16 +256,11 @@ end `plugin_options` and `solver_options` get propagated to the Opti object in CasADi. """ -function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, Symbol}, ode_solver::Symbol = :Default; plugin_options::Dict = Dict(), solver_options::Dict = Dict(), silent = false) +function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, Symbol}, tableau_getter = constructDefault; plugin_options::Dict = Dict(), solver_options::Dict = Dict(), silent = false) model = prob.model + tableau = tableau_getter() opti = model.opti - if ode_solver == :Default - tableau = MTK.constructDefault() - else - tableau_getter = Symbol(:construct, ode_solver) - tableau = @eval Main.tableau_getter() - end solver!(opti, solver, plugin_options, solver_options) add_casadi_solve_constraints!(prob, tableau) solver!(cmodel, "$solver", plugin_options, solver_options) @@ -266,13 +277,13 @@ function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, S ts = value_getter(tₛ) * U.t U_vals = value_getter(U) U_vals = [[U_vals[i][j] for i in 1:length(U_vals)] for j in 1:length(ts)] - sol = DiffEqBase.build_solution(prob, ode_solver, ts, U_vals) + sol = DiffEqBase.build_solution(prob, tableau_getter, ts, U_vals) input_sol = nothing if !isempty(V) V_vals = value_getter(V) V_vals = [[V_vals[i][j] for i in 1:length(V_vals)] for j in 1:length(ts)] - input_sol = DiffEqBase.build_solution(prob, ode_solver, ts, V_vals) + input_sol = DiffEqBase.build_solution(prob, tableau_getter, ts, V_vals) end if failed From 341d1757e9028a7644a7b7b7744563d59dec250a Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 5 May 2025 17:19:15 -0400 Subject: [PATCH 13/28] feat: ode transcription working --- .CondaPkg/.gitattributes | 2 + .CondaPkg/.gitignore | 4 + .CondaPkg/meta | Bin 0 -> 624 bytes .CondaPkg/pixi.lock | 540 ++++++++++++++++++++++++++++++++++ .CondaPkg/pixi.toml | 15 + Project.toml | 2 +- ext/MTKCasADiDynamicOptExt.jl | 139 +++++---- 7 files changed, 640 insertions(+), 62 deletions(-) create mode 100644 .CondaPkg/.gitattributes create mode 100644 .CondaPkg/.gitignore create mode 100644 .CondaPkg/meta create mode 100644 .CondaPkg/pixi.lock create mode 100644 .CondaPkg/pixi.toml diff --git a/.CondaPkg/.gitattributes b/.CondaPkg/.gitattributes new file mode 100644 index 0000000000..887a2c18f0 --- /dev/null +++ b/.CondaPkg/.gitattributes @@ -0,0 +1,2 @@ +# SCM syntax highlighting & preventing 3-way merges +pixi.lock merge=binary linguist-language=YAML linguist-generated=true diff --git a/.CondaPkg/.gitignore b/.CondaPkg/.gitignore new file mode 100644 index 0000000000..740bb7d1ae --- /dev/null +++ b/.CondaPkg/.gitignore @@ -0,0 +1,4 @@ + +# pixi environments +.pixi +*.egg-info diff --git a/.CondaPkg/meta b/.CondaPkg/meta new file mode 100644 index 0000000000000000000000000000000000000000..f82e7e2daa6bea42c0dfd7469940ff53a33ba22d GIT binary patch literal 624 zcmb7BOHRWu5DmWygv1TFKxsWqQcxFAsR|nwm9hX;T@JCEx^ZpGb_!jv=Zsu|OOQ4W zLWnHz#`Ae^-uU?UwL>p(*)Fs*EI~WD=e>52 z#;m}cSxC2Tsbqpez-&634?3PRAQ6d1$39Cdmm6anM1~eAZxG{{G>&^t5S;i(Z`E3T tSAPY~IK5xw)OW{sF&Xu4hvz=eb|2nfD3h}@U+QKxrF-BDe_(wl_yTx^l-mFR literal 0 HcmV?d00001 diff --git a/.CondaPkg/pixi.lock b/.CondaPkg/pixi.lock new file mode 100644 index 0000000000..43d038b1dd --- /dev/null +++ b/.CondaPkg/pixi.lock @@ -0,0 +1,540 @@ +version: 6 +environments: + default: + channels: + - url: https://conda.anaconda.org/conda-forge/ + packages: + osx-arm64: + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/ampl-asl-1.0.0-h286801f_2.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/bzip2-1.0.8-h99b78c6_7.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/ca-certificates-2025.4.26-hbd8a1cb_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/casadi-3.7.0-py312h74b190f_2.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/ipopt-3.14.17-h945cc1c_2.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libblas-3.9.0-31_h10e41b3_openblas.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libblasfeo-0.1.4.2-h37ef02a_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libcblas-3.9.0-31_hb3479ef_openblas.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libcxx-20.1.4-ha82da77_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libexpat-2.7.0-h286801f_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libfatrop-0.0.4-h286801f_1.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libffi-3.4.6-h1da3d7d_1.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libgfortran-5.0.0-14_2_0_h6c33f7e_103.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libgfortran5-14.2.0-h6c33f7e_103.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/liblapack-3.9.0-31_hc9a63f6_openblas.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/liblzma-5.8.1-h39f12f2_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libopenblas-0.3.29-openmp_hf332438_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libosqp-0.6.3-h5833ebf_1.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libqdldl-0.1.7-hb7217d7_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libscotch-7.0.6-he56f69b_1.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libsqlite-3.49.1-h3f77e49_2.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libzlib-1.3.1-h8359307_2.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/llvm-openmp-20.1.4-hdb05f8b_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/metis-5.1.0-h15f6cfe_1007.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/mumps-include-5.7.3-h8c5b6c6_10.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/mumps-seq-5.7.3-h390d176_10.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/ncurses-6.5-h5e97a16_3.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numpy-2.2.5-py312h7c1f314_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/openssl-3.5.0-h81ee809_1.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/ply-3.11-pyhd8ed1ab_3.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/pyomo-6.9.2-py312hd8f9ff3_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/python-3.12.10-hc22306f_0_cpython.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.12-7_cp312.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/readline-8.2-h1d1bf99_2.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.1.0-pyhff2d567_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/tinyxml2-11.0.0-ha1acc90_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/tk-8.6.13-h5083fa2_1.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/tzdata-2025b-h78e105d_0.conda +packages: +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/ampl-asl-1.0.0-h286801f_2.conda + sha256: 47a5da6cd38a32c086017a5d2ed2b67e0cc24e25268a9c4cc91ea7d8f1cb9507 + md5: bb25b8fa2b28474412dda4e1a95853b4 + depends: + - __osx >=11.0 + - libcxx >=18 + constrains: + - ampl-mp >=4.0.0 + arch: arm64 + platform: osx + license: BSD-3-Clause AND SMLNJ + size: 411989 + timestamp: 1732439500161 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/bzip2-1.0.8-h99b78c6_7.conda + sha256: adfa71f158cbd872a36394c56c3568e6034aa55c623634b37a4836bd036e6b91 + md5: fc6948412dbbbe9a4c9ddbbcfe0a79ab + depends: + - __osx >=11.0 + arch: arm64 + platform: osx + license: bzip2-1.0.6 + license_family: BSD + size: 122909 + timestamp: 1720974522888 +- conda: https://conda.anaconda.org/conda-forge/noarch/ca-certificates-2025.4.26-hbd8a1cb_0.conda + sha256: 2a70ed95ace8a3f8a29e6cd1476a943df294a7111dfb3e152e3478c4c889b7ac + md5: 95db94f75ba080a22eb623590993167b + depends: + - __unix + license: ISC + size: 152283 + timestamp: 1745653616541 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/casadi-3.7.0-py312h74b190f_2.conda + sha256: ae008bb0cd026b64d6cbcd7726fbbbb98904bd868fca3b7d53f258f0c80c27e7 + md5: 8431af5d6849a519573d8ed6ffaa8e28 + depends: + - __osx >=11.0 + - ipopt >=3.14.17,<3.14.18.0a0 + - libblas >=3.9.0,<4.0a0 + - libblasfeo >=0.1.4.2,<0.1.5.0a0 + - libcblas >=3.9.0,<4.0a0 + - libcxx >=18 + - libfatrop >=0.0.4,<0.0.5.0a0 + - liblapack >=3.9.0,<4.0a0 + - libosqp >=0.6.3,<0.6.4.0a0 + - numpy >=1.19,<3 + - python >=3.12,<3.13.0a0 + - python >=3.12,<3.13.0a0 *_cpython + - python_abi 3.12.* *_cp312 + - tinyxml2 >=11.0.0,<11.1.0a0 + arch: arm64 + platform: osx + license: LGPL-3.0-or-later + license_family: LGPL + size: 4681272 + timestamp: 1745474347317 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/ipopt-3.14.17-h945cc1c_2.conda + sha256: 379d9b44fc265493a8a42ca3694bb8d7538ed8c7bbaf8626dcf8d75f454317cc + md5: e465f880544e1e6b8d904fb65e25b130 + depends: + - __osx >=11.0 + - ampl-asl >=1.0.0,<1.0.1.0a0 + - libblas >=3.9.0,<4.0a0 + - libcxx >=18 + - liblapack >=3.9.0,<4.0a0 + - mumps-seq >=5.7.3,<5.7.4.0a0 + arch: arm64 + platform: osx + license: EPL-1.0 + size: 727718 + timestamp: 1745419955303 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libblas-3.9.0-31_h10e41b3_openblas.conda + build_number: 31 + sha256: 369586e7688b59b4f92c709b99d847d66d4d095425db327dd32ee5e6ab74697f + md5: 39b053da5e7035c6592102280aa7612a + depends: + - libopenblas >=0.3.29,<0.3.30.0a0 + - libopenblas >=0.3.29,<1.0a0 + constrains: + - liblapacke =3.9.0=31*_openblas + - libcblas =3.9.0=31*_openblas + - blas =2.131=openblas + - mkl <2025 + - liblapack =3.9.0=31*_openblas + arch: arm64 + platform: osx + license: BSD-3-Clause + license_family: BSD + size: 17123 + timestamp: 1740088119350 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libblasfeo-0.1.4.2-h37ef02a_0.conda + sha256: 406885b43c417025be9b2cdb541fa5a0f5426efbbf109d42ac5126efb5d3608d + md5: db8bf9db5c73a85b195f5c973f5514bf + depends: + - __osx >=11.0 + arch: arm64 + platform: osx + license: BSD-2-Clause + license_family: BSD + size: 509787 + timestamp: 1740067975938 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libcblas-3.9.0-31_hb3479ef_openblas.conda + build_number: 31 + sha256: f237486cc9118d09d0f3ff8820280de34365f98ee7b7dc5ab923b04c7cbf25a5 + md5: 7353c2bf0e90834cb70545671996d871 + depends: + - libblas 3.9.0 31_h10e41b3_openblas + constrains: + - liblapacke =3.9.0=31*_openblas + - blas =2.131=openblas + - liblapack =3.9.0=31*_openblas + arch: arm64 + platform: osx + license: BSD-3-Clause + license_family: BSD + size: 17032 + timestamp: 1740088127097 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libcxx-20.1.4-ha82da77_0.conda + sha256: 1837e2c65f8fc8cfd8b240cfe89406d0ce83112ac63f98c9fb3c9a15b4f2d4e1 + md5: 10c809af502fcdab799082d338170994 + depends: + - __osx >=11.0 + arch: arm64 + platform: osx + license: Apache-2.0 WITH LLVM-exception + license_family: Apache + size: 565811 + timestamp: 1745991653948 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libexpat-2.7.0-h286801f_0.conda + sha256: ee550e44765a7bbcb2a0216c063dcd53ac914a7be5386dd0554bd06e6be61840 + md5: 6934bbb74380e045741eb8637641a65b + depends: + - __osx >=11.0 + constrains: + - expat 2.7.0.* + arch: arm64 + platform: osx + license: MIT + license_family: MIT + size: 65714 + timestamp: 1743431789879 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libfatrop-0.0.4-h286801f_1.conda + sha256: ff193f966ee33546a297abe02e49b57ba2b728d4706cec9cf3898dc8ca81e3ee + md5: 00cb2b510ea71f73263e59e8253d0720 + depends: + - __osx >=11.0 + - libblasfeo >=0.1.4.1,<0.1.5.0a0 + - libcxx >=18 + arch: arm64 + platform: osx + license: LGPL-3.0-or-later + license_family: LGPL + size: 195859 + timestamp: 1736273349013 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libffi-3.4.6-h1da3d7d_1.conda + sha256: c6a530924a9b14e193ea9adfe92843de2a806d1b7dbfd341546ece9653129e60 + md5: c215a60c2935b517dcda8cad4705734d + depends: + - __osx >=11.0 + arch: arm64 + platform: osx + license: MIT + license_family: MIT + size: 39839 + timestamp: 1743434670405 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libgfortran-5.0.0-14_2_0_h6c33f7e_103.conda + sha256: 8628746a8ecd311f1c0d14bb4f527c18686251538f7164982ccbe3b772de58b5 + md5: 044a210bc1d5b8367857755665157413 + depends: + - libgfortran5 14.2.0 h6c33f7e_103 + arch: arm64 + platform: osx + license: GPL-3.0-only WITH GCC-exception-3.1 + license_family: GPL + size: 156291 + timestamp: 1743863532821 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libgfortran5-14.2.0-h6c33f7e_103.conda + sha256: 8599453990bd3a449013f5fa3d72302f1c68f0680622d419c3f751ff49f01f17 + md5: 69806c1e957069f1d515830dcc9f6cbb + depends: + - llvm-openmp >=8.0.0 + constrains: + - libgfortran 5.0.0 14_2_0_*_103 + arch: arm64 + platform: osx + license: GPL-3.0-only WITH GCC-exception-3.1 + license_family: GPL + size: 806566 + timestamp: 1743863491726 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/liblapack-3.9.0-31_hc9a63f6_openblas.conda + build_number: 31 + sha256: fe55b9aaf82c6c0192c3d1fcc9b8e884f97492dda9a8de5dae29334b3135fab5 + md5: ff57a55a2cbce171ef5707fb463caf19 + depends: + - libblas 3.9.0 31_h10e41b3_openblas + constrains: + - liblapacke =3.9.0=31*_openblas + - libcblas =3.9.0=31*_openblas + - blas =2.131=openblas + arch: arm64 + platform: osx + license: BSD-3-Clause + license_family: BSD + size: 17033 + timestamp: 1740088134988 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/liblzma-5.8.1-h39f12f2_0.conda + sha256: 4291dde55ebe9868491dc29716b84ac3de21b8084cbd4d05c9eea79d206b8ab7 + md5: ba24e6f25225fea3d5b6912e2ac562f8 + depends: + - __osx >=11.0 + arch: arm64 + platform: osx + license: 0BSD + size: 92295 + timestamp: 1743771392206 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libopenblas-0.3.29-openmp_hf332438_0.conda + sha256: 8989d9e01ec8c9b2d48dbb5efbe70b356fcd15990fb53b64fcb84798982c0343 + md5: 0cd1148c68f09027ee0b0f0179f77c30 + depends: + - __osx >=11.0 + - libgfortran >=5 + - libgfortran5 >=13.2.0 + - llvm-openmp >=18.1.8 + constrains: + - openblas >=0.3.29,<0.3.30.0a0 + arch: arm64 + platform: osx + license: BSD-3-Clause + license_family: BSD + size: 4168442 + timestamp: 1739825514918 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libosqp-0.6.3-h5833ebf_1.conda + sha256: 05740cc58df7a9c966a277331dea1987db485508c0a35d2ac9ef562f8061df06 + md5: 7b4c6d7f35e9fe945e810e213d21b942 + depends: + - __osx >=11.0 + - libcxx >=17 + - libqdldl >=0.1.7,<0.1.8.0a0 + arch: arm64 + platform: osx + license: Apache-2.0 + license_family: APACHE + size: 360361 + timestamp: 1729342406705 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libqdldl-0.1.7-hb7217d7_0.conda + sha256: 39bb718bca8ce87afdc592d857944ba39bfac54c31a1de4d669597b52bc668ea + md5: b4df6812b4ab33c1a520287f46d91220 + depends: + - libcxx >=14.0.6 + arch: arm64 + platform: osx + license: Apache-2.0 + license_family: APACHE + size: 17189 + timestamp: 1680741853527 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libscotch-7.0.6-he56f69b_1.conda + sha256: b616a5943501e35f32699b5f29d4216a3a94bfca87b9ae30feda89a9fa5c2f7a + md5: 291ffb6a18b7e35bf42a317cc6856348 + depends: + - __osx >=11.0 + - bzip2 >=1.0.8,<2.0a0 + - libgfortran >=5 + - libgfortran5 >=13.2.0 + - liblzma >=5.6.3,<6.0a0 + - libzlib >=1.3.1,<2.0a0 + arch: arm64 + platform: osx + license: CECILL-C + size: 273604 + timestamp: 1737537249207 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libsqlite-3.49.1-h3f77e49_2.conda + sha256: 907a95f73623c343fc14785cbfefcb7a6b4f2bcf9294fcb295c121611c3a590d + md5: 3b1e330d775170ac46dff9a94c253bd0 + depends: + - __osx >=11.0 + - libzlib >=1.3.1,<2.0a0 + arch: arm64 + platform: osx + license: Unlicense + size: 900188 + timestamp: 1742083865246 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libzlib-1.3.1-h8359307_2.conda + sha256: ce34669eadaba351cd54910743e6a2261b67009624dbc7daeeafdef93616711b + md5: 369964e85dc26bfe78f41399b366c435 + depends: + - __osx >=11.0 + constrains: + - zlib 1.3.1 *_2 + arch: arm64 + platform: osx + license: Zlib + license_family: Other + size: 46438 + timestamp: 1727963202283 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/llvm-openmp-20.1.4-hdb05f8b_0.conda + sha256: b8e8547116dba85890d7b39bfad1c86ed69a6b923caed1e449c90850d271d4d5 + md5: 00cbae3f2127efef6db76bd423a09807 + depends: + - __osx >=11.0 + constrains: + - openmp 20.1.4|20.1.4.* + arch: arm64 + platform: osx + license: Apache-2.0 WITH LLVM-exception + size: 282599 + timestamp: 1746134861758 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/metis-5.1.0-h15f6cfe_1007.conda + sha256: f54ad3e5d47a0235ba2830848fee590faad550639336fe1e2413ab16fee7ac39 + md5: 7687ec5796288536947bf616179726d8 + depends: + - __osx >=11.0 + arch: arm64 + platform: osx + license: Apache-2.0 + license_family: APACHE + size: 3898314 + timestamp: 1728064659078 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/mumps-include-5.7.3-h8c5b6c6_10.conda + sha256: 9093f47e4a15fd34c90de93ebfa30c63799e2f99123420fd68653274f5574131 + md5: 077b856b81445e47942980465ca52450 + arch: arm64 + platform: osx + license: CECILL-C + size: 20799 + timestamp: 1745406952431 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/mumps-seq-5.7.3-h390d176_10.conda + sha256: 2afb273e8268a0960eb3934edffecfa617c941cccc762a11c86d15cb4f0e17a8 + md5: d211ba253191e04aa341a02203ed42fb + depends: + - mumps-include ==5.7.3 h8c5b6c6_10 + - __osx >=11.0 + - libgfortran 5.* + - libgfortran5 >=13.3.0 + - llvm-openmp >=18.1.8 + - liblapack >=3.9.0,<4.0a0 + - libblas >=3.9.0,<4.0a0 + - metis >=5.1.0,<5.1.1.0a0 + - libscotch >=7.0.6,<7.0.7.0a0 + constrains: + - libopenblas * *openmp* + arch: arm64 + platform: osx + license: CECILL-C + size: 2675502 + timestamp: 1745406952432 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/ncurses-6.5-h5e97a16_3.conda + sha256: 2827ada40e8d9ca69a153a45f7fd14f32b2ead7045d3bbb5d10964898fe65733 + md5: 068d497125e4bf8a66bf707254fff5ae + depends: + - __osx >=11.0 + arch: arm64 + platform: osx + license: X11 AND BSD-3-Clause + size: 797030 + timestamp: 1738196177597 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/numpy-2.2.5-py312h7c1f314_0.conda + sha256: 982aed7df71ae0ca8bc396ae25d43fd9a4f2b15d18faca15d6c27e9efb3955be + md5: 24a41dacf9d624b069d54a6e92594540 + depends: + - __osx >=11.0 + - libblas >=3.9.0,<4.0a0 + - libcblas >=3.9.0,<4.0a0 + - libcxx >=18 + - liblapack >=3.9.0,<4.0a0 + - python >=3.12,<3.13.0a0 + - python >=3.12,<3.13.0a0 *_cpython + - python_abi 3.12.* *_cp312 + constrains: + - numpy-base <0a0 + arch: arm64 + platform: osx + license: BSD-3-Clause + license_family: BSD + size: 6498553 + timestamp: 1745119367238 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/openssl-3.5.0-h81ee809_1.conda + sha256: 73d366c1597a10bcd5f3604b5f0734b31c23225536e03782c6a13f9be9d01bff + md5: 5c7aef00ef60738a14e0e612cfc5bcde + depends: + - __osx >=11.0 + - ca-certificates + arch: arm64 + platform: osx + license: Apache-2.0 + license_family: Apache + size: 3064197 + timestamp: 1746223530698 +- conda: https://conda.anaconda.org/conda-forge/noarch/ply-3.11-pyhd8ed1ab_3.conda + sha256: bae453e5cecf19cab23c2e8929c6e30f4866d996a8058be16c797ed4b935461f + md5: fd5062942bfa1b0bd5e0d2a4397b099e + depends: + - python >=3.9 + license: BSD-3-Clause + license_family: BSD + size: 49052 + timestamp: 1733239818090 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/pyomo-6.9.2-py312hd8f9ff3_0.conda + sha256: 2bec4e5d361ab39d0b06bd026484cdd6399307d4c0e74807c2e2a2cca2c859f2 + md5: a7fc6b79deae458d5522b5103002120c + depends: + - __osx >=11.0 + - libcxx >=18 + - ply + - python >=3.12,<3.13.0a0 + - python >=3.12,<3.13.0a0 *_cpython + - python_abi 3.12.* *_cp312 + - setuptools + arch: arm64 + platform: osx + license: BSD-3-Clause + license_family: BSD + size: 7421923 + timestamp: 1744833268609 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/python-3.12.10-hc22306f_0_cpython.conda + sha256: 69aed911271e3f698182e9a911250b05bdf691148b670a23e0bea020031e298e + md5: c88f1a7e1e7b917d9c139f03b0960fea + depends: + - __osx >=11.0 + - bzip2 >=1.0.8,<2.0a0 + - libexpat >=2.7.0,<3.0a0 + - libffi >=3.4.6,<3.5.0a0 + - liblzma >=5.8.1,<6.0a0 + - libsqlite >=3.49.1,<4.0a0 + - libzlib >=1.3.1,<2.0a0 + - ncurses >=6.5,<7.0a0 + - openssl >=3.5.0,<4.0a0 + - readline >=8.2,<9.0a0 + - tk >=8.6.13,<8.7.0a0 + - tzdata + constrains: + - python_abi 3.12.* *_cp312 + arch: arm64 + platform: osx + license: Python-2.0 + size: 12932743 + timestamp: 1744323815320 +- conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.12-7_cp312.conda + build_number: 7 + sha256: a1bbced35e0df66cc713105344263570e835625c28d1bdee8f748f482b2d7793 + md5: 0dfcdc155cf23812a0c9deada86fb723 + constrains: + - python 3.12.* *_cpython + license: BSD-3-Clause + license_family: BSD + size: 6971 + timestamp: 1745258861359 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/readline-8.2-h1d1bf99_2.conda + sha256: 7db04684d3904f6151eff8673270922d31da1eea7fa73254d01c437f49702e34 + md5: 63ef3f6e6d6d5c589e64f11263dc5676 + depends: + - ncurses >=6.5,<7.0a0 + arch: arm64 + platform: osx + license: GPL-3.0-only + license_family: GPL + size: 252359 + timestamp: 1740379663071 +- conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.1.0-pyhff2d567_0.conda + sha256: 777d34ed359cedd5a5004c930077c101bb3b70e5fbb04d29da5058d75b0ba487 + md5: f6f72d0837c79eaec77661be43e8a691 + depends: + - python >=3.9 + license: MIT + license_family: MIT + size: 778484 + timestamp: 1746085063737 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/tinyxml2-11.0.0-ha1acc90_0.conda + sha256: bbd9294551ff727305f8335819c24d2490d5d79e0f3d90957992c39d2146093a + md5: 6778d917f88222e8f27af8ec5c41f277 + depends: + - __osx >=11.0 + - libcxx >=18 + arch: arm64 + platform: osx + license: Zlib + size: 122269 + timestamp: 1742246179980 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/tk-8.6.13-h5083fa2_1.conda + sha256: 72457ad031b4c048e5891f3f6cb27a53cb479db68a52d965f796910e71a403a8 + md5: b50a57ba89c32b62428b71a875291c9b + depends: + - libzlib >=1.2.13,<2.0.0a0 + arch: arm64 + platform: osx + license: TCL + license_family: BSD + size: 3145523 + timestamp: 1699202432999 +- conda: https://conda.anaconda.org/conda-forge/noarch/tzdata-2025b-h78e105d_0.conda + sha256: 5aaa366385d716557e365f0a4e9c3fca43ba196872abbbe3d56bb610d131e192 + md5: 4222072737ccff51314b5ece9c7d6f5a + license: LicenseRef-Public-Domain + size: 122968 + timestamp: 1742727099393 diff --git a/.CondaPkg/pixi.toml b/.CondaPkg/pixi.toml new file mode 100644 index 0000000000..42e8ed6b17 --- /dev/null +++ b/.CondaPkg/pixi.toml @@ -0,0 +1,15 @@ +[dependencies] +pyomo = "*" +casadi = ">3.3" + + [dependencies.python] + channel = "conda-forge" + build = "*cpython*" + version = ">=3.8,<4" + +[project] +name = ".CondaPkg" +platforms = ["osx-arm64"] +channels = ["conda-forge"] +channel-priority = "strict" +description = "automatically generated by CondaPkg.jl" diff --git a/Project.toml b/Project.toml index 33280a5064..ff97c329a9 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" +CasADi = "c49709b8-5c63-11e9-2fb2-69db5844192f" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" @@ -64,7 +65,6 @@ 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" diff --git a/ext/MTKCasADiDynamicOptExt.jl b/ext/MTKCasADiDynamicOptExt.jl index 4fea4a94e0..6fd8fbdd65 100644 --- a/ext/MTKCasADiDynamicOptExt.jl +++ b/ext/MTKCasADiDynamicOptExt.jl @@ -5,13 +5,27 @@ using DiffEqBase using UnPack const MTK = ModelingToolkit +# 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::Opti + model::CasADiModel kwargs::K function CasADiDynamicOptProblem(f, u0, tspan, p, model, kwargs...) @@ -20,27 +34,13 @@ struct CasADiDynamicOptProblem{uType, tType, isinplace, P, F, K} <: end 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ₛ::Union{Number, MX} -end - function (M::MXLinearInterpolation)(τ) - nt = (τ - M.t) / M.dt + 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.") - M.u[i] + Δ*(M.u[i + 1] - M.u[i]) + M.u[:, i] + Δ*(M.u[:, i + 1] - M.u[:, i]) end """ @@ -66,11 +66,13 @@ function MTK.CasADiDynamicOptProblem(sys::ODESystem, u0map, tspan, pmap; 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, kwargs...) + 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) @@ -85,17 +87,28 @@ function init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t = false) tₛ = variable!(opti) tsteps = LinRange(0, 1, steps) else - tₛ = 1 + 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) - U_interp = MXLinearInterpolation(U, tsteps, tsteps[2]-tsteps[1]) V_interp = MXLinearInterpolation(V, tsteps, tsteps[2]-tsteps[1]) - CasADiModel(opti, U_interp, V_interp, tₛ) + model = CasADiModel(opti, U_interp, V_interp, tₛ) + + set_casadi_bounds!(model, sys, pmap) + add_cost_function!(model, sys, (tspan[1], tspan[2]), pmap) + add_user_constraints!(model, sys, 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) @@ -114,7 +127,7 @@ function set_casadi_bounds!(model, sys, pmap) end end -function add_initial_constraints!(model::CasADiModel, u0, u0_idxs, ts) +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]) @@ -124,19 +137,19 @@ end function add_user_constraints!(model::CasADiModel, sys, pmap; is_free_t = false) @unpack opti, U, V, tₛ = model - iv = get_iv(sys) + 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(sts)]) - pidxmap = Dict([v => i for (i, v) in enumerate(ps)]) + stidxmap = Dict([v => i for (i, v) in enumerate(unknowns(sys))]) cons_unknowns = map(MTK.default_toterm, unknowns(conssys)) for st in cons_unknowns - x = operation(st) - t = only(argments(st)) + x = MTK.operation(st) + t = only(MTK.arguments(st)) idx = stidxmap[x(iv)] - + @show t + MTK.symbolic_type(t) === MTK.NotSymbolic() || continue jconstraints = map(c -> Symbolics.substitute(c, Dict(x(t) => U(t)[idx])), jconstraints) end jconstraints = substitute_casadi_vars(model, sys, pmap, jconstraints) @@ -145,9 +158,9 @@ function add_user_constraints!(model::CasADiModel, sys, pmap; is_free_t = false) if cons isa Equation subject_to!(opti, cons.lhs - cons.rhs==0) elseif cons.relational_op === Symbolics.geq - subject_to!(model, cons.lhs - cons.rhs≥0) + subject_to!(opti, cons.lhs - cons.rhs≥0) else - subject_to!(model, cons.lhs - cons.rhs≤0) + subject_to!(opti, cons.lhs - cons.rhs≤0) end end end @@ -158,7 +171,7 @@ function add_cost_function!(model::CasADiModel, sys, tspan, pmap) consolidate = MTK.get_consolidate(sys) if isnothing(jcosts) || isempty(jcosts) - minimize!(opti, 0) + minimize!(opti, MX(0)) return end stidxmap = Dict([v => i for (i, v) in enumerate(sts)]) @@ -175,8 +188,6 @@ function add_cost_function!(model::CasADiModel, sys, tspan, pmap) end end jcosts = substitute_casadi_vars(model::CasADiModel, sys, pmap, jcosts; auxmap) - jcosts = map( - c -> Symbolics.substitute(c, MTK.∫() => Symbolics.Integral(iv in tspan)), jcosts) dt = U.t[2] - U.t[1] intmap = Dict() @@ -210,64 +221,70 @@ function substitute_casadi_vars(model::CasADiModel, sys, pmap, exprs; auxmap = D exprs end -function add_solve_constraints!(prob, tableau; is_free_t) +function add_solve_constraints(prob, tableau; is_free_t = false) @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ᵤ = length(U) - nᵥ = length(V) + nᵤ = size(U.u, 1) + nᵥ = size(V.u, 1) - if is_explicit(tableau) - K = Any[] + 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 = zeros(nᵤ)) + Δ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!(opti, U.u[:, k] + ΔU == U.u[:, k+1]) + subject_to!(solver_opti, U.u[:, k] + ΔU == U.u[:, k+1]) + empty!(K) end else - ΔU_tot = dt * (K' * α) for k in 1:length(tsteps)-1 - Kᵢ = variable!(opti, length(α), nᵤ) - ΔUs = A * Kᵢ # the stepsize at each stage of the implicit method + τ = 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 = @view ΔUs[i, :] - Uₙ = U.u[:,k] + ΔU + ΔU = ΔUs[i,:]' + Uₙ = U.u[:,k] + ΔU*dt Vₙ = V.u[:,k] - subject_to!(opti, K[i,:] == tₛ * f(Uₙ, Vₙ, p, τ + h*dt)) + subject_to!(solver_opti, Kᵢ[:,i] == tₛ * f(Uₙ, Vₙ, p, τ + h*dt)) end - ΔU_tot = dt*(Kᵢ'*α) - subject_to!(opti, U.u[:, k] + ΔU_tot == U.u[:,k+1]) + Δ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}, tableau_getter = constructDefault; plugin_options::Dict = Dict(), solver_options::Dict = Dict(), silent = false) - model = prob.model +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() - opti = model.opti + @unpack opti, U, V, tₛ = model - solver!(opti, solver, plugin_options, solver_options) - add_casadi_solve_constraints!(prob, tableau) - solver!(cmodel, "$solver", plugin_options, solver_options) + opti = add_solve_constraints(prob, tableau) + solver!(opti, "$solver", plugin_options, solver_options) failed = false + value_getter = nothing try - sol = solve(opti) + sol = CasADi.solve!(opti) value_getter = x -> CasADi.value(sol, x) catch ErrorException value_getter = x -> CasADi.debug_value(opti, x) @@ -275,14 +292,14 @@ function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, S end ts = value_getter(tₛ) * U.t - U_vals = value_getter(U) - U_vals = [[U_vals[i][j] for i in 1:length(U_vals)] for j in 1:length(ts)] + U_vals = value_getter(U.u) + U_vals = [[U_vals[i, j] for i in 1:size(U_vals, 1)] for j in 1:length(ts)] sol = DiffEqBase.build_solution(prob, tableau_getter, ts, U_vals) input_sol = nothing - if !isempty(V) - V_vals = value_getter(V) - V_vals = [[V_vals[i][j] for i in 1:length(V_vals)] for j in 1:length(ts)] + if prod(size(V.u)) != 0 + V_vals = value_getter(V.u) + 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 @@ -292,6 +309,6 @@ function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, S input_sol, SciMLBase.ReturnCode.ConvergenceFailure)) end - DynamicOptSolution(cmodel, sol, input_sol) + DynamicOptSolution(model, sol, input_sol) end end From 65975d4c03c115ebfd915ffdc491d408f042815c Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 5 May 2025 21:27:34 -0400 Subject: [PATCH 14/28] feat: all problems working --- ext/MTKCasADiDynamicOptExt.jl | 131 +++++++++++++++++++++++++--------- ext/MTKInfiniteOptExt.jl | 1 + 2 files changed, 97 insertions(+), 35 deletions(-) diff --git a/ext/MTKCasADiDynamicOptExt.jl b/ext/MTKCasADiDynamicOptExt.jl index 6fd8fbdd65..e8fb93968c 100644 --- a/ext/MTKCasADiDynamicOptExt.jl +++ b/ext/MTKCasADiDynamicOptExt.jl @@ -3,8 +3,16 @@ 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 @@ -40,7 +48,11 @@ function (M::MXLinearInterpolation)(τ) Δ = nt - i + 1 (i > length(M.t) || i < 1) && error("Cannot extrapolate past the tspan.") - M.u[:, i] + Δ*(M.u[:, i + 1] - M.u[:, i]) + if i < length(M.t) + M.u[:, i] + Δ*(M.u[:, i + 1] - M.u[:, i]) + else + M.u[:, i] + end end """ @@ -83,8 +95,16 @@ function init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t = false) 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.") + 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) @@ -93,14 +113,21 @@ function init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t = false) 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]) + 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[1], tspan[2]), pmap) - add_user_constraints!(model, sys, pmap; is_free_t) + 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]) @@ -116,13 +143,15 @@ function set_casadi_bounds!(model, sys, pmap) for (i, u) in enumerate(unknowns(sys)) if MTK.hasbounds(u) lo, hi = MTK.getbounds(u) - subject_to!(opti, lo <= U[i, :] <= hi) + 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, lo <= V[i, :] <= hi) + 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 @@ -134,7 +163,7 @@ function add_initial_constraints!(model::CasADiModel, u0, u0_idxs) end end -function add_user_constraints!(model::CasADiModel, sys, pmap; is_free_t = false) +function add_user_constraints!(model::CasADiModel, sys, tspan, pmap; is_free_t) @unpack opti, U, V, tₛ = model iv = MTK.get_iv(sys) @@ -143,18 +172,29 @@ function add_user_constraints!(model::CasADiModel, sys, pmap; is_free_t = false) (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)) - for st in cons_unknowns - x = MTK.operation(st) - t = only(MTK.arguments(st)) - idx = stidxmap[x(iv)] - @show t - MTK.symbolic_type(t) === MTK.NotSymbolic() || continue - jconstraints = map(c -> Symbolics.substitute(c, Dict(x(t) => U(t)[idx])), jconstraints) - end - jconstraints = substitute_casadi_vars(model, sys, pmap, jconstraints) + 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 @@ -165,29 +205,38 @@ function add_user_constraints!(model::CasADiModel, sys, pmap; is_free_t = false) end end -function add_cost_function!(model::CasADiModel, sys, tspan, pmap) +function add_cost_function!(model::CasADiModel, sys, tspan, pmap; is_free_t) @unpack opti, U, V, tₛ = model - jcosts = MTK.get_costs(sys) + jcosts = copy(MTK.get_costs(sys)) consolidate = MTK.get_consolidate(sys) - if isnothing(jcosts) || isempty(jcosts) minimize!(opti, MX(0)) return end - stidxmap = Dict([v => i for (i, v) in enumerate(sts)]) - pidxmap = Dict([v => i for (i, v) in enumerate(ps)]) + 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) - vars = vars(jcosts[i]) - for st in vars + costvars = MTK.vars(jcosts[i]) + for st in costvars + MTK.iscall(st) || continue x = operation(st) t = only(arguments(st)) - t isa Union{Num, MTK.Symbolic} && continue - idx = stidxmap[x(iv)] - jcosts[i] = Symbolics.substitute(jcosts[i], Dict(x(t) => U(t)[idx])) + 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 - jcosts = substitute_casadi_vars(model::CasADiModel, sys, pmap, jcosts; auxmap) dt = U.t[2] - U.t[1] intmap = Dict() @@ -195,15 +244,17 @@ function add_cost_function!(model::CasADiModel, sys, tspan, pmap) op = MTK.operation(int) arg = only(arguments(MTK.value(int))) lo, hi = (op.domain.domain.left, op.domain.domain.right) - (lo, hi) !== tspan && error("Non-whole interval bounds for integrals are not currently supported.") + !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) - minimize!(opti, consolidate(jcosts)) + jcosts = MTK.value.(jcosts) + minimize!(opti, MX(MTK.value(consolidate(jcosts)))) end -function substitute_casadi_vars(model::CasADiModel, sys, pmap, exprs; auxmap = Dict()) - @unpack opti, U, V = model +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) @@ -213,6 +264,13 @@ function substitute_casadi_vars(model::CasADiModel, sys, pmap, exprs; auxmap = D 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)]; @@ -221,7 +279,7 @@ function substitute_casadi_vars(model::CasADiModel, sys, pmap, exprs; auxmap = D exprs end -function add_solve_constraints(prob, tableau; is_free_t = false) +function add_solve_constraints(prob, tableau) @unpack A, α, c = tableau @unpack model, f, p = prob @unpack opti, U, V, tₛ = model @@ -283,6 +341,7 @@ function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, S failed = false value_getter = nothing + sol = nothing try sol = CasADi.solve!(opti) value_getter = x -> CasADi.value(sol, x) @@ -293,22 +352,24 @@ function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, S 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)] - sol = DiffEqBase.build_solution(prob, tableau_getter, ts, U_vals) + 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 - sol = SciMLBase.solution_new_retcode(sol, SciMLBase.ReturnCode.ConvergenceFailure) + 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, sol, input_sol) + DynamicOptSolution(model, ode_sol, input_sol) end end diff --git a/ext/MTKInfiniteOptExt.jl b/ext/MTKInfiniteOptExt.jl index 5a48f9b96f..27361bd0c1 100644 --- a/ext/MTKInfiniteOptExt.jl +++ b/ext/MTKInfiniteOptExt.jl @@ -173,6 +173,7 @@ function add_jump_cost_function!(model::InfiniteModel, sys, tspan, pmap; is_free return end jcosts = substitute_jump_vars(model, sys, pmap, jcosts; is_free_t) + @show jcosts tₛ = is_free_t ? model[:tf] : 1 # Substitute integral From 583d295a5a6e61cde2f7a9b27ef2a7644ba60af7 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 5 May 2025 21:28:14 -0400 Subject: [PATCH 15/28] rename dynamic optimization test file --- test/downstream/dynamic_opt_systems.jl | 37 --- test/downstream/dynamic_optimization.jl | 353 ++++++++++++++++++++++++ 2 files changed, 353 insertions(+), 37 deletions(-) delete mode 100644 test/downstream/dynamic_opt_systems.jl create mode 100644 test/downstream/dynamic_optimization.jl diff --git a/test/downstream/dynamic_opt_systems.jl b/test/downstream/dynamic_opt_systems.jl deleted file mode 100644 index f13b6446ef..0000000000 --- a/test/downstream/dynamic_opt_systems.jl +++ /dev/null @@ -1,37 +0,0 @@ -function build_lotkavolterra(; with_constraint = false) - @parameters α=1.5 β=1.0 γ=3.0 δ=1.0 - @variables x(..) y(..) - t = M.t_nounits - D = M.D_nounits - - eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y(t), - D(y(t)) ~ -γ * y(t) + δ * x(t) * y(t)] - - tspan = (0.0, 1.0) - parammap = [α => 1.5, β => 1.0, γ => 3.0, δ => 1.0] - - if with_constraint - constr = [x(0.6) ~ 3.5, x(0.3) ~ 7.0] - guess = [x(t) => 4.0, y(t) => 2.0] - u0map = Pair[] - else - constr = nothing - guess = Pair[] - u0map = [x(t) => 4.0, y(t) => 2.0] - end - - @mtkbuild sys = ODESystem(eqs, t; constraints = constr) - sys, u0map, tspan, parammap, guess -end - -function rocket_fft() - -end - -function rocket() - -end - -function cartpole() - -end diff --git a/test/downstream/dynamic_optimization.jl b/test/downstream/dynamic_optimization.jl new file mode 100644 index 0000000000..dfe3e6a36d --- /dev/null +++ b/test/downstream/dynamic_optimization.jl @@ -0,0 +1,353 @@ +using ModelingToolkit +import JuMP, InfiniteOpt +using DiffEqDevTools, DiffEqBase +using SimpleDiffEq +using OrdinaryDiffEqSDIRK, OrdinaryDiffEqVerner, OrdinaryDiffEqTsit5, OrdinaryDiffEqFIRK +using Ipopt +using DataInterpolations +using CasADi + +const M = ModelingToolkit + +probtypes = [JuMPDynamicOptProblem, InfiniteOptDynamicOptProblem, CasADiDynamicOptProblem] +solvers = [Ipopt.Optimizer, Ipopt.Optimizer, "opti"] + +function test_dynamic_opt_probs(tableaus, true_sol, args...; kwargs...) + for (probtype, solver, tableau) in zip(probtypes, solvers tableaus) + prob = probtype(args...; kwargs...) + sol = solve(prob, solver, tableau, silent = true) + end +end + +@testset "ODE Solution, no cost" begin + # Test solving without anything attached. + @parameters α=1.5 β=1.0 γ=3.0 δ=1.0 + @variables x(..) y(..) + t = M.t_nounits + D = M.D_nounits + + eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y(t), + D(y(t)) ~ -γ * y(t) + δ * x(t) * y(t)] + + @mtkbuild sys = ODESystem(eqs, t) + tspan = (0.0, 1.0) + u0map = [x(t) => 4.0, y(t) => 2.0] + parammap = [α => 1.5, β => 1.0, γ => 3.0, δ => 1.0] + + # Test explicit method. + jprob = JuMPDynamicOptProblem(sys, u0map, tspan, parammap, dt = 0.01) + @test JuMP.num_constraints(jprob.model) == 2 # initials + 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(oprob, "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 ≈(isol2.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[] + guess = [x(t) => 4.0, y(t) => 2.0] + constr = [x(0.6) ~ 3.5, x(0.3) ~ 7.0] + @mtkbuild lksys = ODESystem(eqs, t; constraints = constr) + + jprob = JuMPDynamicOptProblem(lksys, u0map, tspan, parammap; guesses = guess, dt = 0.01) + @test JuMP.num_constraints(jprob.model) == 2 + jsol = solve(jprob, Ipopt.Optimizer, constructTsitouras5, silent = true) # 12.190 s, 9.68 GiB + @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, + derivative_method = InfiniteOpt.OrthogonalCollocation(3), silent = true) # 48.564 ms, 9.58 MiB + sol = isol.sol + @test sol(0.6)[1] ≈ 3.5 + @test sol(0.3)[1] ≈ 7.0 + + # Test whole-interval constraints + constr = [x(t) ≳ 1, y(t) ≳ 1] + @mtkbuild lksys = ODESystem(eqs, t; constraints = constr) + iprob = InfiniteOptDynamicOptProblem( + lksys, u0map, tspan, parammap; guesses = guess, dt = 0.01) + isol = solve(iprob, Ipopt.Optimizer, + 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) + for v in 1:(length(input_sol.u[1]) - 1) + all(i -> ≈(i[v], bounds[v]; rtol) || ≈(i[v], bounds[u]; rtol), input_sol.u) || + return false + end + true +end + +function ctrl_to_spline(inputsol, splineType) + us = reduce(vcat, inputsol.u) + ts = reduce(vcat, inputsol.t) + splineType(us, ts) +end + +@testset "Linear systems" begin + # Double integrator + t = M.t_nounits + D = M.D_nounits + @variables x(..) [bounds = (0.0, 0.25)] v(..) + @variables u(..) [bounds = (-1.0, 1.0), input = true] + constr = [v(1.0) ~ 0.0] + cost = [-x(1.0)] # Maximize the final distance. + @named block = ODESystem( + [D(x(t)) ~ v(t), D(v(t)) ~ u(t)], t; costs = cost, constraints = constr) + block, input_idxs = structural_simplify(block, ([u(t)], [])) + + u0map = [x(t) => 0.0, v(t) => 0.0] + tspan = (0.0, 1.0) + parammap = [u(t) => 0.0] + jprob = JuMPDynamicOptProblem(block, u0map, tspan, parammap; dt = 0.01) + jsol = solve(jprob, Ipopt.Optimizer, constructVerner8, silent = true) + # Linear systems have bang-bang controls + @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) + spline = ctrl_to_spline(jsol.input_sol, ConstantInterpolation) + 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) + @test is_bangbang(isol.input_sol, [-1.0], [1.0]) + @test ≈(isol.sol.u[end][2], 0.25, rtol = 1e-5) + osol = solve(oprob, ImplicitEuler(); dt = 0.01, adaptive = false) + @test ≈(isol.sol.u, osol.u, rtol = 0.05) + + ################### + ### Bee example ### + ################### + # From Lawrence Evans' notes + @variables w(..) q(..) α(t) [input = true, bounds = (0, 1)] + @parameters b c μ s ν + + tspan = (0, 4) + eqs = [D(w(t)) ~ -μ * w(t) + b * s * α * w(t), + D(q(t)) ~ -ν * q(t) + c * (1 - α) * s * w(t)] + costs = [-q(tspan[2])] + + @named beesys = ODESystem(eqs, t; costs) + beesys, input_idxs = structural_simplify(beesys, ([α], [])) + u0map = [w(t) => 40, q(t) => 2] + pmap = [b => 1, c => 1, μ => 1, s => 1, ν => 1, α => 1] + + jprob = JuMPDynamicOptProblem(beesys, u0map, tspan, pmap, dt = 0.01) + jsol = solve(jprob, Ipopt.Optimizer, constructTsitouras5, silent = true) + @test is_bangbang(jsol.input_sol, [0.0], [1.0]) + 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), + D(q(t)) ~ -ν * q(t) + c * (1 - α_interp(t)) * s * w(t)] + @mtkbuild beesys_ode = ODESystem(eqs, t) + oprob = ODEProblem(beesys_ode, + u0map, + tspan, + merge(Dict(pmap), + 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 + +@testset "Rocket launch" begin + t = M.t_nounits + D = M.D_nounits + + @parameters h_c m₀ h₀ g₀ D_c c Tₘ m_c + @variables h(..) v(..) m(..) [bounds = (m_c, 1)] T(..) [input = true, bounds = (0, Tₘ)] + drag(h, v) = D_c * v^2 * exp(-h_c * (h - h₀) / h₀) + gravity(h) = g₀ * (h₀ / h) + + eqs = [D(h(t)) ~ v(t), + D(v(t)) ~ (T(t) - drag(h(t), v(t))) / m(t) - gravity(h(t)), + D(m(t)) ~ -T(t) / c] + + (ts, te) = (0.0, 0.2) + costs = [-h(te)] + cons = [T(te) ~ 0, m(te) ~ m_c] + @named rocket = ODESystem(eqs, t; costs, constraints = cons) + rocket, input_idxs = structural_simplify(rocket, ([T(t)], [])) + + u0map = [h(t) => h₀, m(t) => m₀, v(t) => 0] + pmap = [ + g₀ => 1, m₀ => 1.0, h_c => 500, c => 0.5 * √(g₀ * h₀), D_c => 0.5 * 620 * m₀ / g₀, + Tₘ => 3.5 * g₀ * m₀, T(t) => 0.0, h₀ => 1, m_c => 0.6] + 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) + isol = solve(iprob, Ipopt.Optimizer, silent = true) + @test isol.sol.u[end][1] > 1.012 + + # Test solution + @parameters (T_interp::CubicSpline)(..) + eqs = [D(h(t)) ~ v(t), + D(v(t)) ~ (T_interp(t) - drag(h(t), v(t))) / m(t) - gravity(h(t)), + D(m(t)) ~ -T_interp(t) / c] + @mtkbuild rocket_ode = ODESystem(eqs, t) + interpmap = Dict(T_interp => ctrl_to_spline(jsol.input_sol, CubicSpline)) + 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)) + osol1 = solve(oprob1, ImplicitEuler(); adaptive = false, dt = 0.001) + @test ≈(isol.sol.u, osol1.u, rtol = 0.01) +end + +@testset "Free final time problems" begin + t = M.t_nounits + D = M.D_nounits + + @variables x(..) u(..) [input = true, bounds = (0, 1)] + @parameters tf + eqs = [D(x(t)) ~ -2 + 0.5 * u(t)] + # Integral cost function + costs = [-Symbolics.Integral(t in (0, tf))(x(t) - u(t)), -x(tf)] + consolidate(u) = u[1] + u[2] + @named rocket = ODESystem(eqs, t; costs, consolidate) + rocket, input_idxs = structural_simplify(rocket, ([u(t)], [])) + + u0map = [x(t) => 17.5] + pmap = [u(t) => 0.0, tf => 8] + jprob = JuMPDynamicOptProblem(rocket, u0map, (0, tf), pmap; steps = 201) + 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) + + @variables x(..) v(..) + @variables u(..) [bounds = (-1.0, 1.0), input = true] + @parameters tf + constr = [v(tf) ~ 0, x(tf) ~ 0] + cost = [tf] # Minimize time + + @named block = ODESystem( + [D(x(t)) ~ v(t), D(v(t)) ~ u(t)], t; costs = cost, constraints = constr) + block, input_idxs = structural_simplify(block, ([u(t)], [])) + + u0map = [x(t) => 1.0, v(t) => 0.0] + tspan = (0.0, tf) + parammap = [u(t) => 0.0, tf => 1.0] + jprob = JuMPDynamicOptProblem(block, u0map, tspan, parammap; steps = 51) + 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) +end + +@testset "Cart-pole problem" begin + t = M.t_nounits + D = M.D_nounits + # gravity, length, moment of Inertia, drag coeff + @parameters g l mₚ mₖ + @variables x(..) θ(..) u(t) [input = true, bounds = (-10, 10)] + + s = sin(θ(t)) + c = cos(θ(t)) + H = [mₖ+mₚ mₚ*l*c + mₚ*l*c mₚ*l^2] + C = [0 -mₚ*D(θ(t))*l*s + 0 0] + qd = [D(x(t)), D(θ(t))] + G = [0, mₚ * g * l * s] + B = [1, 0] + + tf = 5 + rhss = -H \ Vector(C * qd + G - B * u) + eqs = [D(D(x(t))) ~ rhss[1], D(D(θ(t))) ~ rhss[2]] + cons = [θ(tf) ~ π, x(tf) ~ 0, D(θ(tf)) ~ 0, D(x(tf)) ~ 0] + costs = [Symbolics.Integral(t in (0, tf))(u^2)] + tspan = (0, tf) + + @named cartpole = ODESystem(eqs, t; costs, constraints = cons) + cartpole, input_idxs = structural_simplify(cartpole, ([u], [])) + + u0map = [D(x(t)) => 0.0, D(θ(t)) => 0.0, θ(t) => 0.0, x(t) => 0.0] + pmap = [mₖ => 1.0, mₚ => 0.2, l => 0.5, g => 9.81, u => 0] + jprob = JuMPDynamicOptProblem(cartpole, u0map, tspan, pmap; dt = 0.04) + 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] +end From 3b5ccfff8b86215d2246702cf90e205765ee2e2d Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 5 May 2025 21:45:49 -0400 Subject: [PATCH 16/28] check tests working --- ext/MTKCasADiDynamicOptExt.jl | 3 ++- ext/MTKInfiniteOptExt.jl | 3 +-- test/downstream/dynamic_optimization.jl | 14 ++------------ 3 files changed, 5 insertions(+), 15 deletions(-) diff --git a/ext/MTKCasADiDynamicOptExt.jl b/ext/MTKCasADiDynamicOptExt.jl index e8fb93968c..2c4110b080 100644 --- a/ext/MTKCasADiDynamicOptExt.jl +++ b/ext/MTKCasADiDynamicOptExt.jl @@ -115,7 +115,7 @@ function init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t = false) 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]) - set_initial!(opti, V, DM(repeat(c0, 1, steps))) + !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]) @@ -337,6 +337,7 @@ function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, S @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 diff --git a/ext/MTKInfiniteOptExt.jl b/ext/MTKInfiniteOptExt.jl index 27361bd0c1..7c374cfdaf 100644 --- a/ext/MTKInfiniteOptExt.jl +++ b/ext/MTKInfiniteOptExt.jl @@ -173,7 +173,6 @@ function add_jump_cost_function!(model::InfiniteModel, sys, tspan, pmap; is_free return end jcosts = substitute_jump_vars(model, sys, pmap, jcosts; is_free_t) - @show jcosts tₛ = is_free_t ? model[:tf] : 1 # Substitute integral @@ -338,7 +337,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/test/downstream/dynamic_optimization.jl b/test/downstream/dynamic_optimization.jl index dfe3e6a36d..abb534ea0b 100644 --- a/test/downstream/dynamic_optimization.jl +++ b/test/downstream/dynamic_optimization.jl @@ -9,16 +9,6 @@ using CasADi const M = ModelingToolkit -probtypes = [JuMPDynamicOptProblem, InfiniteOptDynamicOptProblem, CasADiDynamicOptProblem] -solvers = [Ipopt.Optimizer, Ipopt.Optimizer, "opti"] - -function test_dynamic_opt_probs(tableaus, true_sol, args...; kwargs...) - for (probtype, solver, tableau) in zip(probtypes, solvers tableaus) - prob = probtype(args...; kwargs...) - sol = solve(prob, solver, tableau, silent = true) - end -end - @testset "ODE Solution, no cost" begin # Test solving without anything attached. @parameters α=1.5 β=1.0 γ=3.0 δ=1.0 @@ -41,7 +31,7 @@ end oprob = ODEProblem(sys, u0map, tspan, parammap) osol = solve(oprob, SimpleRK4(), dt = 0.01) cprob = CasADiDynamicOptProblem(sys, u0map, tspan, parammap, dt = 0.01) - csol = solve(oprob, "ipopt", constructRK4) + csol = solve(cprob, "ipopt", constructRK4) @test jsol.sol.u ≈ osol.u @test csol.sol.u ≈ osol.u @@ -54,7 +44,7 @@ end isol = solve(iprob, Ipopt.Optimizer, derivative_method = InfiniteOpt.FiniteDifference(InfiniteOpt.Backward()), silent = true) # 11.540 ms, 4.00 MiB - @test ≈(isol2.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) From aa79f37eeb51b375e96a10529e589ec99d5ae5c9 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 5 May 2025 22:05:45 -0400 Subject: [PATCH 17/28] add silen t option --- ext/MTKCasADiDynamicOptExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/MTKCasADiDynamicOptExt.jl b/ext/MTKCasADiDynamicOptExt.jl index 2c4110b080..9d16a3ea5c 100644 --- a/ext/MTKCasADiDynamicOptExt.jl +++ b/ext/MTKCasADiDynamicOptExt.jl @@ -337,7 +337,7 @@ function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, S @unpack opti, U, V, tₛ = model opti = add_solve_constraints(prob, tableau) - silent && solver_options["print_level"] = 0 + silent && (solver_options["print_level"] = 0) solver!(opti, "$solver", plugin_options, solver_options) failed = false From c5d5853df1c83885eabbe090b930fd3a085a0c4d Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 5 May 2025 22:08:11 -0400 Subject: [PATCH 18/28] update gitignore --- .CondaPkg/.gitattributes | 2 - .CondaPkg/.gitignore | 4 - .CondaPkg/meta | Bin 624 -> 0 bytes .CondaPkg/pixi.lock | 540 --------------------------------------- .CondaPkg/pixi.toml | 15 -- .gitignore | 1 + 6 files changed, 1 insertion(+), 561 deletions(-) delete mode 100644 .CondaPkg/.gitattributes delete mode 100644 .CondaPkg/.gitignore delete mode 100644 .CondaPkg/meta delete mode 100644 .CondaPkg/pixi.lock delete mode 100644 .CondaPkg/pixi.toml diff --git a/.CondaPkg/.gitattributes b/.CondaPkg/.gitattributes deleted file mode 100644 index 887a2c18f0..0000000000 --- a/.CondaPkg/.gitattributes +++ /dev/null @@ -1,2 +0,0 @@ -# SCM syntax highlighting & preventing 3-way merges -pixi.lock merge=binary linguist-language=YAML linguist-generated=true diff --git a/.CondaPkg/.gitignore b/.CondaPkg/.gitignore deleted file mode 100644 index 740bb7d1ae..0000000000 --- a/.CondaPkg/.gitignore +++ /dev/null @@ -1,4 +0,0 @@ - -# pixi environments -.pixi -*.egg-info diff --git a/.CondaPkg/meta b/.CondaPkg/meta deleted file mode 100644 index f82e7e2daa6bea42c0dfd7469940ff53a33ba22d..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 624 zcmb7BOHRWu5DmWygv1TFKxsWqQcxFAsR|nwm9hX;T@JCEx^ZpGb_!jv=Zsu|OOQ4W zLWnHz#`Ae^-uU?UwL>p(*)Fs*EI~WD=e>52 z#;m}cSxC2Tsbqpez-&634?3PRAQ6d1$39Cdmm6anM1~eAZxG{{G>&^t5S;i(Z`E3T tSAPY~IK5xw)OW{sF&Xu4hvz=eb|2nfD3h}@U+QKxrF-BDe_(wl_yTx^l-mFR diff --git a/.CondaPkg/pixi.lock b/.CondaPkg/pixi.lock deleted file mode 100644 index 43d038b1dd..0000000000 --- a/.CondaPkg/pixi.lock +++ /dev/null @@ -1,540 +0,0 @@ -version: 6 -environments: - default: - channels: - - url: https://conda.anaconda.org/conda-forge/ - packages: - osx-arm64: - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/ampl-asl-1.0.0-h286801f_2.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/bzip2-1.0.8-h99b78c6_7.conda - - conda: https://conda.anaconda.org/conda-forge/noarch/ca-certificates-2025.4.26-hbd8a1cb_0.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/casadi-3.7.0-py312h74b190f_2.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/ipopt-3.14.17-h945cc1c_2.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libblas-3.9.0-31_h10e41b3_openblas.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libblasfeo-0.1.4.2-h37ef02a_0.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libcblas-3.9.0-31_hb3479ef_openblas.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libcxx-20.1.4-ha82da77_0.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libexpat-2.7.0-h286801f_0.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libfatrop-0.0.4-h286801f_1.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libffi-3.4.6-h1da3d7d_1.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libgfortran-5.0.0-14_2_0_h6c33f7e_103.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libgfortran5-14.2.0-h6c33f7e_103.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/liblapack-3.9.0-31_hc9a63f6_openblas.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/liblzma-5.8.1-h39f12f2_0.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libopenblas-0.3.29-openmp_hf332438_0.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libosqp-0.6.3-h5833ebf_1.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libqdldl-0.1.7-hb7217d7_0.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libscotch-7.0.6-he56f69b_1.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libsqlite-3.49.1-h3f77e49_2.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libzlib-1.3.1-h8359307_2.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/llvm-openmp-20.1.4-hdb05f8b_0.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/metis-5.1.0-h15f6cfe_1007.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/mumps-include-5.7.3-h8c5b6c6_10.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/mumps-seq-5.7.3-h390d176_10.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/ncurses-6.5-h5e97a16_3.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numpy-2.2.5-py312h7c1f314_0.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/openssl-3.5.0-h81ee809_1.conda - - conda: https://conda.anaconda.org/conda-forge/noarch/ply-3.11-pyhd8ed1ab_3.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/pyomo-6.9.2-py312hd8f9ff3_0.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/python-3.12.10-hc22306f_0_cpython.conda - - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.12-7_cp312.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/readline-8.2-h1d1bf99_2.conda - - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.1.0-pyhff2d567_0.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/tinyxml2-11.0.0-ha1acc90_0.conda - - conda: https://conda.anaconda.org/conda-forge/osx-arm64/tk-8.6.13-h5083fa2_1.conda - - conda: https://conda.anaconda.org/conda-forge/noarch/tzdata-2025b-h78e105d_0.conda -packages: -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/ampl-asl-1.0.0-h286801f_2.conda - sha256: 47a5da6cd38a32c086017a5d2ed2b67e0cc24e25268a9c4cc91ea7d8f1cb9507 - md5: bb25b8fa2b28474412dda4e1a95853b4 - depends: - - __osx >=11.0 - - libcxx >=18 - constrains: - - ampl-mp >=4.0.0 - arch: arm64 - platform: osx - license: BSD-3-Clause AND SMLNJ - size: 411989 - timestamp: 1732439500161 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/bzip2-1.0.8-h99b78c6_7.conda - sha256: adfa71f158cbd872a36394c56c3568e6034aa55c623634b37a4836bd036e6b91 - md5: fc6948412dbbbe9a4c9ddbbcfe0a79ab - depends: - - __osx >=11.0 - arch: arm64 - platform: osx - license: bzip2-1.0.6 - license_family: BSD - size: 122909 - timestamp: 1720974522888 -- conda: https://conda.anaconda.org/conda-forge/noarch/ca-certificates-2025.4.26-hbd8a1cb_0.conda - sha256: 2a70ed95ace8a3f8a29e6cd1476a943df294a7111dfb3e152e3478c4c889b7ac - md5: 95db94f75ba080a22eb623590993167b - depends: - - __unix - license: ISC - size: 152283 - timestamp: 1745653616541 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/casadi-3.7.0-py312h74b190f_2.conda - sha256: ae008bb0cd026b64d6cbcd7726fbbbb98904bd868fca3b7d53f258f0c80c27e7 - md5: 8431af5d6849a519573d8ed6ffaa8e28 - depends: - - __osx >=11.0 - - ipopt >=3.14.17,<3.14.18.0a0 - - libblas >=3.9.0,<4.0a0 - - libblasfeo >=0.1.4.2,<0.1.5.0a0 - - libcblas >=3.9.0,<4.0a0 - - libcxx >=18 - - libfatrop >=0.0.4,<0.0.5.0a0 - - liblapack >=3.9.0,<4.0a0 - - libosqp >=0.6.3,<0.6.4.0a0 - - numpy >=1.19,<3 - - python >=3.12,<3.13.0a0 - - python >=3.12,<3.13.0a0 *_cpython - - python_abi 3.12.* *_cp312 - - tinyxml2 >=11.0.0,<11.1.0a0 - arch: arm64 - platform: osx - license: LGPL-3.0-or-later - license_family: LGPL - size: 4681272 - timestamp: 1745474347317 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/ipopt-3.14.17-h945cc1c_2.conda - sha256: 379d9b44fc265493a8a42ca3694bb8d7538ed8c7bbaf8626dcf8d75f454317cc - md5: e465f880544e1e6b8d904fb65e25b130 - depends: - - __osx >=11.0 - - ampl-asl >=1.0.0,<1.0.1.0a0 - - libblas >=3.9.0,<4.0a0 - - libcxx >=18 - - liblapack >=3.9.0,<4.0a0 - - mumps-seq >=5.7.3,<5.7.4.0a0 - arch: arm64 - platform: osx - license: EPL-1.0 - size: 727718 - timestamp: 1745419955303 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libblas-3.9.0-31_h10e41b3_openblas.conda - build_number: 31 - sha256: 369586e7688b59b4f92c709b99d847d66d4d095425db327dd32ee5e6ab74697f - md5: 39b053da5e7035c6592102280aa7612a - depends: - - libopenblas >=0.3.29,<0.3.30.0a0 - - libopenblas >=0.3.29,<1.0a0 - constrains: - - liblapacke =3.9.0=31*_openblas - - libcblas =3.9.0=31*_openblas - - blas =2.131=openblas - - mkl <2025 - - liblapack =3.9.0=31*_openblas - arch: arm64 - platform: osx - license: BSD-3-Clause - license_family: BSD - size: 17123 - timestamp: 1740088119350 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libblasfeo-0.1.4.2-h37ef02a_0.conda - sha256: 406885b43c417025be9b2cdb541fa5a0f5426efbbf109d42ac5126efb5d3608d - md5: db8bf9db5c73a85b195f5c973f5514bf - depends: - - __osx >=11.0 - arch: arm64 - platform: osx - license: BSD-2-Clause - license_family: BSD - size: 509787 - timestamp: 1740067975938 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libcblas-3.9.0-31_hb3479ef_openblas.conda - build_number: 31 - sha256: f237486cc9118d09d0f3ff8820280de34365f98ee7b7dc5ab923b04c7cbf25a5 - md5: 7353c2bf0e90834cb70545671996d871 - depends: - - libblas 3.9.0 31_h10e41b3_openblas - constrains: - - liblapacke =3.9.0=31*_openblas - - blas =2.131=openblas - - liblapack =3.9.0=31*_openblas - arch: arm64 - platform: osx - license: BSD-3-Clause - license_family: BSD - size: 17032 - timestamp: 1740088127097 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libcxx-20.1.4-ha82da77_0.conda - sha256: 1837e2c65f8fc8cfd8b240cfe89406d0ce83112ac63f98c9fb3c9a15b4f2d4e1 - md5: 10c809af502fcdab799082d338170994 - depends: - - __osx >=11.0 - arch: arm64 - platform: osx - license: Apache-2.0 WITH LLVM-exception - license_family: Apache - size: 565811 - timestamp: 1745991653948 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libexpat-2.7.0-h286801f_0.conda - sha256: ee550e44765a7bbcb2a0216c063dcd53ac914a7be5386dd0554bd06e6be61840 - md5: 6934bbb74380e045741eb8637641a65b - depends: - - __osx >=11.0 - constrains: - - expat 2.7.0.* - arch: arm64 - platform: osx - license: MIT - license_family: MIT - size: 65714 - timestamp: 1743431789879 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libfatrop-0.0.4-h286801f_1.conda - sha256: ff193f966ee33546a297abe02e49b57ba2b728d4706cec9cf3898dc8ca81e3ee - md5: 00cb2b510ea71f73263e59e8253d0720 - depends: - - __osx >=11.0 - - libblasfeo >=0.1.4.1,<0.1.5.0a0 - - libcxx >=18 - arch: arm64 - platform: osx - license: LGPL-3.0-or-later - license_family: LGPL - size: 195859 - timestamp: 1736273349013 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libffi-3.4.6-h1da3d7d_1.conda - sha256: c6a530924a9b14e193ea9adfe92843de2a806d1b7dbfd341546ece9653129e60 - md5: c215a60c2935b517dcda8cad4705734d - depends: - - __osx >=11.0 - arch: arm64 - platform: osx - license: MIT - license_family: MIT - size: 39839 - timestamp: 1743434670405 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libgfortran-5.0.0-14_2_0_h6c33f7e_103.conda - sha256: 8628746a8ecd311f1c0d14bb4f527c18686251538f7164982ccbe3b772de58b5 - md5: 044a210bc1d5b8367857755665157413 - depends: - - libgfortran5 14.2.0 h6c33f7e_103 - arch: arm64 - platform: osx - license: GPL-3.0-only WITH GCC-exception-3.1 - license_family: GPL - size: 156291 - timestamp: 1743863532821 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libgfortran5-14.2.0-h6c33f7e_103.conda - sha256: 8599453990bd3a449013f5fa3d72302f1c68f0680622d419c3f751ff49f01f17 - md5: 69806c1e957069f1d515830dcc9f6cbb - depends: - - llvm-openmp >=8.0.0 - constrains: - - libgfortran 5.0.0 14_2_0_*_103 - arch: arm64 - platform: osx - license: GPL-3.0-only WITH GCC-exception-3.1 - license_family: GPL - size: 806566 - timestamp: 1743863491726 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/liblapack-3.9.0-31_hc9a63f6_openblas.conda - build_number: 31 - sha256: fe55b9aaf82c6c0192c3d1fcc9b8e884f97492dda9a8de5dae29334b3135fab5 - md5: ff57a55a2cbce171ef5707fb463caf19 - depends: - - libblas 3.9.0 31_h10e41b3_openblas - constrains: - - liblapacke =3.9.0=31*_openblas - - libcblas =3.9.0=31*_openblas - - blas =2.131=openblas - arch: arm64 - platform: osx - license: BSD-3-Clause - license_family: BSD - size: 17033 - timestamp: 1740088134988 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/liblzma-5.8.1-h39f12f2_0.conda - sha256: 4291dde55ebe9868491dc29716b84ac3de21b8084cbd4d05c9eea79d206b8ab7 - md5: ba24e6f25225fea3d5b6912e2ac562f8 - depends: - - __osx >=11.0 - arch: arm64 - platform: osx - license: 0BSD - size: 92295 - timestamp: 1743771392206 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libopenblas-0.3.29-openmp_hf332438_0.conda - sha256: 8989d9e01ec8c9b2d48dbb5efbe70b356fcd15990fb53b64fcb84798982c0343 - md5: 0cd1148c68f09027ee0b0f0179f77c30 - depends: - - __osx >=11.0 - - libgfortran >=5 - - libgfortran5 >=13.2.0 - - llvm-openmp >=18.1.8 - constrains: - - openblas >=0.3.29,<0.3.30.0a0 - arch: arm64 - platform: osx - license: BSD-3-Clause - license_family: BSD - size: 4168442 - timestamp: 1739825514918 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libosqp-0.6.3-h5833ebf_1.conda - sha256: 05740cc58df7a9c966a277331dea1987db485508c0a35d2ac9ef562f8061df06 - md5: 7b4c6d7f35e9fe945e810e213d21b942 - depends: - - __osx >=11.0 - - libcxx >=17 - - libqdldl >=0.1.7,<0.1.8.0a0 - arch: arm64 - platform: osx - license: Apache-2.0 - license_family: APACHE - size: 360361 - timestamp: 1729342406705 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libqdldl-0.1.7-hb7217d7_0.conda - sha256: 39bb718bca8ce87afdc592d857944ba39bfac54c31a1de4d669597b52bc668ea - md5: b4df6812b4ab33c1a520287f46d91220 - depends: - - libcxx >=14.0.6 - arch: arm64 - platform: osx - license: Apache-2.0 - license_family: APACHE - size: 17189 - timestamp: 1680741853527 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libscotch-7.0.6-he56f69b_1.conda - sha256: b616a5943501e35f32699b5f29d4216a3a94bfca87b9ae30feda89a9fa5c2f7a - md5: 291ffb6a18b7e35bf42a317cc6856348 - depends: - - __osx >=11.0 - - bzip2 >=1.0.8,<2.0a0 - - libgfortran >=5 - - libgfortran5 >=13.2.0 - - liblzma >=5.6.3,<6.0a0 - - libzlib >=1.3.1,<2.0a0 - arch: arm64 - platform: osx - license: CECILL-C - size: 273604 - timestamp: 1737537249207 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libsqlite-3.49.1-h3f77e49_2.conda - sha256: 907a95f73623c343fc14785cbfefcb7a6b4f2bcf9294fcb295c121611c3a590d - md5: 3b1e330d775170ac46dff9a94c253bd0 - depends: - - __osx >=11.0 - - libzlib >=1.3.1,<2.0a0 - arch: arm64 - platform: osx - license: Unlicense - size: 900188 - timestamp: 1742083865246 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libzlib-1.3.1-h8359307_2.conda - sha256: ce34669eadaba351cd54910743e6a2261b67009624dbc7daeeafdef93616711b - md5: 369964e85dc26bfe78f41399b366c435 - depends: - - __osx >=11.0 - constrains: - - zlib 1.3.1 *_2 - arch: arm64 - platform: osx - license: Zlib - license_family: Other - size: 46438 - timestamp: 1727963202283 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/llvm-openmp-20.1.4-hdb05f8b_0.conda - sha256: b8e8547116dba85890d7b39bfad1c86ed69a6b923caed1e449c90850d271d4d5 - md5: 00cbae3f2127efef6db76bd423a09807 - depends: - - __osx >=11.0 - constrains: - - openmp 20.1.4|20.1.4.* - arch: arm64 - platform: osx - license: Apache-2.0 WITH LLVM-exception - size: 282599 - timestamp: 1746134861758 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/metis-5.1.0-h15f6cfe_1007.conda - sha256: f54ad3e5d47a0235ba2830848fee590faad550639336fe1e2413ab16fee7ac39 - md5: 7687ec5796288536947bf616179726d8 - depends: - - __osx >=11.0 - arch: arm64 - platform: osx - license: Apache-2.0 - license_family: APACHE - size: 3898314 - timestamp: 1728064659078 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/mumps-include-5.7.3-h8c5b6c6_10.conda - sha256: 9093f47e4a15fd34c90de93ebfa30c63799e2f99123420fd68653274f5574131 - md5: 077b856b81445e47942980465ca52450 - arch: arm64 - platform: osx - license: CECILL-C - size: 20799 - timestamp: 1745406952431 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/mumps-seq-5.7.3-h390d176_10.conda - sha256: 2afb273e8268a0960eb3934edffecfa617c941cccc762a11c86d15cb4f0e17a8 - md5: d211ba253191e04aa341a02203ed42fb - depends: - - mumps-include ==5.7.3 h8c5b6c6_10 - - __osx >=11.0 - - libgfortran 5.* - - libgfortran5 >=13.3.0 - - llvm-openmp >=18.1.8 - - liblapack >=3.9.0,<4.0a0 - - libblas >=3.9.0,<4.0a0 - - metis >=5.1.0,<5.1.1.0a0 - - libscotch >=7.0.6,<7.0.7.0a0 - constrains: - - libopenblas * *openmp* - arch: arm64 - platform: osx - license: CECILL-C - size: 2675502 - timestamp: 1745406952432 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/ncurses-6.5-h5e97a16_3.conda - sha256: 2827ada40e8d9ca69a153a45f7fd14f32b2ead7045d3bbb5d10964898fe65733 - md5: 068d497125e4bf8a66bf707254fff5ae - depends: - - __osx >=11.0 - arch: arm64 - platform: osx - license: X11 AND BSD-3-Clause - size: 797030 - timestamp: 1738196177597 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/numpy-2.2.5-py312h7c1f314_0.conda - sha256: 982aed7df71ae0ca8bc396ae25d43fd9a4f2b15d18faca15d6c27e9efb3955be - md5: 24a41dacf9d624b069d54a6e92594540 - depends: - - __osx >=11.0 - - libblas >=3.9.0,<4.0a0 - - libcblas >=3.9.0,<4.0a0 - - libcxx >=18 - - liblapack >=3.9.0,<4.0a0 - - python >=3.12,<3.13.0a0 - - python >=3.12,<3.13.0a0 *_cpython - - python_abi 3.12.* *_cp312 - constrains: - - numpy-base <0a0 - arch: arm64 - platform: osx - license: BSD-3-Clause - license_family: BSD - size: 6498553 - timestamp: 1745119367238 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/openssl-3.5.0-h81ee809_1.conda - sha256: 73d366c1597a10bcd5f3604b5f0734b31c23225536e03782c6a13f9be9d01bff - md5: 5c7aef00ef60738a14e0e612cfc5bcde - depends: - - __osx >=11.0 - - ca-certificates - arch: arm64 - platform: osx - license: Apache-2.0 - license_family: Apache - size: 3064197 - timestamp: 1746223530698 -- conda: https://conda.anaconda.org/conda-forge/noarch/ply-3.11-pyhd8ed1ab_3.conda - sha256: bae453e5cecf19cab23c2e8929c6e30f4866d996a8058be16c797ed4b935461f - md5: fd5062942bfa1b0bd5e0d2a4397b099e - depends: - - python >=3.9 - license: BSD-3-Clause - license_family: BSD - size: 49052 - timestamp: 1733239818090 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/pyomo-6.9.2-py312hd8f9ff3_0.conda - sha256: 2bec4e5d361ab39d0b06bd026484cdd6399307d4c0e74807c2e2a2cca2c859f2 - md5: a7fc6b79deae458d5522b5103002120c - depends: - - __osx >=11.0 - - libcxx >=18 - - ply - - python >=3.12,<3.13.0a0 - - python >=3.12,<3.13.0a0 *_cpython - - python_abi 3.12.* *_cp312 - - setuptools - arch: arm64 - platform: osx - license: BSD-3-Clause - license_family: BSD - size: 7421923 - timestamp: 1744833268609 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/python-3.12.10-hc22306f_0_cpython.conda - sha256: 69aed911271e3f698182e9a911250b05bdf691148b670a23e0bea020031e298e - md5: c88f1a7e1e7b917d9c139f03b0960fea - depends: - - __osx >=11.0 - - bzip2 >=1.0.8,<2.0a0 - - libexpat >=2.7.0,<3.0a0 - - libffi >=3.4.6,<3.5.0a0 - - liblzma >=5.8.1,<6.0a0 - - libsqlite >=3.49.1,<4.0a0 - - libzlib >=1.3.1,<2.0a0 - - ncurses >=6.5,<7.0a0 - - openssl >=3.5.0,<4.0a0 - - readline >=8.2,<9.0a0 - - tk >=8.6.13,<8.7.0a0 - - tzdata - constrains: - - python_abi 3.12.* *_cp312 - arch: arm64 - platform: osx - license: Python-2.0 - size: 12932743 - timestamp: 1744323815320 -- conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.12-7_cp312.conda - build_number: 7 - sha256: a1bbced35e0df66cc713105344263570e835625c28d1bdee8f748f482b2d7793 - md5: 0dfcdc155cf23812a0c9deada86fb723 - constrains: - - python 3.12.* *_cpython - license: BSD-3-Clause - license_family: BSD - size: 6971 - timestamp: 1745258861359 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/readline-8.2-h1d1bf99_2.conda - sha256: 7db04684d3904f6151eff8673270922d31da1eea7fa73254d01c437f49702e34 - md5: 63ef3f6e6d6d5c589e64f11263dc5676 - depends: - - ncurses >=6.5,<7.0a0 - arch: arm64 - platform: osx - license: GPL-3.0-only - license_family: GPL - size: 252359 - timestamp: 1740379663071 -- conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.1.0-pyhff2d567_0.conda - sha256: 777d34ed359cedd5a5004c930077c101bb3b70e5fbb04d29da5058d75b0ba487 - md5: f6f72d0837c79eaec77661be43e8a691 - depends: - - python >=3.9 - license: MIT - license_family: MIT - size: 778484 - timestamp: 1746085063737 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/tinyxml2-11.0.0-ha1acc90_0.conda - sha256: bbd9294551ff727305f8335819c24d2490d5d79e0f3d90957992c39d2146093a - md5: 6778d917f88222e8f27af8ec5c41f277 - depends: - - __osx >=11.0 - - libcxx >=18 - arch: arm64 - platform: osx - license: Zlib - size: 122269 - timestamp: 1742246179980 -- conda: https://conda.anaconda.org/conda-forge/osx-arm64/tk-8.6.13-h5083fa2_1.conda - sha256: 72457ad031b4c048e5891f3f6cb27a53cb479db68a52d965f796910e71a403a8 - md5: b50a57ba89c32b62428b71a875291c9b - depends: - - libzlib >=1.2.13,<2.0.0a0 - arch: arm64 - platform: osx - license: TCL - license_family: BSD - size: 3145523 - timestamp: 1699202432999 -- conda: https://conda.anaconda.org/conda-forge/noarch/tzdata-2025b-h78e105d_0.conda - sha256: 5aaa366385d716557e365f0a4e9c3fca43ba196872abbbe3d56bb610d131e192 - md5: 4222072737ccff51314b5ece9c7d6f5a - license: LicenseRef-Public-Domain - size: 122968 - timestamp: 1742727099393 diff --git a/.CondaPkg/pixi.toml b/.CondaPkg/pixi.toml deleted file mode 100644 index 42e8ed6b17..0000000000 --- a/.CondaPkg/pixi.toml +++ /dev/null @@ -1,15 +0,0 @@ -[dependencies] -pyomo = "*" -casadi = ">3.3" - - [dependencies.python] - channel = "conda-forge" - build = "*cpython*" - version = ">=3.8,<4" - -[project] -name = ".CondaPkg" -platforms = ["osx-arm64"] -channels = ["conda-forge"] -channel-priority = "strict" -description = "automatically generated by CondaPkg.jl" 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 From f3128613277b1fa4e3ed721626fc59a68dda8c4d Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 5 May 2025 22:09:42 -0400 Subject: [PATCH 19/28] cleanup loose ends --- Project.toml | 2 +- src/systems/systems.jl | 5 --- test/extensions/Project.toml | 1 - test/extensions/jump_control.jl | 57 --------------------------------- 4 files changed, 1 insertion(+), 64 deletions(-) delete mode 100644 test/extensions/jump_control.jl diff --git a/Project.toml b/Project.toml index ff97c329a9..33280a5064 100644 --- a/Project.toml +++ b/Project.toml @@ -8,7 +8,6 @@ ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" -CasADi = "c49709b8-5c63-11e9-2fb2-69db5844192f" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" @@ -65,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" diff --git a/src/systems/systems.jl b/src/systems/systems.jl index 9da7249300..52f93afb9b 100644 --- a/src/systems/systems.jl +++ b/src/systems/systems.jl @@ -163,11 +163,6 @@ function __structural_simplify( end end -function toterm_auxsystems(system::ODESystem) - constraints = system.constraintsystem.constraints - -end - """ $(TYPEDSIGNATURES) diff --git a/test/extensions/Project.toml b/test/extensions/Project.toml index bd51ea5203..0e018b4a22 100644 --- a/test/extensions/Project.toml +++ b/test/extensions/Project.toml @@ -6,7 +6,6 @@ DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" diff --git a/test/extensions/jump_control.jl b/test/extensions/jump_control.jl deleted file mode 100644 index c0a071e31d..0000000000 --- a/test/extensions/jump_control.jl +++ /dev/null @@ -1,57 +0,0 @@ -using ModelingToolkit -using JuMP, InfiniteOpt -using DiffEqDevTools, DiffEqBase -using SimpleDiffEq -using HiGHS -const M = ModelingToolkit - -@testset "ODE Solution, no cost" begin - # Test solving without anything attached. - @parameters α=1.5 β=1.0 γ=3.0 δ=1.0 - M.@variables x(..) y(..) - t = M.t_nounits; D = M.D_nounits - - eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y(t), - D(y(t)) ~ -γ * y(t) + δ * x(t) * y(t)] - - tspan = (0.0, 1.0) - u0map = [x(t) => 4.0, y(t) => 2.0] - parammap = [α => 7.5, β => 4, γ => 8.0, δ => 5.0] - @mtkbuild sys = ODESystem(eqs, t) - - # Test explicit method. - jprob = JuMPControlProblem(sys, u0map, tspan, parammap, dt = 0.01) - @test num_constraints(jprob.model) == 2 # initials - jsol = solve(jprob, Ipopt.Optimizer, :Tsitouras5) - oprob = ODEProblem(sys, u0map, tspan, parammap) - osol = solve(oprob, SimpleTsit5(), adaptive = false) - @test jsol.sol.u ≈ osol.u - - # Implicit method. - jsol2 = solve(prob, Ipopt.Optimizer, :RK4) - osol2 = solve(oprob, SimpleRK4(), adaptive = false) - @test jsol2.sol.u ≈ osol2.u - - # With a constraint - @parameters α=1.5 β=1.0 γ=3.0 δ=1.0 - @variables x(..) y(..) - - eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y(t), - D(y(t)) ~ -γ * y(t) + δ * x(t) * y(t)] - - u0map = [] - tspan = (0.0, 1.0) - guess = [x(t) => 4.0, y(t) => 2.0] - constr = [x(0.6) ~ 3.5, x(0.3) ~ 7.0] - @mtkbuild lksys = ODESystem(eqs, t; constraints = constr) - - jprob = JuMPControlProblem(sys, u0map, tspan, parammap; guesses, dt = 0.01) - @test num_constraints(jprob.model) == 2 == num_variables(jprob.model) == 2 - jsol = solve(prob, HiGHS.Optimizer, :Tsitouras5) - sol = jsol.sol - @test sol(0.6)[1] ≈ 3.5 - @test sol(0.3)[1] ≈ 7.0 -end - -@testset "Optimal control problem" begin -end From ab2c777e24e1331ebabb2bc4e502d5ec0980ab5a Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 5 May 2025 22:11:17 -0400 Subject: [PATCH 20/28] bump CasADi compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 33280a5064..ccbe65c6ab 100644 --- a/Project.toml +++ b/Project.toml @@ -88,7 +88,7 @@ BifurcationKit = "0.4" BlockArrays = "1.1" BoundaryValueDiffEqAscher = "1.1.0" BoundaryValueDiffEqMIRK = "1.4.0" -CasADi = "1.0.0" +CasADi = "1.0.2" ChainRulesCore = "1" Combinatorics = "1" CommonSolve = "0.2.4" From b2845736c0ee0dc57e8d3ac13da63849f935ac27 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 6 May 2025 09:47:09 -0400 Subject: [PATCH 21/28] Update Project.toml --- test/downstream/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml index f64ed17de6..8de7753346 100644 --- a/test/downstream/Project.toml +++ b/test/downstream/Project.toml @@ -1,4 +1,5 @@ [deps] +CasADi = "c49709b8-5c63-11e9-2fb2-69db5844192f" ControlSystemsMTK = "687d7614-c7e5-45fc-bfc3-9ee385575c88" DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" From 189f8330c1937beb6d48c5905aab637105be0648 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 7 May 2025 11:29:56 -0400 Subject: [PATCH 22/28] bump Project toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ccbe65c6ab..f94d84b000 100644 --- a/Project.toml +++ b/Project.toml @@ -88,7 +88,7 @@ BifurcationKit = "0.4" BlockArrays = "1.1" BoundaryValueDiffEqAscher = "1.1.0" BoundaryValueDiffEqMIRK = "1.4.0" -CasADi = "1.0.2" +CasADi = "1.0.3" ChainRulesCore = "1" Combinatorics = "1" CommonSolve = "0.2.4" From 28a5e74c845595e5dce884c0637786b911900304 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 7 May 2025 12:31:07 -0400 Subject: [PATCH 23/28] remove is_explicit from jump --- ext/MTKInfiniteOptExt.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/ext/MTKInfiniteOptExt.jl b/ext/MTKInfiniteOptExt.jl index 7c374cfdaf..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) From 54c01f6de0b4d6e326fc4c98e7daf4f136e69628 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 7 May 2025 14:24:06 -0400 Subject: [PATCH 24/28] bump downstream casadi compat --- test/downstream/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml index 8de7753346..73e9b0becc 100644 --- a/test/downstream/Project.toml +++ b/test/downstream/Project.toml @@ -16,3 +16,4 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" [compat] ModelingToolkitStandardLibrary = "2.19" +CasADi = "1.0.3" From 47d9de4ed5b69761ce8f8bf75ee304f32d992800 Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 8 May 2025 17:03:04 -0400 Subject: [PATCH 25/28] fix: bump CasADi compat to 1.0.5 --- Project.toml | 2 +- test/downstream/Project.toml | 2 +- test/downstream/dynamic_optimization.jl | 1 + 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index f94d84b000..a38fdc99cc 100644 --- a/Project.toml +++ b/Project.toml @@ -88,7 +88,7 @@ BifurcationKit = "0.4" BlockArrays = "1.1" BoundaryValueDiffEqAscher = "1.1.0" BoundaryValueDiffEqMIRK = "1.4.0" -CasADi = "1.0.3" +CasADi = "1.0.5" ChainRulesCore = "1" Combinatorics = "1" CommonSolve = "0.2.4" diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml index 73e9b0becc..9b3f79487b 100644 --- a/test/downstream/Project.toml +++ b/test/downstream/Project.toml @@ -16,4 +16,4 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" [compat] ModelingToolkitStandardLibrary = "2.19" -CasADi = "1.0.3" +CasADi = "1.0.5" diff --git a/test/downstream/dynamic_optimization.jl b/test/downstream/dynamic_optimization.jl index abb534ea0b..4b3b4edeb6 100644 --- a/test/downstream/dynamic_optimization.jl +++ b/test/downstream/dynamic_optimization.jl @@ -7,6 +7,7 @@ using Ipopt using DataInterpolations using CasADi +import DiffEqBase: solve const M = ModelingToolkit @testset "ODE Solution, no cost" begin From a9f6aea5e381299dd60c7958b37c755d2719861d Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 8 May 2025 17:28:18 -0400 Subject: [PATCH 26/28] rebase: move tests --- test/downstream/Project.toml | 2 -- test/extensions/Project.toml | 4 ++++ test/{downstream => extensions}/dynamic_optimization.jl | 0 3 files changed, 4 insertions(+), 2 deletions(-) rename test/{downstream => extensions}/dynamic_optimization.jl (100%) diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml index 9b3f79487b..f64ed17de6 100644 --- a/test/downstream/Project.toml +++ b/test/downstream/Project.toml @@ -1,5 +1,4 @@ [deps] -CasADi = "c49709b8-5c63-11e9-2fb2-69db5844192f" ControlSystemsMTK = "687d7614-c7e5-45fc-bfc3-9ee385575c88" DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -16,4 +15,3 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" [compat] ModelingToolkitStandardLibrary = "2.19" -CasADi = "1.0.5" diff --git a/test/extensions/Project.toml b/test/extensions/Project.toml index 0e018b4a22..e4f59dbec7 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.5" diff --git a/test/downstream/dynamic_optimization.jl b/test/extensions/dynamic_optimization.jl similarity index 100% rename from test/downstream/dynamic_optimization.jl rename to test/extensions/dynamic_optimization.jl From 7d561a3a3a703bf51a50854f326626fe476e7ffd Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 8 May 2025 19:25:53 -0400 Subject: [PATCH 27/28] rename file in runtests --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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") From e1e0b90a16399bf9ccf74b89b8649e07f29710e5 Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 8 May 2025 20:18:18 -0400 Subject: [PATCH 28/28] bump Casadi compat --- Project.toml | 2 +- test/extensions/Project.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index a38fdc99cc..9538bb0894 100644 --- a/Project.toml +++ b/Project.toml @@ -88,7 +88,7 @@ BifurcationKit = "0.4" BlockArrays = "1.1" BoundaryValueDiffEqAscher = "1.1.0" BoundaryValueDiffEqMIRK = "1.4.0" -CasADi = "1.0.5" +CasADi = "1.0.6" ChainRulesCore = "1" Combinatorics = "1" CommonSolve = "0.2.4" diff --git a/test/extensions/Project.toml b/test/extensions/Project.toml index e4f59dbec7..9f43e6f4a4 100644 --- a/test/extensions/Project.toml +++ b/test/extensions/Project.toml @@ -28,4 +28,4 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] -CasADi = "1.0.5" +CasADi = "1.0.6"