快速开始

本节通过求解最优时间双积分问题，展示 pockit 的基本使用方法。本节假设已按照下载安装中的说明安装了 pockit。对于多段最优控制问题和 pockit 组件的详细介绍，请参阅用户指南和 API 文档。

简介

最优控制问题的优化对象为无穷维的函数，约束也为对函数的无穷维约束。为了实现问题的数值求解，一种方法为给出问题的有限维近似。主要包含如下两个方面：

对状态和控制变量（即时间区间上的函数）进行有限维离散和近似
进一步对目标函数、约束进行有限维近似

然而再利用非线性规划（NLP）求解器求解离散问题。

Pockit 主要提供如下功能：

输入最优控制问题和离散参数，生成连续问题
与 NLP 求解器便捷交互的接口，包括提供初值和提取计算结果
离散化解只是原始连续问题的解的近似。Pockit 提供方法检验当前解是否达到了误差容许范围，并可以根据目标误差范围自动调整离散网格。上一轮计算结果可以自动插值为适合新网格点的初值并进行下一轮迭代，从而实现问题的离散求解

问题描述

在这篇简介中，利用 pockit 求解最优时间双积分器问题。

考虑以下最优控制问题，其中有 2 个状态变量 $x(t), v(t)$， 1 个控制变量 $u(t)$ 和自由的结束时间 $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} $$

问题是一个最简单的 bang-bang 控制问题，即最优控制为间断函数。最优解为 $\hat t_f=2$，

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

解析解的计算方法可参阅最优控制教材。

生成离散问题

Pockit 为了更加通用，被设计用于解决多段最优控制问题，其由单个系统对象和一个或多个阶段对象构成。但其对常见的单段问题也同样保持高效。将连续问题输入至 pockit 即可自动生成编译离散问题。

首先导入模块。

from pockit.radau import System, linear_guess
from pockit.optimizer import ipopt
# 若使用 Scipy 后端：
# from openoc.optimizer import scipy

由于问题为单段问题，仅需一个阶段对象。

S = System(0) # 参数 0 表示没有系统参数
P = S.new_phase(['x', 'v'], ['u']) # x, v, u 为变量名

Pockit 自动为状态和控制变量生成 Sympy 符号，可以访问并提取

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

这些符号后续用于定义问题中的所有函数关系。我们首先定义系统的动力学

P.set_dynamics([v, u])
# 参数是状态变量动力学表达式列表，其长度必须等于状态变量的数量

在多阶段问题中，目标在系统级别定义。在每个阶段，我们定义可以被系统访问和用来构成目标函数的积分。在这个问题中，只有一个关注的积分，即最终时间 $t_f$

P.set_integral([1]) 
# 参数是关注的积分的表达式列表

然后，我们定义阶段约束（也称路径约束），约束在阶段时间间隔内的任何时间都满足。

P.set_phase_constraint([u], [-1], [1], [True])
# 第一个参数是阶段约束函数的表达式列表
# 第二个参数是约束函数的下界列表
# 第三个参数是约束函数的上界列表
# 第四个参数是一个布尔值列表，表示约束函数是否是 bang-bang 控制
# 这四个参数的长度必须相等

接下来，我们定义阶段的边界条件。

P.set_boundary_condition([1, 0], [0, 0], 0, None)
# 第一和第二个参数是状态变量的初始和终端边界条件，这些列表的长度必须等于状态变量的数量
# 第三和第四个参数是阶段的开始和结束时间
# 这些参数中的自由量设为 None

最后，我们指定阶段的离散参数

P.set_discretization(3, 5)
# 第一个参数是子区间的数量，
# 第二个参数是每个子区间中的插值节点数量。

当前定义完成了阶段级别的所有关系。可以将阶段设置到系统对象。

S.set_phase([P])

然后，我们定义系统的目标函数

S.set_objective(P.I[0])
# 在阶段上定义的积分自动生成了 Sympy 符号，通过 P.I 访问该列表

至此，我们已经完成了问题的设置，且 pockit 已自动生成离散问题。总结一下，我们通过以下步骤定义了问题

创建系统和阶段对象
设置每个阶段的动力学、积分、阶段约束、边界条件和离散参数。（其中动力学、边界条件和离散参数是必需的）
为系统对象设置阶段、目标函数和系统约束。（其中阶段和目标函数是必需的）

给出初始猜测

为了将离散化问题传递给 NLP 求解器，我们还需要提供一个解的初始猜测。初始猜测的质量对于优化是否成功至关重要。Pockit 使用 Variable 对象来处理每个阶段的离散优化变量。我们首先让 pockit 基于边界条件的线性插值生成一个猜测，作为起点。

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.t_x 中，控制变量的时间插值节点存储在属性 v.t_u 中。注意：

变量的插值节点不会均匀分布在区间内
状态变量和控制变量可能具有不同的插值节点

求解离散问题

最后，我们通过以下方式调用 NLP 求解器求解离散问题：

v, info = ipopt.solve(S, v)
print(info['status_msg'].decode())
# 或者使用 Scipy 后端:
# v, info = scipy.solve(S, v)
# print(info['success'])

绘图数值解

现在，可以通过以下方式访问解：

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

但是，这些数字很难阅读。正如上面提到的，插值节点存储在属性 v.t_x 和 v.t_u 中。绘图的正确方法是将变量与插值节点进行绘图，如下所示：

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()

生成的图如下所示：

非精确解，控制不是 Bang-Bang 控制

调整网格

不幸的是，从图中可以看出，解并不准确。如上所述，最优控制是 bang-bang 控制。Pockit 采用分段多项式插值的方法近似函数，因此只能在网格处允许间断。在本例中，理论分析表明控制的切换时间为 $\hat t_s = 1$，其恰好位于子区间的中间。对于一般问题，我们无法通过解析方法找到切换时间。Pockit 提供了一个利用多次迭代数值逼近切换时间的方法：

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"])
    # 或者使用 Scipy 后端:
    # v, info = scipy.solve(S, v)
    # assert info["success"]
    # print(info["fun"])

在代码片段中，我们首先将最大迭代次数设置为 10。然后，我们使用 System.check() 方法来检查解是否足够准确。容许误差范围是可调的，但为了简单起见，我们在这里使用默认值。如果解通过测试，则中断循环。否则，我们使用 System.refine()方法来调整离散参数。该方法根据当前优化变量调整每个子区间中的网格点和插值节点的数量。网格调整后，该方法返回一个新的Variable 对象，其结构与新的离散参数兼容，并且值从旧解中插值得到。然后，我们再次使用这个新变量作为初始猜测来求解问题。重复该过程，直到解通过测试或达到最大迭代次数。最终解为：

精确解，控制为 bang-bang 控制

解非常准确，目标函数取值为 t_f = 2.000000044988834，与理论解非常接近，此外，控制具有 bang-bang 控制特征。（图中 + 表示插值节点）

总结

这个简短的教程展示了如何使用 pockit 来解决一个简单的最优控制问题。我们通过将问题定义输入 pockit 来离散问题。然后，我们生成了解的初始猜测并求解了问题。最后，我们通过迭代调整网格来提高解的精度。同时，我们还展示了如何访问解并绘制解的图形。

进一步了解请参阅用户指南和 API 文档。

完整代码

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()