Skip to content

feat: PyomoDynamicOptProblem, docs #3676

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 22 commits into from
Jun 9, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
22 commits
Select commit Hold shift + click to select a range
2ccdbba
refactor: add definitions of shared functions to optimal_control_inte…
vyudu May 15, 2025
df1a0e4
refactor: aduse optimal control itnerface
vyudu May 16, 2025
912d71b
refactor: add interface functions for CasADi
vyudu May 16, 2025
4ca744a
correctly implement interface
vyudu May 17, 2025
8d268e9
refactor: move set_variable_bounds to interface function
vyudu May 19, 2025
b5b0f44
refactor: centralize substitute_differentiasl and substitute_integral
vyudu May 19, 2025
072824c
chore: move the docstrings to the interface file
vyudu May 20, 2025
7fda5b0
test: add Pyomo tests
vyudu May 20, 2025
a275b33
refactor: use lowered_var instead of manual substitution
vyudu May 22, 2025
3c591c0
feat: PyomoDynamicOtpPRoblem
vyudu May 26, 2025
33ad60a
fix: fix constraint and cost transcription for pyomo
vyudu May 27, 2025
5988c81
fix: pyomo needs one fewer finite element than # tsteps
vyudu May 28, 2025
5f33573
docs: dynamic optimization docs
vyudu May 29, 2025
f7d246f
docs: make examples dynamic
vyudu May 29, 2025
10719ae
fix: update ext/tests to v10
vyudu May 29, 2025
7daecfb
docs: update docs to v10
vyudu May 29, 2025
f33f8c7
fix: error when using collocation on fft problems
vyudu Jun 2, 2025
d5a49b0
reemove solver info from tutorial
vyudu Jun 2, 2025
8750986
update Project.toml, fix substitution bug
vyudu Jun 2, 2025
97bb8ff
add Pyomo to extensions
vyudu Jun 2, 2025
a818178
fix: fix pyomo problems on lts
vyudu Jun 5, 2025
639ed6a
fix: use symbolic indexing for solutions of pprobs
vyudu Jun 5, 2025
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,7 @@ DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
FMI = "14a09403-18e3-468f-ad8a-74f8dda2d9ac"
InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57"
LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
Pyomo = "0e8e1daf-01b5-4eba-a626-3897743a3816"

[extensions]
MTKBifurcationKitExt = "BifurcationKit"
Expand All @@ -81,6 +82,7 @@ MTKDeepDiffsExt = "DeepDiffs"
MTKFMIExt = "FMI"
MTKInfiniteOptExt = "InfiniteOpt"
MTKLabelledArraysExt = "LabelledArrays"
MTKPyomoDynamicOptExt = "Pyomo"

[compat]
ADTypes = "1.14.0"
Expand Down Expand Up @@ -140,6 +142,7 @@ OrdinaryDiffEqCore = "1.15.0"
OrdinaryDiffEqDefault = "1.2"
OrdinaryDiffEqNonlinearSolve = "1.5.0"
PrecompileTools = "1"
Pyomo = "0.1.0"
REPL = "1"
RecursiveArrayTools = "3.26"
Reexport = "0.2, 1"
Expand Down
3 changes: 3 additions & 0 deletions docs/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,14 @@ BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
FMI = "14a09403-18e3-468f-ad8a-74f8dda2d9ac"
FMIZoo = "724179cf-c260-40a9-bd27-cccc6fe2f195"
InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57"
Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
Expand Down
41 changes: 41 additions & 0 deletions docs/src/API/dynamic_opt.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,41 @@
### Solvers

Currently 4 backends are exposed for solving dynamic optimization problems using collocation: JuMP, InfiniteOpt, CasADi, and Pyomo.

Please note that there are differences in how to construct the collocation solver for the different cases. For example, the Python based ones, CasADi and Pyomo, expect the solver to be passed in as a string (CasADi and Pyomo come pre-loaded with Ipopt, but other solvers may need to be manually installed using `pip` or `conda`), while JuMP/InfiniteOpt expect the optimizer object to be passed in directly:

```
JuMPCollocation(Ipopt.Optimizer, constructRK4())
CasADiCollocation("ipopt", constructRK4())
```

