The solve function

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

solve(
    ocp::Model,
    description::Symbol...;
    display,
    kwargs...
) -> Solution{TimeGridModelType, TimesModelType, StateModelType, ControlModelType, VariableModelType, CostateModelType, Float64, DualModelType, CTModels.SolverInfos{Dict{Symbol, Any}}, ModelType} where {TimeGridModelType<:CTModels.TimeGridModel, TimesModelType<:CTModels.TimesModel, StateModelType<:Union{CTModels.StateModelSolution{TS} where TS<:CTModels.var"#114#136", CTModels.StateModelSolution{TS} where TS<:CTModels.var"#115#137"}, ControlModelType<:Union{CTModels.ControlModelSolution{TS} where TS<:CTModels.var"#116#138", CTModels.ControlModelSolution{TS} where TS<:CTModels.var"#117#139"}, VariableModelType<:Union{CTModels.VariableModelSolution{Vector{Float64}}, CTModels.VariableModelSolution{Float64}}, CostateModelType<:Union{CTModels.var"#118#140", CTModels.var"#119#141"}, DualModelType<:(CTModels.DualModel{PC_Dual, Vector{Float64}, SC_LB_Dual, SC_UB_Dual, CC_LB_Dual, CC_UB_Dual, Vector{Float64}, Vector{Float64}} where {PC_Dual<:Union{CTModels.var"#121#143"{CTModels.var"#120#142"}, CTModels.var"#122#144"{CTModels.var"#120#142"}}, SC_LB_Dual<:Union{CTModels.var"#124#146"{CTModels.var"#123#145"}, CTModels.var"#125#147"{CTModels.var"#123#145"}}, SC_UB_Dual<:Union{CTModels.var"#127#149"{CTModels.var"#126#148"}, CTModels.var"#128#150"{CTModels.var"#126#148"}}, CC_LB_Dual<:Union{CTModels.var"#130#152"{CTModels.var"#129#151"}, CTModels.var"#131#153"{CTModels.var"#129#151"}}, CC_UB_Dual<:Union{CTModels.var"#133#155"{CTModels.var"#132#154"}, CTModels.var"#134#156"{CTModels.var"#132#154"}}}), ModelType<:Model}

julia> available_methods()

julia> sol = solve(ocp)

julia> sol = solve(ocp, :direct, :adnlp, :ipopt)

julia> sol = solve(ocp, :direct)

julia> sol = solve(ocp, :direct, :ipopt, max_iter=100)

julia> sol = solve(ocp, init=(state=[-0.5, 0.2], control=0.5))

Basic usage

Let 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]
    ẋ(t)  == [v(t), u(t)]
    0.5∫( u(t)^2 ) → min
end

We can now solve the problem:

using NLPModelsIpopt
solve(ocp)

▫ This is OptimalControl version v1.1.0 running with: direct, adnlp, ipopt.

▫ The optimal control problem is solved with CTDirect version v0.15.2.

   ┌─ The NLP is modelled with ADNLPModels and solved with NLPModelsIpopt.
   │
   ├─ Number of time steps⋅: 250
   └─ Discretisation scheme: trapeze

▫ This is Ipopt version 3.14.17, running with linear solver MUMPS 5.8.0.

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 3.11e-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.0003828724303228e+00    6.0003828724303228e+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                               = 1.869

EXIT: Optimal Solution Found.

Note that we must import NLPModelsIpopt.jl before calling solve. This is because the default method uses a direct approach, which transforms the optimal control problem into a nonlinear program (NLP) of the form:

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

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

Resolution methods and algorithms

OptimalControl.jl offers a list of methods. To get it, simply call available_methods.

available_methods()

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

Each line is a method, with priority going from top to bottom. This means that

solve(ocp)

is equivalent to

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

1. The first symbol refers to the general class of method. The only possible value is:
   • :direct: currently, only the so-called direct approach is implemented. Direct methods discretise the original optimal control problem and solve the resulting NLP. In this case, the main solve method redirects to CTDirect.solve.
