Skip to content

Commit 0c3f363

Browse files
committed
reemove solver info from tutorial
1 parent d0f2243 commit 0c3f363

File tree

1 file changed

+19
-46
lines changed

1 file changed

+19
-46
lines changed
Lines changed: 19 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
# Solving Dynamic Optimization Problems
2+
23
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.
34

45
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.
@@ -10,8 +11,8 @@ D = ModelingToolkit.D_nounits
1011
1112
@parameters h_c m₀ h₀ g₀ D_c c Tₘ m_c
1213
@variables begin
13-
h(..)
14-
v(..)
14+
h(..)
15+
v(..)
1516
m(..), [bounds = (m_c, 1)]
1617
T(..), [input = true, bounds = (0, Tₘ)]
1718
end
@@ -20,8 +21,8 @@ drag(h, v) = D_c * v^2 * exp(-h_c * (h - h₀) / h₀)
2021
gravity(h) = g₀ * (h₀ / h)
2122
2223
eqs = [D(h(t)) ~ v(t),
23-
D(v(t)) ~ (T(t) - drag(h(t), v(t))) / m(t) - gravity(h(t)),
24-
D(m(t)) ~ -T(t) / c]
24+
D(v(t)) ~ (T(t) - drag(h(t), v(t))) / m(t) - gravity(h(t)),
25+
D(m(t)) ~ -T(t) / c]
2526
2627
(ts, te) = (0.0, 0.2)
2728
costs = [-h(te)]
@@ -35,24 +36,28 @@ pmap = [
3536
g₀ => 1, m₀ => 1.0, h_c => 500, c => 0.5 * √(g₀ * h₀), D_c => 0.5 * 620 * m₀ / g₀,
3637
Tₘ => 3.5 * g₀ * m₀, T(t) => 0.0, h₀ => 1, m_c => 0.6]
3738
```
39+
3840
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.
3941

4042
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).
4143

4244
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.
45+
4346
```@example dynamic_opt
4447
using InfiniteOpt, Ipopt, DiffEqDevTools
4548
jprob = JuMPDynamicOptProblem(rocket, [u0map; pmap], (ts, te); dt = 0.001)
4649
jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRadauIIA5()));
4750
```
51+
4852
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.
4953

5054
Let's plot the final solution and the controller here:
55+
5156
```@example dynamic_opt
5257
using CairoMakie
5358
fig = Figure(resolution = (800, 400))
54-
ax1 = Axis(fig[1,1], title = "Rocket trajectory", xlabel = "Time")
55-
ax2 = Axis(fig[1,2], title = "Control trajectory", xlabel = "Time")
59+
ax1 = Axis(fig[1, 1], title = "Rocket trajectory", xlabel = "Time")
60+
ax2 = Axis(fig[1, 2], title = "Control trajectory", xlabel = "Time")
5661
5762
for u in unknowns(rocket)
5863
lines!(ax1, jsol.sol.t, jsol.sol[u], label = string(u))
@@ -64,12 +69,14 @@ fig
6469
```
6570

6671
### Free final time problems
72+
6773
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.
6874

6975
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.
76+
7077
```@example dynamic_opt
7178
@variables begin
72-
x(..)
79+
x(..)
7380
v(..)
7481
u(..), [bounds = (-1.0, 1.0), input = true]
7582
end
@@ -88,6 +95,7 @@ u0map = [x(t) => 1.0, v(t) => 0.0]
8895
tspan = (0.0, tf)
8996
parammap = [u(t) => 0.0, tf => 1.0]
9097
```
98+
9199
The `tf` mapping in the parameter map is treated as an initial guess.
92100

93101
Please note that, at the moment, free final time problems cannot support constraints defined at definite time values, like `x(3) ~ 2`.
@@ -97,16 +105,18 @@ Please note that, at the moment, free final time problems cannot support constra
97105
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.
98106

99107
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:
108+
100109
```@example dynamic_opt
101110
iprob = InfiniteOptDynamicOptProblem(block, [u0map; parammap], tspan; steps = 100)
102111
isol = solve(iprob, InfiniteOptCollocation(Ipopt.Optimizer));
103112
```
104113

105114
Let's plot the final solution and the controller here:
115+
106116
```@example dynamic_opt
107117
fig = Figure(resolution = (800, 400))
108-
ax1 = Axis(fig[1,1], title = "Block trajectory", xlabel = "Time")
109-
ax2 = Axis(fig[1,2], title = "Control trajectory", xlabel = "Time")
118+
ax1 = Axis(fig[1, 1], title = "Block trajectory", xlabel = "Time")
119+
ax2 = Axis(fig[1, 2], title = "Control trajectory", xlabel = "Time")
110120
111121
for u in unknowns(block)
112122
lines!(ax1, isol.sol.t, isol.sol[u], label = string(u))
@@ -116,40 +126,3 @@ axislegend(ax1)
116126
axislegend(ax2)
117127
fig
118128
```
119-
120-
### Solvers
121-
Currently 4 backends are exposed for solving dynamic optimization problems using collocation: JuMP, InfiniteOpt, CasADi, and Pyomo.
122-
123-
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:
124-
```
125-
JuMPCollocation(Ipopt.Optimizer, constructRK4())
126-
CasADiCollocation("ipopt", constructRK4())
127-
```
128-
129-
**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.
130-
131-
**Pyomo** and **InfiniteOpt** each have their own built-in collocation methods.
132-
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.
133-
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)`.
134-
135-
Some examples of the latter two collocations:
136-
```julia
137-
PyomoCollocation("ipopt", LagrangeRadau(2))
138-
InfiniteOptCollocation(Ipopt.Optimizer, OrthogonalCollocation(3))
139-
```
140-
141-
```@docs; canonical = false
142-
JuMPCollocation
143-
InfiniteOptCollocation
144-
CasADiCollocation
145-
PyomoCollocation
146-
solve(::AbstractDynamicOptProblem)
147-
```
148-
149-
### Problem constructors
150-
```@docs; canonical = false
151-
JuMPDynamicOptProblem
152-
InfiniteOptDynamicOptProblem
153-
CasADiDynamicOptProblem
154-
PyomoDynamicOptProblem
155-
```

0 commit comments

Comments
 (0)