The solve function

In this tutorial, we explain the solve function from OptimalControl.jl package.

Basic usage

Les us define a basic optimal control problem.

using OptimalControl

t0 = 0
tf = 1
x0 = [-1, 0]

ocp = @def begin

    t ∈ [ t0, tf ], time
    x = (q, v) ∈ R², state
    u ∈ R, control

    x(t0) == x0
    x(tf) == [ 0, 0 ]

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

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

end

Let us try to solve the problem:

julia> solve(ocp)
ERROR: ExtensionError. Please make: julia> using NLPModelsIpopt

As you can see, an error occured since we need the package NLPModelsIpopt.jl.

Actually, the default solving method is what we call a direct method. In a direct method, the optimal control problem is transcribed to a nonlinear optimization problem (NLP) of the form

\[\text{minimize}\quad F(y), \quad\text{subject to the constraints}\quad g(y)=0, \quad h(y)\le 0. \]

OptimalControl.jl package makes the transcription but it needs a package to model the NLP problem and another one to solve it. NLPModelsIpopt.jl package provides an interface to the well-known solver Ipopt that can be used to solve general nonlinear programming problems.

using NLPModelsIpopt

solve(ocp)

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.

Solvers

OptimalControl.jl offers a list of methods to solve your optimal control problem. To get the list of methods, simply call available_methods.

available_methods()

(:direct, :adnlp, :ipopt)
(:direct, :adnlp, :madnlp)

Each line is a method. The priority is given from top to bottom. This means that

solve(ocp)

is equivalent to

solve(ocp, :direct, :adnlp, :ipopt)

Let us detail the three symbols. As you can see, there are only direct methods. The symbol :adnlp is for the choice of modeler. As said before, the NLP problem needs to be modelled in Julia code. We use ADNLPModels.jl which provides automatic differentiation (AD)-based model implementations that conform to the NLPModels.jl API.

The last symbol is what distinguish the two available methods. It corresponds to the NLP solver. By default, we use the NLPModelsIpopt.jl package but you can choose MadNLP.jl which is an open-source nonlinear programming solver, purely implemented in Julia and which implements a filter line-search interior-point algorithm, as used in Ipopt.

Note that you are not compelled to provide the full description of the method to use it. A partial description, if not ambiguous, will work. Just remember that if several full descriptions contain the partial one, then, the priority is given to first one in the list. Hence, these calls are all equivalent:

solve(ocp)
solve(ocp, :direct                )
solve(ocp,          :adnlp        )
solve(ocp,                  :ipopt)
solve(ocp, :direct, :adnlp        )
solve(ocp, :direct,         :ipopt)
solve(ocp, :direct, :adnlp, :ipopt)

Let us try MadNLP.jl.

using MadNLP

solve(ocp, :direct, :adnlp, :madnlp; display=true)

The options set during resolve may not have an effect
This is MadNLP version v0.8.4, running with umfpack

Number of nonzeros in constraint Jacobian............:     3005
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 2.45e-13  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -5.0000000e-03 7.36e-02 3.46e-12  -1.0 6.08e+00    -  1.00e+00 1.00e+00h  1
   2  6.0003829e+00 1.78e-15 2.66e-15  -2.5 6.01e+00    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 2

                                   (scaled)                 (unscaled)
Objective...............:   6.0003828724303387e+00    6.0003828724303387e+00
Dual infeasibility......:   2.6645352591003757e-15    2.6645352591003757e-15
Constraint violation....:   1.7763568394002505e-15    1.7763568394002505e-15
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   2.6645352591003757e-15    2.6645352591003757e-15

Number of objective function evaluations             = 3
Number of objective gradient evaluations             = 3
Number of constraint evaluations                     = 3
Number of constraint Jacobian evaluations            = 3
Number of Lagrangian Hessian evaluations             = 2
Total wall-clock secs in solver (w/o fun. eval./lin. alg.)  =  1.690
Total wall-clock secs in linear solver                      =  0.004
Total wall-clock secs in NLP function evaluations           =  0.007
Total wall-clock secs                                       =  1.702

EXIT: Optimal Solution Found (tol = 1.0e-08).

Again, there are several equivalent manners to use MadNLP.jl.

solve(ocp,                  :madnlp)
solve(ocp, :direct,         :madnlp)
solve(ocp,          :adnlp, :madnlp)
solve(ocp, :direct, :adnlp, :madnlp)

Options

Direct method

The options common to all the direct methods may be seen directly in the direct_solve keywords. There are init, grid_size and time_grid.

The init option can be used to set an initial guess for the solver. See the initial guess tutorial.
The grid_size option corresponds to the size of the (uniform) time discretization grid. More precisely, it is the number of steps, that is if N = grid_size and if the initial and final times are denoted respectively t0 and tf, then we have:

Δt = (tf - t0) / N

The time_grid option is the grid of times: t0, t1, ..., tf. If the initial and/or the final times are free, then you can provide a normalised grid between 0 and 1. Note that you can set either grid_size or time_grid but not both.

sol = solve(ocp; grid_size=10, display=false)
time_grid(sol)

11-element Vector{Float64}:
 0.0
 0.1
 0.2
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0.9
 1.0

Or with MadNLP.jl:

sol = solve(ocp, :madnlp; grid_size=10, display=false)
time_grid(sol)

11-element Vector{Float64}:
 0.0
 0.1
 0.2
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0.9
 1.0

NLPModelsIpopt

You can provide any option of NLPModelsIpopt.jl or Ipopt with a pair keyword=value. Please check the list of Ipopt options and the NLPModelsIpopt.jl documentation.

solve(ocp, :ipopt; max_iter=0)

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

Number of Iterations....: 0

                                   (scaled)                 (unscaled)
Objective...............:   1.0000000000000001e-01    1.0000000000000001e-01
Dual infeasibility......:   1.9885771589406212e-14    1.9885771589406212e-14
Constraint violation....:   1.1000000000000001e+00    1.1000000000000001e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   1.1000000000000001e+00    1.1000000000000001e+00


Number of objective function evaluations             = 1
Number of objective gradient evaluations             = 1
Number of equality constraint evaluations            = 1
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             = 0
Total seconds in IPOPT                               = 0.075

EXIT: Maximum Number of Iterations Exceeded.

MadNLP

If you use the MadNLP.jl solver, then you can provide any option of it. Please check the MadNLP.jl documentation and the list of MadNLP.jl options.

solve(ocp, :madnlp; max_iter=1, display=true)

The options set during resolve may not have an effect
This is MadNLP version v0.8.4, running with umfpack

Number of nonzeros in constraint Jacobian............:     3005
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 2.45e-13  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -5.0000000e-03 7.36e-02 3.46e-12  -1.0 6.08e+00    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:  -4.9999999999990191e-03   -4.9999999999990191e-03
Dual infeasibility......:   3.4621194799910882e-12    3.4621194799910882e-12
Constraint violation....:   7.3553711558341242e-02    7.3553711558341242e-02
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   7.3553711558341242e-02    7.3553711558341242e-02

Number of objective function evaluations             = 2
Number of objective gradient evaluations             = 2
Number of constraint evaluations                     = 2
Number of constraint Jacobian evaluations            = 2
Number of Lagrangian Hessian evaluations             = 1
Total wall-clock secs in solver (w/o fun. eval./lin. alg.)  =  0.004
Total wall-clock secs in linear solver                      =  0.003
Total wall-clock secs in NLP function evaluations           =  0.005
Total wall-clock secs                                       =  0.011

EXIT: Maximum Number of Iterations Exceeded.