Getting Started

This section demonstrates the basic usage of pockit by solving the optimal time double integrator problem. It assumes that pockit has already been installed following the instructions in installation guide. For a detailed introduction to multi-phase optimal control problems and pockit components, please refer to the User Guide and API documentation.

Introduction

The optimization variables of an optimal control problem are infinite-dimensional functions, and the constraints on the function are infinite-dimensional either. To solve such problems numerically, one approach is to give a finite-dimensional approximation of the problem. This typically involves two steps:

• Finite-dimensional discretization and approximation of the state and control variables (i.e., functions over the time interval).
• Further finite-dimensional approximation of the objective function and constraints.

Then, the discretized problems can be solved using nonlinear programming (NLP) solvers.

Pockit provides the following main features:

• Input the optimal control problem and discretization parameters to generate the discrete NLP problem
• Convenient interface for interacting with NLP solvers, including providing initial guesses and extracting computation results
• Discretized solutions are only approximations of the original continuous problem's solution. Pockit provides methods to check whether the current solution is within the allowable error tolerance and can automatically adjust the discretization parameters according to the target error range. Results from previous iterations can be automatically interpolated as initial guesses for the new grid points, allowing for solving the continuous problem iteratively.

Problem Description

In this example, we will use pockit to solve the optimal time double integrator problem.

Consider the following optimal control problem, which has 2 state variables \(x(t)\), \(v(t)\), 1 control variable \(u(t)\), and a free final time \(t_f\):

$$ \begin{array}{cl} \min & \displaystyle t_f = \int_0^{t_f} 1 \,\mathrm{d}t\\ \text{s.t.} & \dot{x}(t) = v(t), \quad x(0) = 1, \quad x(t_f) = 0\\ & \dot{v}(t) = u(t), \quad v(0) = 0, \quad v(t_f) = 0\\ & -1 \leq u(t) \leq 1, \quad \forall\, t \in [0, t_f] \end{array} $$

This is a simple bang-bang control problem, where the optimal control is a discontinuous function. The optimal solution is \( \hat t_f = 2 \), with the optimal control being:

$$ \hat u(t) = \begin{cases}-1,\quad t<1\\ \text{any value}\in[-1,1],\quad t=1\\ 1,\quad t>1\end{cases} $$

The method for acquiring the analytical solution can be found in optimal control textbooks.

Discretizing the Problem

Pockit is designed to solve multi-phase optimal control problems, which consist of a single system object and one or more phase objects. However, it is also efficient for common single-phase problems. By inputting the continuous problem into pockit, it will automatically generate the discretized problem.

Start by importing the required modules.

from pockit.radau import System, linear_guess
from pockit.optimizer import ipopt
# If using Scipy backend:
# from openoc.optimizer import scipy

Since the problem is a single-phase problem, only one phase object is needed.

S = System(0) # Parameter 0 means no system parameters
P = S.new_phase(['x', 'v'], ['u']) # x, v, u are the variable names

Pockit automatically generates Sympy symbols for the state and control variables, which can be accessed and extracted:

x, v = P.x
u, = P.u

These symbols will be used to define all the functional relationships in the problem. First, define the dynamics of the state variables:

P.set_dynamics([v, u])
# The argument is a list of state variable dynamics expressions. The length must equal the number of state variables.

In multi-phase problems, the objective is defined at the system level. In each phase, we define the integrals that can be accessed by the system and used to form the objective function. In this problem, only one integral is of interest, the final time \(t_f\):

P.set_integral([1]) 
# The argument is a list of expressions for the integrals of interest.

Next, define the phase constraints (also called path constraints), which must hold for any time within the phase time interval:

P.set_phase_constraint([u], [-1], [1], [True])
# The first argument is a list of expressions for the phase constraint functions.
# The second argument is a list of lower bounds for the constraint functions.
# The third argument is a list of upper bounds for the constraint functions.
# The fourth argument is a boolean list, indicating whether the constraint functions are bang-bang controls.
# The length of these four arguments must be the same.

Now, define the boundary conditions for the phase:

P.set_boundary_condition([1, 0], [0, 0], 0, None)
# The first and second arguments are lists for the initial and terminal boundary conditions for the state variables.
# The third and fourth arguments are the start and end times of the phase.
# The free parameter is set to None.

Finally, specify the discretization parameters for the phase:

P.set_discretization(3, 5)
# The first parameter is the number of subintervals,
# and the second parameter is the number of interpolation nodes per subinterval.

At this point, all phase-level relationships have been defined. Now, set the phase in the system object:

S.set_phase([P])

Next, define the objective function for the system:

S.set_objective(P.I[0])
# The integrals defined in the phase automatically generate Sympy symbols, which can be accessed through P.I.

We have now completed the problem setup, and pockit has automatically generated the discretized problem. To summarize, we defined the problem through the following steps:

1. Create system and phase objects.
2. Set dynamics, integrals, phase constraints, boundary conditions, and discretization parameters for each phase (dynamics, boundary conditions, and discretization are mandatory).
3. Set phases, objectives, and system constraints for the system object (phases and objectives are mandatory).

Providing Initial Guess

