NLP and DOCP manipulations

We describe here some more advanced operations related to the discretized optimal control problem. When calling solve(ocp) three steps are performed internally:

first, the OCP is discretized into a DOCP (a nonlinear optimization problem) with direct_transcription,
then, this DOCP is solved (with the internal function solve_docp),
finally, a functional solution of the OCP is rebuilt from the solution of the discretized problem, with OptimalControlSolution.

These steps can also be done separately, for instance if you want to use your own NLP solver.

Let us load the packages.

using OptimalControl
using Plots

Definition of the optimal control problem

We define a test problem

ocp = @def begin

    t ∈ [0, 1], time
    x ∈ R², state
    u ∈ R, control

    x(0) == [ -1, 0 ]
    x(1) == [ 0, 0 ]

    ẋ(t) == [ x₂(t), u(t) ]

    ∫( 0.5u(t)^2 ) → min

end

Discretization and NLP problem

We discretize the problem.

docp, nlp = direct_transcription(ocp)

The DOCP contains information related to the transcription, including a copy of the original OCP, and the NLP is the resulting discretized problem, in our case an ADNLPModel.

We can now use the solver of our choice to solve it.

Resolution of the NLP problem

For a first example we use the ipopt solver from NLPModelsIpopt.jl package to solve the NLP problem.

using NLPModelsIpopt

nlp_sol = ipopt(nlp; print_level=5, mu_strategy="adaptive", tol=1e-8, sb="yes")

This is Ipopt version 3.14.16, running with linear solver MUMPS 5.7.3.

Number of nonzeros in equality constraint Jacobian...:     3005
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      251

Total number of variables............................:     1004
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:      755
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-01 1.10e+00 1.99e-14   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -5.0000000e-03 7.36e-02 2.66e-15 -11.0 6.08e+00    -  1.00e+00 1.00e+00h  1
   2  6.0003829e+00 8.88e-16 1.78e-15 -11.0 6.01e+00    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 2

                                   (scaled)                 (unscaled)
Objective...............:   6.0003828724303254e+00    6.0003828724303254e+00
Dual infeasibility......:   1.7763568394002505e-15    1.7763568394002505e-15
Constraint violation....:   8.8817841970012523e-16    8.8817841970012523e-16
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   1.7763568394002505e-15    1.7763568394002505e-15


Number of objective function evaluations             = 3
Number of objective gradient evaluations             = 3
Number of equality constraint evaluations            = 3
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 3
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 2
Total seconds in IPOPT                               = 0.031

EXIT: Optimal Solution Found.

Then we can rebuild and plot an optimal control problem solution (note that the multipliers are optional, but the OCP costate will not be retrieved if the multipliers are not provided).

sol = OptimalControlSolution(docp; primal=nlp_sol.solution, dual=nlp_sol.multipliers)
plot(sol)

Change the NLP solver

Alternatively, we can use MadNLP.jl to solve anew the NLP problem:

using MadNLP

nlp_sol = madnlp(nlp)

"Execution stats: Optimal Solution Found (tol = 1.0e-08)."

Another possible NLP solver is Percival.jl.

using Percival

nlp_sol = percival(nlp; verbose=1)

"Execution stats: maximum elapsed time"

Initial guess

An initial guess, including warm start, can be passed to direct_transcription the same way as for solve.

docp, nlp = direct_transcription(ocp; init=sol)

It can also be changed after the transcription is done, with set_initial_guess.

set_initial_guess(docp, nlp, sol)