2. The second symbol refers to the NLP modeler. The possible values are:
   • :adnlp: the NLP problem is modeled by a ADNLPModels.ADNLPModel. It provides automatic differentiation (AD)-based models that follow the NLPModels.jl API.
   • :exa: the NLP problem is modeled by a ExaModels.ExaModel. It provides automatic differentiation and SIMD abstraction.
3. The third symbol specifies the NLP solver. Possible values are:
   • :ipopt: calls NLPModelsIpopt.ipopt to solve the NLP problem.
   • :madnlp: creates a MadNLP.MadNLPSolver instance from the NLP problem and solve it. MadNLP.jl is an open-source solver in Julia implementing a filter line-search interior-point algorithm like Ipopt.
   • :knitro: uses the Knitro solver (license required).

For instance, let us try MadNLP solver with ExaModel modeller.

using MadNLP

ocp = @def begin
    t ∈ [ t0, tf ], time
    x = (q, v) ∈ R², state
    u ∈ R, control
    x(t0) == x0
    x(tf) == [0, 0]
    ∂(q)(t) == v(t)
    ∂(v)(t) == u(t)
    0.5∫( u(t)^2 ) → min
end

solve(ocp, :exa, :madnlp)

▫ This is OptimalControl version v1.1.0 running with: direct, exa, madnlp.

▫ The optimal control problem is solved with CTDirect version v0.15.2.

   ┌─ The NLP is modelled with ExaModels and solved with MadNLP.
   │
   ├─ Number of time steps⋅: 250
   └─ Discretisation scheme: trapeze

▫ This is MadNLP version v0.8.7, running with umfpack

Number of nonzeros in constraint Jacobian............:     2004
Number of nonzeros in Lagrangian Hessian.............:     1751

Total number of variables............................:      753
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:      504
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  5.0000000e-03 1.10e+00 5.55e-17  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  6.0003829e+00 1.85e-15 2.15e-13  -1.0 6.08e+00    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:   6.0003828724303929e+00    6.0003828724303929e+00
Dual infeasibility......:   2.1493917756743031e-13    2.1493917756743031e-13
Constraint violation....:   1.8492152253912764e-15    1.8492152253912764e-15
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   2.1493917756743031e-13    2.1493917756743031e-13

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.)  = 13.004
Total wall-clock secs in linear solver                      =  0.307
Total wall-clock secs in NLP function evaluations           =  0.000
Total wall-clock secs                                       = 13.311

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

Note that you can provide a partial description. If multiple full descriptions contain it, priority is given to the first one in the list. Hence, all of the following calls are 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)

Direct method

The main options for the direct method, with their [default] values, are:

• display ([true], false): setting display = false disables output.
• init: information for the initial guess. It can be given as numerical values, functions, or an existing solution. See how to set an initial guess.
• grid_size ([250]): number of time steps in the (uniform) time discretization grid. More precisely, if N = grid_size and the initial and final times are t0 and tf, then the step length Δt = (tf - t0) / N.
• time_grid ([nothing]): explicit time grid (can be non-uniform). If time_grid = nothing, a uniform grid of length grid_size is used.
• disc_method ([:trapeze], :midpoint, :euler, :euler_implicit, :gauss_legendre_2, :gauss_legendre_3): the discretisation scheme to transform the dynamics into nonlinear equations. See the discretization method tutorial for more details.
• adnlp_backend ([:optimized], :manual, :default): backend used for automatic differentiation to create the ADNLPModels.ADNLPModel.

For advanced usage, see:

• discrete continuation tutorial,
• NLP manipulation tutorial.

NLP solvers specific options

In addition to these options, all remaining keyword arguments passed to solve will be transmitted to the NLP solver used.

Please check the list of Ipopt options and the NLPModelsIpopt.jl documentation.

sol = solve(ocp; max_iter=0, tol=1e-6, display=false)
iterations(sol)

0

Similarly, please check the MadNLP.jl documentation and the list of MadNLP.jl options.

sol = solve(ocp, :madnlp; max_iter=0, tol=1e-6, display=false)
iterations(sol)

0