To pass the discretized problem to the NLP solver, we need to provide an initial guess of the solution. The quality of the initial guess is crucial for the success of the optimization. Pockit uses Variable objects to handle the discretized optimization variables for each phase. First, we let pockit generate a guess based on linear interpolation of the boundary conditions, which will serve as the starting point:

v = linear_guess(P)

Linear interpolations are rough but provide a starting point that is compatible with the phase's discretization scheme. We can manually adjust the guess to make it more accurate. First, set the final time:

v.t_f = 5

Then, modify the guesses for the state and control variables:

v.x[1] = np.zeros_like(v.t_x)
v.u[0] = np.array([-1. if t < 2 else 1. for t in v.t_u])

The time interpolation nodes for the state variables are stored in the attribute v.t_x, and the time interpolation nodes for the control variables are stored in v.t_u. Note:

• The interpolation nodes are not uniformly distributed over the interval.
• The state variables and control variables may have different interpolation nodes.

Solving the Discretized Problem

Finally, we call the NLP solver to solve the discretized problem:

v, info = ipopt.solve(S, v)
print(info['status_msg'].decode())
# Or use the Scipy backend:
# v, info = scipy.solve(S, v)
# print(info['success'])

Plotting the Numerical Solution

Now, we can access the solution as follows:

print(v.t_f)
print(v.x[0])
print(v.x[1])
print(v.u[0])

However, these numbers are hard to interpret. As mentioned earlier, the interpolation nodes are stored in v.t_x and v.t_u. The correct way to plot the variables is to plot them against the interpolation nodes, as shown below:

import matplotlib.pyplot as plt
fig, axs = plt.subplots(3, 1, sharex=True)
axs[0].plot(v.t_x, v.x[0], marker='+')
axs[1].plot(v.t_x, v.x[1], marker='+')
axs[2].plot(v.t_u, v.u[0], marker='+')
plt.show()

The resulting plot looks like this:

Adjusting the Grid

Unfortunately, as seen in the plot, the solution is not accurate. As mentioned earlier, the optimal control is bang-bang control. Pockit uses piecewise polynomial interpolation to approximate the functions, so discontinuities are only allowed at mesh points. In this case, the theoretical analysis shows that the control switching time is \( \hat t_s = 1 \), which happens to be exactly in the middle of the middle subinterval. For general problems, we cannot determine the switching time analytically. Pockit provides a method to approximate the switching time through multiple iterations of numerical refinement:

max_iter = 10
for it in range(max_iter):
    if S.check(v):
        break
    v = S.refine(v)
    v, info = ipopt.solve(S, v)
    assert info["status"] == 0
    print(it, info["obj_val"])
    # Or use Scipy backend:
    # v, info = scipy.solve(S, v)
    # assert info["success"]
    # print(info["fun"])

In the code snippet, we first set the maximum number of iterations to 10. We then use System.check() to check if the solution is accurate enough. The allowable error range is adjustable, but for simplicity, we use the default value here. If the solution passes the test, the loop breaks. Otherwise, we use System.refine() to adjust the discretization parameters. This method adjusts the number of mesh points and interpolation nodes in each subinterval based on the current optimization variables. After adjusting the grid, the method returns a new Variable object that is compatible with the new discretization parameters, and its values are interpolated from the old solution. We then use this new variable as the initial guess to solve the problem. This process is repeated until the solution passes the test or the maximum iterations are reached. The final solution looks like this:

The solution is now very accurate, with the objective function value \( t_f = 2.000000044988834 \), which is very close to the theoretical solution. Moreover, the control exhibits bang-bang behavior (the + markers in the plot represent the interpolation nodes).

Summary

This brief tutorial demonstrated how to use pockit to solve a simple optimal control problem. We discretized the problem by inputting the problem definition into pockit. Then, we generated an initial guess for the solution and solved the problem. Finally, we improved the accuracy of the solution by iterating and adjusting the mesh. We also showed how to access the solution and plot the results.

For further details, refer to the User Guide and API documentation.

Complete Code

import numpy as np
import matplotlib.pyplot as plt
from pockit.optimizer import ipopt
from pockit.radau import System, linear_guess

S = System(0)
P = S.new_phase(['x', 'v'], ['u'])
x, v = P.x
u, = P.u
P.set_dynamics([v, u])
P.set_integral([1])
P.set_phase_constraint([u], [-1], [1], [True])
P.set_boundary_condition([1, 0], [0, 0], 0, None)
P.set_discretization(3, 5)
S.set_phase([P])
S.set_objective(P.I[0])

v = linear_guess(P)
v.t_f = 5
v.x[1] = np.zeros_like(v.t_x)
v.u[0] = np.array([-1. if t < 2 else 1. for t in v.t_u])

v, info = ipopt.solve(S, v)
print(info['status_msg'].decode())

max_iter = 10
for it in range(max_iter):
    if S.check(v):
        break
    v = S.refine(v)
    v, info = ipopt.solve(S, v)
    assert info["status"] == 0
    print(it, info["obj_val"])

fig, axs = plt.subplots(3, 1, sharex=True)
axs[0].plot(v.t_x, v.x[0], marker='+')
axs[1].plot(v.t_x, v.x[1], marker='+')
axs[2].plot(v.t_u, v.u[0], marker='+')
plt.show()