Solve a problem
In this tutorial, we consider the Beam problem. First, import the OptimalControlProblems package to access the problem:
using OptimalControlProblemsSolving from an OptimalControl model
First, import the problem:
using OptimalControl
docp = beam(OptimalControlBackend())
nlp = nlp_model(docp)The nlp model represents the nonlinear programming problem (NLP) obtained after discretising the optimal control problem (OCP). See the Introduction page for details. The model is an ADNLPModels.ADNLPModel, which provides automatic differentiation (AD)-based models that follow the NLPModels.jl API.
We can then solve the problem using, for instance, NLPModelsIpopt.ipopt:
using NLPModelsIpopt
# Solve the model
sol = NLPModelsIpopt.ipopt(
nlp;
print_level=5,
tol=1e-8,
mu_strategy="adaptive",
sb="yes",
)This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.1.
Number of nonzeros in equality constraint Jacobian...: 4004
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 501
Total number of variables............................: 1503
variables with only lower bounds: 0
variables with lower and upper bounds: 1002
variables with only upper bounds: 0
Total number of equality constraints.................: 1004
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.0000000e-02 1.10e+00 2.62e-13 0.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 3.6382396e-01 7.50e-01 1.72e+00 -0.9 3.39e+00 - 5.94e-01 3.18e-01h 1
2 1.1944708e+01 7.50e-03 1.82e+00 -0.8 7.62e+00 - 6.98e-01 9.90e-01f 1
3 1.2010093e+01 7.48e-05 2.84e-01 -1.9 2.53e-01 - 9.98e-01 9.90e-01h 1
4 1.0698046e+01 5.30e-07 1.46e+00 -2.4 1.80e+00 - 9.94e-01 9.93e-01f 1
5 9.2863700e+00 1.11e-16 2.37e+01 -3.2 9.06e-01 - 9.99e-01 1.00e+00f 1
6 8.9473960e+00 1.11e-16 2.67e+03 -4.2 2.48e-01 - 1.00e+00 7.69e-01f 1
7 8.9328684e+00 1.11e-16 6.58e+03 -5.1 2.35e-01 - 1.00e+00 2.34e-01f 1
8 8.9059261e+00 1.11e-16 2.68e+03 -4.7 2.86e-01 - 3.16e-01 6.56e-01f 1
9 8.9052608e+00 1.11e-16 2.64e+03 -10.6 1.81e-01 - 9.64e-02 3.95e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 8.8997178e+00 1.11e-16 1.53e+03 -6.7 1.67e-01 - 1.98e-01 3.58e-01f 1
11 8.8995100e+00 2.22e-16 1.57e+03 -10.8 1.74e-01 - 1.05e-01 2.18e-02f 1
12 8.8968682e+00 2.22e-16 1.34e+03 -6.0 9.89e-02 - 5.06e-01 2.49e-01f 1
13 8.8965717e+00 2.22e-16 1.31e+03 -11.0 1.65e-01 - 9.30e-02 4.03e-02f 1
14 8.8940168e+00 2.22e-16 8.08e+02 -6.4 1.13e-01 - 1.82e-01 3.42e-01f 1
15 8.8938704e+00 2.22e-16 8.31e+02 -11.0 1.29e-01 - 1.86e-01 3.28e-02f 1
16 8.8927053e+00 2.22e-16 7.10e+02 -5.9 6.43e-02 - 7.11e-01 2.26e-01f 1
17 8.8900801e+00 1.11e-16 5.70e+01 -5.6 9.74e-02 - 3.96e-01 9.85e-01f 1
18 8.8900330e+00 1.11e-16 1.10e+02 -11.0 4.06e-02 - 3.60e-01 6.68e-02f 1
19 8.8898907e+00 2.22e-16 2.70e+01 -5.5 2.95e-02 - 6.93e-01 3.13e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 8.8896270e+00 1.11e-16 2.79e+01 -6.0 3.55e-02 - 4.27e-01 5.86e-01f 1
21 8.8895382e+00 1.11e-16 1.02e+01 -5.7 1.09e-02 - 9.98e-01 7.56e-01f 1
22 8.8894231e+00 1.11e-16 5.46e+01 -6.3 2.89e-02 - 9.68e-01 6.10e-01f 1
23 8.8893916e+00 1.11e-16 7.11e-15 -6.0 6.02e-03 - 1.00e+00 1.00e+00f 1
24 8.8893248e+00 1.11e-16 1.78e+00 -6.5 1.48e-02 - 9.95e-01 9.70e-01f 1
25 8.8893119e+00 2.22e-16 2.19e+01 -7.7 1.43e-02 - 1.00e+00 3.06e-01f 1
26 8.8892891e+00 1.11e-16 4.32e+00 -8.0 1.69e-02 - 1.00e+00 8.10e-01f 1
27 8.8892829e+00 1.11e-16 1.73e-03 -8.7 8.14e-03 - 1.00e+00 1.00e+00f 1
28 8.8892819e+00 1.11e-16 1.54e-02 -9.6 3.13e-03 - 1.00e+00 9.96e-01f 1
29 8.8892817e+00 1.11e-16 7.11e-15 -10.8 1.42e-03 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 29
(scaled) (unscaled)
Objective...............: 8.8892816510267849e+00 8.8892816510267849e+00
Dual infeasibility......: 7.1054273576010019e-15 7.1054273576010019e-15
Constraint violation....: 1.1102230246251565e-16 1.1102230246251565e-16
Variable bound violation: 9.9998963970504562e-09 9.9998963970504562e-09
Complementarity.........: 1.1068094190397985e-09 1.1068094190397985e-09
Overall NLP error.......: 1.1068094190397985e-09 1.1068094190397985e-09
Number of objective function evaluations = 30
Number of objective gradient evaluations = 30
Number of equality constraint evaluations = 30
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 30
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 29
Total seconds in IPOPT = 0.079
EXIT: Optimal Solution Found.To get the number of iterations from the NLP solution:
sol.iter29and the objective value:
sol.objective8.889281651026785To recover the state, control, and costate, we recommend building an optimal control solution and using the associated getters (you can also retrieve the number of iterations and the objective value from the OCP solution):
ocp_sol = build_ocp_solution(docp, sol)
t = time_grid(ocp_sol) # t0, ..., tN = tf
x = state(ocp_sol) # function of time
u = control(ocp_sol) # function of time
p = costate(ocp_sol) # function of time
o = objective(ocp_sol) # scalar objective value
i = iterations(ocp_sol) # number of iteration
tf = t[end]
println("tf = ", tf)
println("x(tf) = ", x(tf))
println("u(tf) = ", u(tf))
println("p(tf) = ", p(tf))
println("objective: ", o)
println("iterations: ", i)tf = 1.0
x(tf) = [6.983483423835114e-36, -1.0]
u(tf) = -6.644778101442487
p(tf) = [44.44721421847331, -13.37845063459404]
objective: 8.889281651026785
iterations: 29If the problem includes additional optimisation variables, such as the final time, you can retrieve them with:
v = variable(ocp_sol)From ocp_sol you can also plot the state, control, and costate trajectories. For more details about plotting optimal control problems, see the Plot Manual.
using Plots
plt = plot(ocp_sol; color=1, size=(800, 700), control_style=(label="OptimalControl", ))
Solving from a JuMP model
First, import the JuMP model:
using JuMP
nlp = beam(JuMPBackend())A JuMP Model
├ solver: none
├ objective_sense: MIN_SENSE
│ └ objective_function_type: QuadExpr
├ num_variables: 1503
├ num_constraints: 3008
│ ├ AffExpr in MOI.EqualTo{Float64}: 1004
│ ├ VariableRef in MOI.GreaterThan{Float64}: 1002
│ └ VariableRef in MOI.LessThan{Float64}: 1002
└ Names registered in the model
└ :u, :x1, :x2, :∂x1, :∂x2We can then solve the problem using the JuMP.optimize! function:
using Ipopt
# Set the optimiser
set_optimizer(nlp, Ipopt.Optimizer)
# Set optimiser attributes
set_optimizer_attribute(nlp, "tol", 1e-8)
set_optimizer_attribute(nlp, "mu_strategy", "adaptive")
set_optimizer_attribute(nlp, "linear_solver", "mumps")
set_optimizer_attribute(nlp, "sb", "yes")
# Solve the model
optimize!(nlp)This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.1.
Number of nonzeros in equality constraint Jacobian...: 4004
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 501
Total number of variables............................: 1503
variables with only lower bounds: 0
variables with lower and upper bounds: 1002
variables with only upper bounds: 0
Total number of equality constraints.................: 1004
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.0000000e-02 1.10e+00 2.54e-13 0.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 3.6382396e-01 7.50e-01 1.72e+00 -0.9 3.39e+00 - 5.94e-01 3.18e-01h 1
2 1.1944708e+01 7.50e-03 1.82e+00 -0.8 7.62e+00 - 6.98e-01 9.90e-01f 1
3 1.2010093e+01 7.48e-05 2.84e-01 -1.9 2.53e-01 - 9.98e-01 9.90e-01h 1
4 1.0698046e+01 5.30e-07 1.46e+00 -2.4 1.80e+00 - 9.94e-01 9.93e-01f 1
5 9.2863700e+00 9.45e-17 2.37e+01 -3.2 9.06e-01 - 9.99e-01 1.00e+00f 1
6 8.9473960e+00 9.02e-17 2.67e+03 -4.2 2.48e-01 - 1.00e+00 7.69e-01f 1
7 8.9328684e+00 1.41e-16 6.58e+03 -5.1 2.35e-01 - 1.00e+00 2.34e-01f 1
8 8.9059261e+00 1.08e-16 2.68e+03 -4.7 2.86e-01 - 3.16e-01 6.56e-01f 1
9 8.9052608e+00 1.48e-16 2.64e+03 -10.6 1.81e-01 - 9.64e-02 3.95e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 8.8997178e+00 1.77e-16 1.53e+03 -6.7 1.67e-01 - 1.98e-01 3.58e-01f 1
11 8.8995100e+00 2.03e-16 1.57e+03 -10.8 1.74e-01 - 1.05e-01 2.18e-02f 1
12 8.8968682e+00 1.80e-16 1.34e+03 -6.0 9.89e-02 - 5.06e-01 2.49e-01f 1
13 8.8965717e+00 2.28e-16 1.31e+03 -11.0 1.65e-01 - 9.30e-02 4.03e-02f 1
14 8.8940168e+00 1.65e-16 8.08e+02 -6.4 1.13e-01 - 1.82e-01 3.42e-01f 1
15 8.8938704e+00 2.53e-16 8.31e+02 -11.0 1.29e-01 - 1.86e-01 3.28e-02f 1
16 8.8927053e+00 2.46e-16 7.10e+02 -5.9 6.43e-02 - 7.11e-01 2.26e-01f 1
17 8.8900778e+00 1.04e-16 5.51e+01 -5.6 9.75e-02 - 3.97e-01 9.86e-01f 1
18 8.8900312e+00 1.47e-16 1.08e+02 -11.0 4.05e-02 - 3.59e-01 6.63e-02f 1
19 8.8898898e+00 1.45e-16 2.71e+01 -5.5 2.97e-02 - 6.90e-01 3.10e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 8.8896196e+00 1.28e-16 3.13e+01 -6.0 3.56e-02 - 4.26e-01 6.05e-01f 1
21 8.8895320e+00 1.04e-16 9.66e+00 -5.7 1.08e-02 - 1.00e+00 7.59e-01f 1
22 8.8894202e+00 1.42e-16 5.53e+01 -6.3 2.88e-02 - 9.79e-01 6.08e-01f 1
23 8.8893851e+00 1.08e-16 1.47e-14 -6.0 6.56e-03 - 1.00e+00 1.00e+00f 1
24 8.8893214e+00 8.93e-17 2.60e-01 -6.5 1.41e-02 - 9.96e-01 1.00e+00f 1
25 8.8893107e+00 1.66e-16 2.21e+01 -8.4 1.51e-02 - 1.00e+00 2.72e-01f 1
26 8.8892856e+00 1.03e-16 1.82e-02 -8.1 1.71e-02 - 1.00e+00 9.99e-01f 1
27 8.8892828e+00 1.03e-16 1.22e+00 -9.5 5.42e-03 - 1.00e+00 8.33e-01f 1
28 8.8892818e+00 1.03e-16 7.11e-15 -8.9 2.98e-03 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 28
(scaled) (unscaled)
Objective...............: 8.8892818464565817e+00 8.8892818464565817e+00
Dual infeasibility......: 7.1054273576010019e-15 7.1054273576010019e-15
Constraint violation....: 1.0321604682062002e-16 1.0321604682062002e-16
Variable bound violation: 9.9671177561377888e-09 9.9671177561377888e-09
Complementarity.........: 4.1493706604063998e-09 4.1493706604063998e-09
Overall NLP error.......: 4.1493706604063998e-09 4.1493706604063998e-09
Number of objective function evaluations = 29
Number of objective gradient evaluations = 29
Number of equality constraint evaluations = 29
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 1
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 1
Total seconds in IPOPT = 0.056
EXIT: Optimal Solution Found.To get the number of iterations:
barrier_iterations(nlp)28To get the objective value:
objective_value(nlp)8.889281846456582To get the time grid, state, control, and costate, OptimalControlProblems provides the following getters (you have also similar getters to retrieve the number of iterations and the objective value):
problem = :beam
t = time_grid(problem, nlp) # t0, ..., tN = tf
x = state(problem, nlp) # function of time
u = control(problem, nlp) # function of time
p = costate(problem, nlp) # function of time
o = objective(problem, nlp) # scalar objective value
i = iterations(problem, nlp) # number of iteration
tf = t[end]
println("tf = ", tf)
println("x(tf) = ", x(tf))
println("u(tf) = ", u(tf))
println("p(tf) = ", p(tf))
println("objective: ", o)
println("iterations: ", i)tf = 1.0
x(tf) = [-4.454407807932683e-34, -1.0]
u(tf) = -6.644786948253476
p(tf) = [-44.447394877959965, 13.378468969723418]
objective: 8.889281846456582
iterations: 28If the problem includes additional optimisation variables, such as the final time, you can retrieve them with:
v = variable(problem, nlp)We can add the state, costate, and control to the plot:
n = length(metadata[problem][:state_name]) # dimension of the state
m = length(metadata[problem][:control_name]) # dimension of the control
for i in 1:n # state
plot!(plt[i], t, t -> x(t)[i]; color=2, linestyle=:dash, label=:none)
end
for i in 1:n # costate
plot!(plt[n+i], t, t -> -p(t)[i]; color=2, linestyle=:dash, label=:none)
end
for i in 1:m # control
plot!(plt[2n+i], t, t -> u(t)[i]; color=2, linestyle=:dash, label="JuMP")
end