Setting Up Problem in Pockit

Introduction

After establishing a multi-phase optimal control model for the problem, you only need to input the expressions into pockit to automatically generate a discrete optimization problem. For the definition of a multi-phase optimal control problem, please refer to Multi-Phase Optimal Control Problem.

Pockit is entirely based on constructing expressions with Sympy. When objects are created, the corresponding Sympy symbols are automatically generated. Use these symbols to construct the required expressions and input them into pockit. This allows you to conveniently manipulate expressions using Sympy library functions, and they can be displayed as nicely formatted formulas in environments like Jupyter Notebook.

Creating Objects

Pockit has two types of interpolation nodes: Radau and Lobatto, both having identical programming interfaces. The import methods are:

from pockit.radau import System, Phase
from pockit.labotto import System, Phase
# Choose one

For the choice between Radau and Lobatto nodes, refer to Discretization and Interpolation.

To create a system object, use the System class constructor

S = System(n_s)
# n_s is the number of system variables, or
S = System(['s_0', 's_1', ..., 's_{n_s - 1}'])
# The list contains the names of the system variables

Use the System.new_phase method to create phase objects

P = S.new_phase(n_x, n_u)
# n_x, n_u are the numbers of state and control variables, or
P = S.new_phase(['name_x_0', ..., 'name_x_{n_x - 1}'],
                ['name_u_0', ..., 'name_u_{n_u - 1}'])
# The list contains the names of the state and control variables

Setting Function Relationships

After creating the objects, the corresponding Sympy symbols are automatically generated. We can use these symbols to construct the function relationships of the problem and input the expressions to build the problem. Note: set_* methods may be slow because they need to calculate symbolic derivatives and compile them into Numba functions. It is best to execute each of them in a separate Jupyter Notebook cell. If any of them are modified (e.g., changing the boundary conditions for the same system), you do not need to re-execute the other cells, just execute System.update() at the end to update the system object data (see below). A lot of time can be saved by doing this.

Accessing Variables

System variables can be accessed through the System.s and Phase.s properties, both yielding the same result. The Sympy symbol lists of state and control variables can be accessed through the Phase.x and Phase.u properties. Additionally, each Phase automatically generates an extra Sympy symbol P.t representing the time variable.

Setting Phases

We need to first configure the function relationships and discretization methods for all phases, and then integrate them at the system level. All input mathematical functions should be constructed as Sympy expressions composed of symbols from the respective properties.

Set the dynamic differential equations: Phase.set_dynamics

P.set_dynamics([expr_0, ..., expr_{n_x - 1}])
# expr_* are expressions composed of symbols from P.x, P.u, P.s, and P.t

Set the integrals: Phase.set_integral

P.set_integral([expr_0, ..., expr_{n_i - 1}])
# expr_* are expressions composed of symbols from P.x, P.u, P.s, and P.t

Set phase (path) constraints: Phase.set_phase_constraint

P.set_phase_constraint(
    [expr_0, ..., expr_{n_c - 1}],
    [lb_0, ..., lb_{n_c - 1}],
    [ub_0, ..., ub_{n_c - 1}],
    [b_0, ..., b_{n_c - 1}],
)
# expr_* are expressions composed of symbols from P.x, P.u, P.s, and P.t
# lb_*, ub_* are floating-point numbers (can be +/- inf), lb_i <= ub_i
# b_* are True/False, indicating whether the phase constraint is bang-bang control

Set boundary conditions: Phase.set_boundary_condition

P.set_boundary_condition(
    [x_0_0, ..., x_0_{n_x - 1}],
    [x_f_0, ..., x_f_{n_x - 1}],
    t_0,
    t_f,
)
# x_0_*, x_f_*, t_0, t_f are None or floating-point numbers or expressions composed of symbols from S.s
# Representing free, fixed to a constant, and fixed to a function respectively

Set discretization parameters: Phase.set_discretization

P.set_discretization(
  [mesh_0, ..., mesh_{n_i}],
  [num_points_0, ..., num_points_{n_i - 1}],
)
# mesh_* are the positions of the mesh points, which are monotonically increasing floating-point numbers
# (can also be simplified to an integer n_i, representing n_i uniform subintervals)
# num_points_* are the number of interpolation points in each subinterval
# (can also be simplified to an integer n_p, representing n_p interpolation points in each subinterval)

Setting System Functions

After setting the phase relationships, combine them at the system level to construct the final optimization problem.

Set phases: System.set_phase
```
S.set_phase(P_0, ..., P_{n_p - 1})
```
If the functions or discretization parameters of any phase are changed after setting the phases (through set_* methods), you need to call System.update() at the system level to update data in the system object.

Set the objective function: System.set_objective

S.set_objective(expr)
# expr is an expression composed of symbols from P_0.I, ..., P_{n_p - 1}.I and S.s

Set system (algebraic) constraints: System.set_system_constraint

S.set_system_constraint(
  [expr_0, ..., expr_{n_l - 1}],
  [lb_0, ..., lb_{n_l - 1}],
  [ub_0, ..., ub_{n_l - 1}])
)
# expr_* are expressions composed of symbols from P_1.I, ..., P_{n_p - 1}.I and S.s
# lb_*, ub_* are floating-point numbers (can be +/- inf), lb_i <= ub_i