**JuMP** and **CasADi** collocation require an ODE tableau to be passed in. These can be constructed by calling the `constructX()` functions from DiffEqDevTools. The list of tableaus can be found [here](https://docs.sciml.ai/DiffEqDevDocs/dev/internals/tableaus/). If none is passed in, both solvers will default to using Radau second-order with five collocation points.

**Pyomo** and **InfiniteOpt** each have their own built-in collocation methods.

1. **InfiniteOpt**: The list of InfiniteOpt collocation methods can be found [in the table on this page](https://infiniteopt.github.io/InfiniteOpt.jl/stable/guide/derivative/). If none is passed in, the solver defaults to `FiniteDifference(Backward())`, which is effectively implicit Euler.
2. **Pyomo**: The list of Pyomo collocation methods can be found [at the bottom of this page](https://github.com/SciML/Pyomo.jl). If none is passed in, the solver defaults to a `LagrangeRadau(3)`.

Some examples of the latter two collocations:

```julia
PyomoCollocation("ipopt", LagrangeRadau(2))
InfiniteOptCollocation(Ipopt.Optimizer, OrthogonalCollocation(3))
```

```@docs; canonical = false
JuMPCollocation
InfiniteOptCollocation
CasADiCollocation
PyomoCollocation
solve(::AbstractDynamicOptProblem)
```

### Problem constructors

```@docs; canonical = false
JuMPDynamicOptProblem
InfiniteOptDynamicOptProblem
CasADiDynamicOptProblem
PyomoDynamicOptProblem
```
128 changes: 128 additions & 0 deletions docs/src/tutorials/dynamic_optimization.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,128 @@
# Solving Dynamic Optimization Problems

Systems in ModelingToolkit.jl can be directly converted to dynamic optimization or optimal control problems. In such systems, one has one or more input variables that are externally controlled to control the dynamics of the system. A dynamic optimization solves for the optimal time trajectory of the input variables in order to maximize or minimize a desired objective function. For example, a car driver might like to know how to step on the accelerator if the goal is to finish a race while using the least gas.

To begin, let us take a rocket launch example. The input variable here is the thrust exerted by the engine. The rocket state is described by its current height, mass, and velocity. The mass decreases as the rocket loses fuel while thrusting.

```@example dynamic_opt
using ModelingToolkit
t = ModelingToolkit.t_nounits
D = ModelingToolkit.D_nounits

@parameters h_c m₀ h₀ g₀ D_c c Tₘ m_c
@variables begin
h(..)
v(..)
m(..), [bounds = (m_c, 1)]
T(..), [input = true, bounds = (0, Tₘ)]
end

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 = System(eqs, t; costs, constraints = cons)
rocket = mtkcompile(rocket, inputs = [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]
```

What we would like to optimize here is the final height of the rocket. We do this by providing a vector of expressions corresponding to the costs. By default, the sense of the optimization is to minimize the provided cost. So to maximize the rocket height at the final time, we write `-h(te)` as the cost.

Now we can construct a problem and solve it. Let us use JuMP as our backend here. Note that the package trigger is actually [InfiniteOpt](https://infiniteopt.github.io/InfiniteOpt.jl/stable/), and not JuMP - this package includes JuMP but is designed for optimization on function spaces. Additionally we need to load the solver package - we will use [Ipopt](https://github.com/jump-dev/Ipopt.jl) here (a good choice in general).

Here we have also loaded DiffEqDevTools because we will need to construct the ODE tableau. This is only needed if one desires a custom ODE tableau for the collocation - by default the solver will use RadauIIA5.

```@example dynamic_opt
using InfiniteOpt, Ipopt, DiffEqDevTools
jprob = JuMPDynamicOptProblem(rocket, [u0map; pmap], (ts, te); dt = 0.001)
jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRadauIIA5()));
```

The solution has three fields: `jsol.sol` is the ODE solution for the states, `jsol.input_sol` is the ODE solution for the inputs, and `jsol.model` is the wrapped model that we can use to query things like objective and constraint residuals.

Let's plot the final solution and the controller here:

```@example dynamic_opt
using CairoMakie
fig = Figure(resolution = (800, 400))
ax1 = Axis(fig[1, 1], title = "Rocket trajectory", xlabel = "Time")
ax2 = Axis(fig[1, 2], title = "Control trajectory", xlabel = "Time")

for u in unknowns(rocket)
lines!(ax1, jsol.sol.t, jsol.sol[u], label = string(u))
end
lines!(ax2, jsol.input_sol, label = "Thrust")
axislegend(ax1)
axislegend(ax2)
fig
```

### Free final time problems

There are additionally a class of dynamic optimization problems where we would like to know how to control our system to achieve something in the least time. Such problems are called free final time problems, since the final time is unknown. To model these problems in ModelingToolkit, we declare the final time as a parameter.

Below we have a model system called the double integrator. We control the acceleration of a block in order to reach a desired destination in the least time.

```@example dynamic_opt
@variables begin
x(..)
v(..)
u(..), [bounds = (-1.0, 1.0), input = true]
end

@parameters tf

constr = [v(tf) ~ 0, x(tf) ~ 0]
cost = [tf] # Minimize time

@named block = System(
[D(x(t)) ~ v(t), D(v(t)) ~ u(t)], t; costs = cost, constraints = constr)

block = mtkcompile(block; inputs = [u(t)])

u0map = [x(t) => 1.0, v(t) => 0.0]
tspan = (0.0, tf)
parammap = [u(t) => 0.0, tf => 1.0]
```

The `tf` mapping in the parameter map is treated as an initial guess.

Please note that, at the moment, free final time problems cannot support constraints defined at definite time values, like `x(3) ~ 2`.

!!! warning

The Pyomo collocation methods (LagrangeRadau, LagrangeLegendre) currently are bugged for free final time problems. Strongly suggest using BackwardEuler() for such problems when using Pyomo as the backend.

When declaring the problem in this case we need to provide the number of steps, since dt can't be known in advanced. Let's solve plot our final solution and the controller for the block, using InfiniteOpt as the backend:

```@example dynamic_opt
iprob = InfiniteOptDynamicOptProblem(block, [u0map; parammap], tspan; steps = 100)
isol = solve(iprob, InfiniteOptCollocation(Ipopt.Optimizer));
```

Let's plot the final solution and the controller here:

```@example dynamic_opt
fig = Figure(resolution = (800, 400))
ax1 = Axis(fig[1, 1], title = "Block trajectory", xlabel = "Time")
ax2 = Axis(fig[1, 2], title = "Control trajectory", xlabel = "Time")

for u in unknowns(block)
lines!(ax1, isol.sol.t, isol.sol[u], label = string(u))
end
lines!(ax2, isol.input_sol, label = "Acceleration")
axislegend(ax1)
axislegend(ax2)
fig
```
Loading
Loading