Problem Definition

This section outlines the general structure of an optimal control problem and its implementation in YAPSS. See the following sections for details of how to define callback function, set variable constraint bounds, define the initial guess, select options for evaluating derivatives, select solver options, configure the Ipopt NLP solver source, and access the solution.

General Formulation

YAPSS implements in Python a pseudospectral method for solving optimal control problems, similar to the MATLAB algorithm GPOPS-II [1] described by Patterson and Rao. YAPSS generalizes the GPOPS-II algorithm by allowing Legendre-Gauss (LG) and Legendre-Gauss-Loabatto (LGL) collocation points in addition to the Legendre-Gauss-Radau (LGR) points used by GPOPS-II. In this formulation, the optimal control problem is defined over multiple phases, each with its own dynamics, path constraints, and integrals. In addition, the problem may depend on a vector of static parameters that are to be optimized over as well.

The general multistage optimal control problem is defined as follows: There are \(n_{p}\) phases, each with a state vector \(x^{(p)}\) with dimension \(n_{x}^{(p)}\), and a control vector \(u^{(p)}\) with dimension \(n_{u}^{(p)}\). In addition, there is a static parameter vector \(s\) with dimension \(n_{s}\) that applies to all phases. The dynamics of the problem for each phase are then given by

\[\dot{x}^{(p)} = f^{(p)}(x^{(p)}, u^{(p)}, t, s), \quad p=0, \ldots, n_{p}-1\]

over the interval \(t \in [t_{0}^{(p)}, t_{f}^{(p)}]\). The trajectory of the system over each phase is subject to the path constraints

\[h_{\text{min}}^{(p)} \leq h^{(p)} ( x^{(p)}, u^{(p)}, t, s ) \leq h_{\text{max}}^{(p)}, \quad p=0, \ldots, n_{p}-1\]

where \(h^{(p)}\) is a vector-valued function with dimension \(n_{h}^{(p)}\). In addition, each phase may have integrals associated with it of the form

\[q^{(p)} = \int_{t_{0}^{(p)}}^{t_{f}^{(p)}} g^{(p)}( x^{(p)}, u^{(p)}, t, s ) \,dt, \quad p=0, \ldots, n_{p}-1\]

where \(g^{(p)}\) is a vector-valued function with dimension \(n_{q}^{(p)}\). The integrals may appear as a term in the cost function (a Lagrangian term) or as a perimetric constraint.

The cost to be minimized is given by a function of all the discrete variables in problem:

\[\begin{split}\begin{aligned} J=\phi\Big[ & x^{(0)} (t_{0}^{(0)}),\ldots,x^{(n_{p}-1)}(t_{0}^{(n_{p}-1)}), t_{0}^{(0)},\ldots,t_{0}^{(n_{p}-1)}, \\ & x^{(0)}(t_{f}^{(0)}),\ldots,x^{(n_{p}-1)}(t_{f}^{(n_{p}-1)}), t_{f}^{(0)},\ldots,t_{f}^{(n_{p}-1)},{q}^{(0)},\ldots,{q}^{(n_{p}-1)},{s}\Big] \end{aligned}\end{split}\]

subject to additional constraints on the discrete variables,

\[\begin{split}\begin{aligned} d_\text{min} \le d \Big[ & x^{(0)} (t_{0}^{(0)}),\ldots,x^{(n_{p}-1)}(t_{0}^{(n_{p}-1)}), t_{0}^{(0)},\ldots,t_{0}^{(n_{p}-1)}, \\ & x^{(0)}(t_{f}^{(0)}),\ldots,x^{(n_{p}-1)}(t_{f}^{(n_{p}-1)}), t_{f}^{(0)},\ldots,t_{f}^{(n_{p}-1)},{q}^{(0)},\ldots,{q}^{(n_{p}-1)},{s} \Big] \le d_\text{max} \end{aligned}\end{split}\]

where the function \(d\) is a vector-valued with dimension \(n_{d}\).

In addition, upper and lower bounds may be specified on all the decision variables:

  • The state \(x^{(p)}\) and the control \(u^{(p)}\) vectors

  • The initial and final state vectors \(x^{(p)}(t_{0}^{(p)})\) and \(x^{(p)}(t_{f}^{(p)})\)

  • The initial and final times \(t_{0}^{(p)}\) and \(t_{f}^{(p)}\)

  • The parameter vector \(s\)

  • The integrals \(q^{(p)}\)

