The solve function

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

CommonSolve.solveMethod
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}

Solve the optimal control problem ocp by the method given by the (optional) description. The get the list of available methods:

julia> available_methods()

The higher in the list, the higher is the priority. The keyword arguments are specific to the chosen method and represent the options of the solver.

Arguments

  • ocp::OptimalControlModel: the optimal control problem to solve.
  • description::Symbol...: the description of the method used to solve the problem.
  • kwargs...: the options of the method.

Examples

The simplest way to solve the optimal control problem is to call the function without any argument.

julia> sol = solve(ocp)

The method description is a list of Symbols. The default is

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

You can provide a partial description, the function will find the best match.

julia> sol = solve(ocp, :direct)
Note

See the resolution methods section for more details about the available methods.

The keyword arguments are specific to the chosen method and correspond to the options of the different solvers. For example, the keyword max_iter is an Ipopt option that may be used to set the maximum number of iterations.

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

See the direct method section for more details about associated options. These options also detailed in the CTDirect.solve documentation. This main solve method redirects to CTDirect.solve when the :direct Symbol is given in the description. See also the NLP solvers section for more details about Ipopt or MadNLP options.

To help the solve converge, an initial guess can be provided within the keyword init. You can provide the initial guess for the state, control, and variable.

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

See how to set an initial guess for more details.

source

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.1.

   ┌─ 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 2.73e-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 1.78e-15 1.78e-15 -11.0 6.01e+00    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 2

                                   (scaled)                 (unscaled)
Objective...............:   6.0003828724303263e+00    6.0003828724303263e+00
Dual infeasibility......:   1.7763568394002505e-15    1.7763568394002505e-15
Constraint violation....:   1.7763568394002505e-15    1.7763568394002505e-15
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.840

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. \]

Caveat

Calling solve without loading a NLP solver package first will notify the user:

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:
  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).
Warning

The dynamics must be defined coordinatewise to use ExaModels.jl (:exa).

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.1.

   ┌─ 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.144
Total wall-clock secs in linear solver                      =  0.311
Total wall-clock secs in NLP function evaluations           =  0.000
Total wall-clock secs                                       = 13.455

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:

Note

The main solve method from OptimalControl.jl simply redirects to CTDirect.solve in that case.

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