(In this formulation, the integrals are treated as decision variables, subject to the constraint that they are the integrals of the integrand over the phase.)

Problem Instantiation

Given a formulation of the problem as described above, the problem can be implemented in YAPSS and solved. In this section, we describe the instantiation of a YAPSS problem object.

Consider for example, the Goddard problem, a classic optimal control problem to maximize the altitude of a sounding rocket launched vertically from the surface of the Earth, taking into account the forces of gravity, drag, and thrust. (See the JupyterLab notebooks for the one phase Goddard problem and the three phase Goddard problem with a singular arc.) For the three phase problem, there are three states (altitude, velocity, and mass), one control (thrust), and one path constraint (the singular arc constraint). In addition, there are eight discrete constraints that ensure the continuity of the state variables across phases, and that the final time of each phase is equal to the initial time of the next phase.

The problem is instantiated as an instance of the Problem class. To initialize, the user specifies the name of the problem; the relevant dimensions of the state, control, and path constraint for each phase; and the dimensions of the static parameters and discrete constraints. In this case, the problem is instantiated as follows:

>>> import yapss.numpy as np
>>> from yapss import Problem, Solution
>>>
>>> problem = Problem(
...     name="Goddard Rocket Problem with Singular Arc",
...     nx=[3, 3, 3],  # Number of states in each phase
...     nu=[1, 1, 1],  # Number of controls in each phase
...     nh=[0, 1, 0],  # Number of path constraints in each phase
...     nd=8,          # Number of discrete constraints
... )

Note that all arguments are keyword arguments, and the name and nx arguments are required. The number of phases is determined by the length of the nx argument. If an optional argument is not provided, the assumed value is either zero or a tuple of zeros, as appropriate.

The string representation of the problem object provides a summary of the problem:

>>> print(problem)
Problem(
    name='Goddard Rocket Problem with Singular Arc',
    nx=(3, 3, 3),
    nu=(1, 1, 1),
    nq=(0, 0, 0),
    nh=(0, 1, 0),
    nd=8,
    ns=0
)

Other sections of this reference describe the remaining steps required to so solve an optimal control problem using YAPSS:

Problem Class Reference

class Problem(*, name: str, nx: Sequence[int], nu: Sequence[int] | None = None, nq: Sequence[int] | None = None, nh: Sequence[int] | None = None, ns: int | None = 0, nd: int | None = 0)[source]

Instances of the Problem class define the optimal control problem.

Parameters:
namestr

Name of the optimal control problem

nxSequence[int]

Number of states in each phase

nuOptional[Sequence[int]]

Number of controls in each phase

nqSequence[int], optional

Number of integrals in each phase

nhSequence[int], optional

Number of path constraints in each phase

nsint, optional

Number of parameters in problem

ndint, optional

Number of discrete constraints in problem

Attributes:
namestr

The name of the problem.

npint

The number of phases in the problem.

nxtuple[int, …]

The number of state variables in each phase.

nutuple[int, …]

The number of control variables in each phase.

nqtuple[int, …]

The number of integrals in each phase.

nhtuple[int, …]

The number of path constraints in each phase.

ndint

The number of discrete constraints in the problem.

nsint

The number of parameters in the problem.

auxdataAuxdata

The auxiliary data for the problem.

boundsBounds

The bounds object structure for the problem.

derivativesDerivatives

The derivative options for the problem. Attributes: method, order.

functionsUserFunctions

The user-defined functions for the problem. Attributes: objective, continuous, discrete, objective_gradient, objective_hessian, continuous_jacobian, continuous_hessian, discrete_jacobian, discrete_hessian.

guessGuess

The initial guess for the problem. Attributes: parameter, phase[p].time, phase[p].state, phase[p].control, phase[p].integral.

ipopt_optionsIpoptOptions

The user-selected Ipopt options for the problem. To select a particular Ipopt option, use the ipopt_options attribute of the Problem class. For example, for a Problem instance problem, set the tol option to 1e-6 by setting problem.ipopt_options.tol = 1e-6.

scaleScale

The scaling data structure for the problem.

meshMesh

The mesh data structure for the problem.

spectral_method{“lg”, “lgr”, “lgl”}

The type of interpolation used for the problem.

solve() Solution[source]

Solve the optimal control problem.

Returns:
solutionSolution

The solution to the optimal control problem.

validate() None[source]

Validate the optimal control problem input.

Raises:
ValueError

If the problem is invalid.

References