A rapid trajectory optimization method based on parallel computing

Tianyi WANG; Jisong ZHAO

doi:10.1051/jnwpu/20254340659

All issues

Volume 43 / No 4 (August 2025)

JNWPU, 43 4 (2025) 659-667

Full HTML

Open Access

Issue		JNWPU Volume 43, Number 4, August 2025


Page(s)		659 - 667
DOI		https://doi.org/10.1051/jnwpu/20254340659
Published online		07 October 2025

JNWPU 2025, 43(4): 659–667

A rapid trajectory optimization method based on parallel computing

基于并行计算的轨迹快速优化方法研究

Tianyi WANG (王天一) and Jisong ZHAO (赵吉松)

School of Astronautics, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China

Received: 24 August 2024

Abstract

The direct collocation method transforms a trajectory optimization problem into a nonlinear programming (NLP) problem by discretizing both control and state variables. During the NLP solution process, repeated calculations of the first and second derivatives of the NLP and the values of the dynamic system at each discrete point are required, leading to great computational complexities. Therefore, this paper proposes the following method: First, the hyper-dual number method is introduced to accurately identify the sparsity of the second-derivative matrix of the NLP and to determine the locations of the non-zero elements. Then, a multi-core parallel approach is used to rapidly compute the non-zero elements of the first and second derivatives of the NLP as well as the values of the dynamic system at each discrete point. Finally, OpenMP is employed for programming calculation in the C++ environment to further enhance computational efficiency from the perspective of programming language. Simulation results demonstrate that the proposed method effectively improves the efficiency of trajectory optimization and its computational efficiency without compromising accuracy.

摘要

直接配点法通过对控制变量和状态变量都进行离散将轨迹优化问题转化为非线性规划(nonlinear programming, NLP)进行求解。在求解NLP时, 需要反复计算NLP的一阶/二阶偏导数和动力学系统在各离散点处的值, 计算量比较大。针对该问题, 提出如下解决策略: 引入超对偶数方法准确识别NLP二阶偏导数矩阵的稀疏型, 确定其中非零元素的位置; 采用多核并行方式快速计算NLP的一阶/二阶偏导数的非零元素以及动力学系统在各离散点处的值; 在C++环境下采用OpenMP方式进行编程计算, 从编程语言角度进一步提高计算效率。仿真结果表明, 文中方法给出的策略在不影响精度的情况下, 均能显著提高轨迹优化效率。

Key words: trajectory optimization / nonlinear programming / sparsity / hyper-dual numbers / parallel computing

关键字 : 轨迹优化 / 非线性规划 / 稀疏 / 超对偶数 / 并行计算

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

轨迹优化在航空航天领域应用广泛, 具有非常重要的意义和实际应用价值[1]。轨迹优化本质上是一个最优控制问题, 目前广泛使用的求解方法有2种, 即直接法和间接法。间接法基于变分法或极大值原理, 而直接法是通过对状态变量和控制变量的离散化构造最优控制问题的近似, 从而将最优控制问题转化为非线性规划问题。直接法中的配点法[2]由于不需要推导最优性必要条件, 对初值的敏感性较低且易收敛, 近年来得到了广泛研究[3–4]。

目前广泛使用的NLP求解器SNOPT[5]和IPOPT[6]都需要反复计算NLP的一阶/二阶偏导数, 因此, 提高NLP一阶/二阶偏导数的计算效率对于提升轨迹优化的效率至关重要。Agamawi和Patterson等[7–8]进行了深入研究, 发现轨迹优化问题离散得到的NLP具有高度的稀疏性, 即NLP的一阶/二阶偏导数矩阵中包含大量零元素。赵吉松[9]研究了局部配点法离散得到的NLP一阶偏导数的稀疏特性, 建立了一阶偏导数的高效计算方法。此外, 为提高偏导数的计算效率, Fike等[10]提出了基于超对偶数的偏导数计算方法(简称超对偶数方法), 相比于目前广泛使用的有限差分法、自动微分法[11]等偏导数计算方法, 该方法可以迅速准确地计算一阶和二阶偏导数。Agamawi等[12]使用该方法为NLP求解器提供准确的偏导数, 但他们在对偏导数矩阵进行稀疏性分析时是通过一阶偏导数矩阵的稀疏型近似得到二阶偏导数矩阵的稀疏型[13], 该方法准确度较低, 会错误识别出额外的二阶偏导数, 增加偏导数矩阵非零元素的个数, 从而增大计算量, 最终降低轨迹优化问题求解的速度和精度。因此, 本文改进了NLP二阶偏导数矩阵的稀疏性分析方法, 与Rao等人使用超对偶数方法进行实际偏导数计算不同, 本文侧重于使用该方法去准确识别NLP二阶偏导数矩阵的稀疏型。

此外, 通过研究发现, NLP的一阶/二阶偏导数矩阵中偏导数的计算以及动力学系统在各离散点处的计算都是相互独立的, 非常适合进行并行计算。并行计算是相对于串行计算而言, 通过同时使用多个计算资源来解决计算问题的过程, 是提高计算速度和处理能力的一种有效手段[14]。目前广泛使用的并行计算方法有多种, 例如基于GPU的CUDA、基于CPU的MPI、OpenMP, 这些方法在轨迹优化领域也有较为广泛的应用。Betts等[15]提出了一种使用并行处理器来减少这些计算成本的技术, 将轨迹分解成多个阶段, 可以并行模拟, 从而降低单个轨迹的成本。Antony等[16]基于CUDA平台, 提出了一种基于GPU加速的间接弹道优化方法, 由此产生了快速最优控制框架, 能够构建高质量的最优轨迹。钟新玉等[17]基于MPI消息传递机制实现了弹道仿真计算模型的并行编程, 获得了较高的加速比和计算效率, 与单机运行相比, 计算机集群可以大大缩短计算时间, 提高计算效率。孟振华[18]以组合导航、轨迹规划、图像制导等为背景进行基于共享存储的OpenMP并行编程模型的并行设计和优化设计, 显著缩短了各类算法的计算时间。与MPI和CUDA相比, OpenMP在便于实现的同时兼具较高的计算效率, 在处理规模不是特别大的计算任务时比前面2种方法更有优势。Agamawi和Rao在文献[12]中使用串行计算偏导数, 而并行计算相比于串行计算可以进一步提升计算效率。因此本文选择使用OpenMP并行计算NLP一阶/二阶偏导数矩阵中的非零元素。具体而言, 在得到NLP一阶/二阶偏导数矩阵的稀疏型后, 通过OpenMP将稀疏型中的非零元素分配到不同的线程中进行并行计算, 最后再将计算完成的各非零元素放回其原来的位置, 组装得到NLP的一阶/二阶偏导数矩阵。此外, 通过使用OpenMP并行计算动力学系统在各离散点处的值, 可进一步提升计算效率。

本文结合超对偶数方法和OpenMP并行计算的综合策略不仅提高了计算效率, 还保留了计算的精确性, 为解决较为复杂的轨迹优化问题提供了一种高效且可靠的途径。

1 轨迹优化问题的数学描述

轨迹优化问题本质上属于最优控制问题, 以Bolza型最优控制问题为例, 可描述为: 求解控制变量u(t)∈R^m使得(1)式所示目标函数最小化。

$J=M\left(\boldsymbol{x}\left(t_0\right), t_0, \boldsymbol{x}\left(t_{\mathrm{f}}\right), t_{\mathrm{f}}\right)+\int_{t_0}^{t_{\mathrm{f}}} L(\boldsymbol{x}(t), \boldsymbol{u}(t), t) \mathrm{d} t$ (1)

式中，M: Rⁿ×R×Rⁿ×R→R, L: Rⁿ×R^m×R→R, x∈Rⁿ, u∈R^m, t∈[t₀, t_f]。

状态方程为

$\dot{\boldsymbol{x}}=f(\boldsymbol{x}(t), \boldsymbol{u}(t), t), t \in\left[t_0, t_{\mathrm{f}}\right]$ (2)

端点条件为

$\boldsymbol{E}\left(\boldsymbol{x}\left(t_0\right), t_0, \boldsymbol{x}\left(t_{\mathrm{f}}\right), t_{\mathrm{f}}\right)=0$ (3)

路径约束为

$\boldsymbol{C}(\boldsymbol{x}(t), \boldsymbol{u}(t), t) \leqslant 0, t \in\left[t_0, t_{\mathrm{f}}\right]$ (4)

式中f: Rⁿ×R^m×R→Rⁿ, E: Rⁿ×R×R^m×R→R^e, C: Rⁿ×R^m×R→R^c。(1)~(4)式所描述的问题称为连续Bolza型最优控制问题。

2 基于超对偶数的稀疏性分析方法

2.1 NLP二阶偏导数矩阵

为构造NLP二阶偏导数矩阵, 定义函数F, F包括状态方程f、端点条件E和路径约束C三部分。定义总二阶偏导数矩阵如(5)式所示。

$\begin{aligned} & \boldsymbol{H}=\left[\begin{array}{lll}\boldsymbol{H}_1 & & \\ \boldsymbol{H}_2 & \boldsymbol{H}_3 & \\ \boldsymbol{H}_4 & \boldsymbol{H}_5 & \boldsymbol{H}_6\end{array}\right] \\ & \boldsymbol{H}_1=\left[\begin{array}{cccc}\frac{\partial^2 F}{\partial x_{1 i}^2} & & & \\ \frac{\partial^2 F}{\partial x_{1 i} \partial x_{2 i}} & \frac{\partial^2 F}{\partial x_{2 i}^2} & & \\ \vdots & \vdots & & \\ \frac{\partial^2 F}{\partial x_{1 i} \partial x_{n i}} & \frac{\partial^2 F}{\partial x_{2 i} \partial x_{n i}} & \cdots & \frac{\partial^2 F}{\partial x_{n i}^2}\end{array}\right] \\ & \boldsymbol{H}_2=\left[\begin{array}{cccc}\frac{\partial^2 F}{\partial x_{1 i} \partial u_{1 i}} & \frac{\partial^2 F}{\partial x_{2 i} \partial u_{1 i}} & \cdots & \frac{\partial^2 F}{\partial x_{n i} \partial u_{1 i}} \\ \frac{\partial^2 F}{\partial x_{1 i} \partial u_{2 i}} & \frac{\partial^2 F}{\partial x_{2 i} \partial u_{2 i}} & \cdots & \frac{\partial^2 F}{\partial x_{n i} \partial u_{2 i}} \\ \vdots & \vdots & & \vdots \\ \frac{\partial^2 F}{\partial x_{1 i} \partial u_{m i}} & \frac{\partial^2 F}{\partial x_{2 i} \partial u_{m i}} & \cdots & \frac{\partial^2 F}{\partial x_{n i} \partial u_{m i}}\end{array}\right] \\& \boldsymbol{H}_3 =\left[\begin{array}{cccc} \frac{\partial^2 F}{\partial u_{1 i}^2} & & \\ \frac{\partial^2 F}{\partial u_{1 i} \partial u_{2 i}} & \frac{\partial^2 F}{\partial u_{2 i}^2} & \\ \vdots & \vdots & \\ \frac{\partial^2 F}{\partial u_{1 i} \partial u_{m i}} & \frac{\partial^2 F}{\partial u_{2 i} \partial u_{m i}} & \cdots & \frac{\partial^2 F}{\partial u_{m i}^2} \end{array}\right] \\ &\boldsymbol{H}_4 =\left[\begin{array}{llll} \frac{\partial^2 F}{\partial x_{1 i} \partial t_{0 i}} & \frac{\partial^2 F}{\partial x_{2 i} \partial t_{0 i}} & \cdots & \frac{\partial^2 F}{\partial x_{n i} \partial t_{0 i}} \\ \frac{\partial^2 F}{\partial x_{1 i} \partial t_{{\rm f} i}} & \frac{\partial^2 F}{\partial x_{2 i} \partial t_{{\rm f} i}} & \cdots & \frac{\partial^2 F}{\partial x_{n i} \partial t_{{\rm f} i}} \end{array}\right] \\ &\boldsymbol{H}_5 =\left[\begin{array}{llll} \frac{\partial^2 F}{\partial u_{1 i} \partial t_{0 i}} & \frac{\partial^2 F}{\partial u_{2 i} \partial t_{0 i}} & \cdots & \frac{\partial^2 F}{\partial u_{m i} \partial t_{0 i}} \\ \frac{\partial^2 F}{\partial u_{1 i} \partial t_{{\rm f} i}} & \frac{\partial^2 F}{\partial u_{2 i} \partial t_{{\rm f} i}} & \cdots & \frac{\partial^2 F}{\partial u_{m i} \partial t_{{\rm f} i}} \end{array}\right] \\& \boldsymbol{H}_6 =\left[\begin{array}{cc} \frac{\partial^2 F}{\partial t_{0 i}^2} & \\ \frac{\partial^2 F}{\partial t_{0 i} \partial t_{{\rm f} i}} & \frac{\partial^2 F}{\partial t_{{\rm f} i}^2} \end{array}\right] \end{aligned}$ (5)

式中：x_1i, x_2i, …, x_ni为第i个节点处的n个状态变量; u_1i, u_2i, …, u_mi为第i个节点处的m个控制变量; t_0i为第i个节点处的初始时刻; t_fi为第i个节点处的终端时刻。以x_1i为例, 1表示第1个状态变量x₁, i表示第i个节点, 其余变量同理。上述二阶偏导数矩阵在各个节点处的形式都是相同的。

提供准确的NLP二阶偏导数矩阵可以显著提高NLP求解器的计算效率, 然而得到精确的二阶偏导数矩阵却非常困难。考虑到NLP二阶偏导数矩阵具有高度的稀疏性, 如果准确识别其稀疏型, 可以有效提升NLP问题求解的效率和准确性。目前流行的通用轨迹优化软件GPOPS-Ⅱ已经可以准确地识别NLP一阶偏导数矩阵稀疏型, 但在识别NLP二阶偏导数矩阵的稀疏型时还不够准确。

以函数z=x+y为例, GPOPS-Ⅱ计算一阶导数

$\frac{\partial z}{\partial x}=1, \quad \frac{\partial z}{\partial y}=1$ (6)

由于GPOPS-Ⅱ对各变量存在一阶导数, 即使z对变量x和y的二阶偏导数都为0, GPOPS-Ⅱ也会推导二阶偏导数存在且不为0, 这就导致GPOPS-Ⅱ在识别二阶偏导数矩阵的稀疏型时过于保守, 将零元素识别成非零元素。虽然GPOPS-Ⅱ的稀疏型识别方法会识别出额外的非零元素, 但在进行非零元素的计算时, 该元素的计算结果还是为0, 所以并不会影响最终的寻优结果, 只会影响寻优速度。

GPOPS-Ⅱ通过一阶偏导数矩阵的稀疏型近似得到二阶偏导数矩阵的稀疏型, 该方法准确度较低, 会识别出额外的二阶偏导数, 增加偏导数矩阵非零元素的个数, 从而增大计算量, 最终降低轨迹优化问题求解的速度和精度。

2.2 基于超对偶数的偏导数计算方法

对于多变量函数, 基于超对偶数的二阶偏导数计算如(7)式所示。

$\frac{\partial^2 \boldsymbol{f}(\boldsymbol{x})}{\partial x_i \partial x_j}=\frac{\varepsilon_1 \varepsilon_2 \operatorname{part}\left[\boldsymbol{f}\left(\boldsymbol{x}+h_1 \varepsilon_1 \boldsymbol{e}_i+h_2 \varepsilon_2 \boldsymbol{e}_j+0 \varepsilon_1 \varepsilon_2\right)\right]}{h_1 h_2}$ (7)

式中: f(x)为问题的函数，f: Rⁿ→R^m, x为n维向量; h₁, h₂为步长, 本文设置为1;e_i, e_j为2个单位向量。ε₁ε₂part[·]表示提取函数中ε₁ε₂的系数, 即f(x)对x_ix_j的m维偏导数向量。

下面演示偏导数的计算过程。

给定一个n维自变量x和2个n维单位向量e_i, e_j，如(8)式所示。

$\boldsymbol{x}=\left[\begin{array}{c}x_1 \\ x_2 \\ \vdots \\ x_n\end{array}\right], \boldsymbol{e}_i=\left[\begin{array}{c}\mathbf{0}_{(i-1) \times 1} \\ 1 \\ \mathbf{0}_{(n-i) \times 1}\end{array}\right], \boldsymbol{e}_j=\left[\begin{array}{c}\mathbf{0}_{(j-1) \times 1} \\ 1 \\ \mathbf{0}_{(n-j) \times 1}\end{array}\right]$ (8)

式中，单位向量e_i的第i个分量为1，其余分量均为0；单位向量e_j的第j个分量为1，其余分量均为0。

利用上述向量与零向量构造的超对偶数向量为新的函数自变量向量，包括4个部分：①n维自变量x；②单位向量e_i乘以复数单位ε₁；③单位向量e_j乘以复数单位ε₂；④零向量乘以复数单位ε₁ε₂，如(9)式所示。

$\boldsymbol{x}=\left[\begin{array}{c}x_1 \\ x_2 \\ \vdots \\ x_i \\ \vdots \\ x_j \\ \vdots \\ x_n\end{array}\right]+\left[\begin{array}{c}\mathbf{0}_{(i-1) \times 1} \\ 1 \\ \mathbf{0}_{(n-i) \times 1}\end{array}\right] \boldsymbol{\varepsilon}_1+\left[\begin{array}{c}\mathbf{0}_{(j-1) \times 1} \\ 1 \\ \mathbf{0}_{(n-j) \times 1}\end{array}\right] \boldsymbol{\varepsilon}_2+\left[\begin{array}{c}0 \\ 0 \\ \vdots \\ 0 \\ \vdots \\ 0 \\ \vdots \\ 0\end{array}\right] \boldsymbol{\varepsilon}_1 \boldsymbol{\varepsilon}_2$ (9)

将新的函数自变量向量代入函数f中得到

$\boldsymbol{f}(\boldsymbol{x})=\left[\begin{array}{c}f_1(\boldsymbol{x}) \\ f_2(\boldsymbol{x}) \\ \vdots \\ f_m(\boldsymbol{x})\end{array}\right]$ (10)

式中, f(x)为得到的m维函数向量, 该向量由m个超对偶数向量分量组成。

提取函数输出中ε₁ε₂的系数, 便可得到函数f对x_ix_j的二阶偏导数，如(11)式所示。

$\frac{\partial^2 \boldsymbol{f}(\boldsymbol{x})}{\partial x_i \partial x_j}=\left[\begin{array}{c}\frac{\partial^2 f_1(\boldsymbol{x})}{\partial x_i \partial x_j} \\ \frac{\partial^2 f_2(\boldsymbol{x})}{\partial x_i \partial x_j} \\ \vdots \\ \frac{\partial^2 f_m(\boldsymbol{x})}{\partial x_i \partial x_j}\end{array}\right]=\varepsilon_1 \varepsilon_2 \mathrm{part}\left[\begin{array}{c}f_1(\boldsymbol{x}) \\ f_2(\boldsymbol{x}) \\ \vdots \\ f_m(\boldsymbol{x})\end{array}\right]$ (11)

2.3 稀疏型对比

为了描述H的稀疏型, 定义struct函数[9]

$\operatorname{struct}(x)= \begin{cases}1 & x \neq 0 \\ 0 & x=0\end{cases}$ (12)

记作

$\boldsymbol{S}=\operatorname{struct}(\boldsymbol{H})$ (13)

式中, struct(H)表示对H的每个元素进行struct运算。S表示H的稀疏型。为了得到S, 不需要计算H的每个元素具体值, 只需要判断是否为0。本文使用的超对偶数方法是通过随机选取多个实部值进行偏导数计算来判断某个元素是否为0, 进而识别出偏导数矩阵的稀疏型, 并非计算元素的准确值, 在进行非零元素计算时依然基于中心差分法。

下面以轨道转移问题为例, 将GPOPS-Ⅱ方法与超对偶数方法识别出的稀疏型进行对比, 轨道转移问题的动力学方程在4.1节中给出。表 1为使用不同方法识别出的二阶偏导数矩阵稀疏型的对比。

表1

轨道转移问题NLP二阶偏导数矩阵稀疏型比较

显然, 使用超对偶数方法, 能够更为精确地识别二阶偏导数矩阵的稀疏型, 与手工解析方法得到的稀疏型完全相同。从表 1中可以看出, 使用GPOPS-Ⅱ方法得到的稀疏型中非零元素个数为30, 而使用超对偶数方法得到的稀疏型中非零元素个数仅为14, 与GPOPS-Ⅱ方法相比减少了16个非零元素。综上所述, 使用超对偶数方法减少了冗余非零元素的计算, 提高了轨迹优化的效率。

3 基于OpenMP的并行计算

如图 1所示, NLP偏导数矩阵中包含大量元素。

图1

NLP一阶偏导数矩阵示意图

其中, 空白元素表示恒为零, “●”表示非零元素, “×”表示可能不为零的元素。

本文使用OpenMP作为并行计算方法。OpenMP是一种基于共享内存的多线程并行编程模型, 专为多核处理器设计, 通过使用OpenMP, 程序能够同时对独立的计算子块进行并行处理, 从而显著提高计算速度和程序的整体运行效率。

具体来说, 本文首先准确识别NLP一阶/二阶偏导数矩阵的稀疏型, 然后使用OpenMP将其中的非零元素分配到多个线程中进行计算。如图 2所示, 假设需要计算的非零元素的个数为m, 通过OpenMP将这m个非零元素分配到k个线程中并行计算, 在每个线程计算完成后收集结果, 组装得到NLP一阶/二阶偏导数矩阵。

图2

非零元素并行计算流程示意图

此外, 除了NLP的一阶/二阶偏导数, 动力学系统在各离散点处的值也可以进行并行计算。同理, 将各离散点处的值也分配到这k个线程中进行计算, 可以进一步提升计算效率。

4 仿真与分析

4.1 轨道转移问题

轨道转移问题是航空航天领域的经典问题, 需要在满足物理约束的前提下找到最优的转移策略。

该轨迹优化问题的状态变量分别为矢径r、轨道角度θ、径向速度v_r、切向速度v_θ, 控制变量为u_r和u_θ。该轨迹优化问题的描述如下: 求解1个最优控制u^*(t)=(u_r(t), u_θ(t))^T, 使系统在末端时刻矢径最大化。因此, 目标函数为

$J=\min \left(-r_{t_{\mathrm{f}}}\right)$ (14)

满足动力学方程

$\left\{\begin{array}{l}\dot{r}=v_r \\ \dot{\theta}=\frac{v_\theta}{r} \\ \dot{v}_r=\frac{v_\theta^2}{r}-\frac{\mu}{r^2}+\frac{T}{m} \mu_r \\ \dot{v}_\theta=\frac{v_r v_\theta}{r}+\frac{T}{m} \mu_\theta\end{array}\right.$ (15)

式中：μ为太阳引力常数；T为发动机推力大小；m为飞行器质量，它会因为燃料和其他部件的消耗而随时间逐渐变化。m可以表示为

$m=m\left(t_0\right)-m_{\mathrm{s}} \cdot t$ (16)

式中：t₀为初始时间；m_s为飞行器质量变化速度, 为简化研究, 设其为常值。

边界条件为

$\boldsymbol{x}\left(t_0\right)=x_0, \boldsymbol{x}\left(t_{\mathrm{f}}\right)=x_{\mathrm{f}}$ (17)

在轨道转移的过程中, 状态变量和控制变量需要满足约束

$\boldsymbol{x}_{\min } \leqslant \boldsymbol{x}(t) \leqslant \boldsymbol{x}_{\max }$ (18)

$\boldsymbol{u}_{\min } \leqslant \boldsymbol{u}(t) \leqslant \boldsymbol{u}_{\max }$ (19)

同时, 由于推力受限, 在轨道转移中需要满足路径约束

$\|\boldsymbol{u}\|_2=\frac{T}{m\left(t_0\right)}$ (20)

初始条件是t₀=0, m(t₀)=1 000 kg, r(t₀)=1 AU, θ(t₀)=0, v_r(t₀)=0, v_θ(t₀)=1 AU/TU, T=0.140 5 m₀V₀²/r₀, μ=1 AU³/TU³。终端约束是t_f=3.32 TU, v_r(t_f)=0。其中AU为天文单位, 指地球到太阳的平均距离, 在这里取1.495 978 7×10¹¹ m。在本算例中, TU(time unit)为时间单位，其取值为一年，即365.25 d。

项目组前期开发的轨迹优化软件[19]使用的优化算法为梯形格式的配点法, 使用GPOPS-Ⅱ中的偏导数矩阵稀疏型识别方法, 并且不使用OpenMP并行计算(简称对比方法)。本文方法使用的优化算法为梯形格式的配点法, 使用超对偶数方法识别偏导数矩阵的稀疏型, 并使用OpenMP并行计算一阶/二阶偏导数矩阵中的非零元素以及动力学系统在各离散点处的值。

本文所使用的计算平台的配置为Intel Core i3-10100 3.60 GHz处理器, 8 GB RAM内存, Windows 10家庭版64位操作系统。本算例使用的优化方法均取101个节点进行求解, OpenMP采用4线程并行。为确保结果的可靠性和准确性, 最终的优化耗时结果是基于10次优化的平均耗时。仿真结果如表 2所示, 算例的状态变量随时间变化曲线如图 3所示, 最优控制方案输出结果如图 4所示。

表2

不同优化方法的优化结果

图3

状态变量随时间变化曲线

图4

控制量随时间变化曲线

由优化结果可以看出, 每种优化方法的目标函数都是相同的, 均为1.525 22 AU, 这表明3种方法在优化结果的精度上是等效的。但值得注意的是, 本文方法的优化耗时明显短于GPOPS-Ⅱ和对比方法, 求解时间仅为715 ms, 相比于GPOPS-Ⅱ大幅减少了约85.8%, 相比于对比方法减少了35.7%。

4.2 超声速飞行器再入问题

超声速飞行器再入问题的状态变量包括飞行器到地心的距离r, 速度v, 航迹倾角γ, 航迹偏角ψ, 经度θ, 纬度ϕ; 控制变量包括攻角α和侧倾角σ。

超声速飞行器再入问题的目标函数为

$J=\min \left(\phi_{\mathrm{f}}\right)$ (21)

满足动力学方程

$\left\{\begin{array}{l}\frac{\mathrm{d} r}{\mathrm{~d} t}=v \sin \gamma \\ \frac{\mathrm{~d} v}{\mathrm{~d} t}=-\frac{D}{m}-g \sin \gamma \\ \frac{\mathrm{~d} \gamma}{\mathrm{~d} t}=\frac{L \cos \sigma}{m v}-\frac{g \cos \gamma}{v}+\frac{v \cos \gamma}{r} \\ \frac{\mathrm{~d} \psi}{\mathrm{~d} t}=\frac{L \sin \sigma}{m v \cos \gamma}-\frac{v \cos \gamma \cos \psi \tan \phi}{r} \\ \frac{\mathrm{~d} \theta}{\mathrm{~d} t}=\frac{v \cos \gamma \cos \psi}{r \cos \phi} \\ \frac{\mathrm{~d} \phi}{\mathrm{~d} t}=\frac{v \cos \gamma \sin \psi}{r}\end{array}\right.$ (22)

式中：m为质量；L和D分别为升力和阻力；g为重力加速度；r=h+R_e, h为高度；R_e为地球半径。

升力和阻力的计算公式为

$\left\{\begin{array}{l}L=\frac{1}{2} \rho v^2 S C_L \\ D=\frac{1}{2} \rho v^2 S C_D\end{array}\right.$ (23)

式中：C_L为升力系数；C_D为阻力系数；ρ为大气密度；S为参考面积。

升力系数和阻力系数计算公式如(24)式所示。

$\left\{\begin{array}{l}C_L=C_{L 0}+C_{L 1} \cdot \alpha \\ C_D=C_{D 0}+C_{D 1} \cdot \alpha+C_{D 2} \cdot \alpha^2\end{array}\right.$ (24)

重力加速度满足

$g=\frac{\mu}{r^2}$ (25)

式中, μ为地球引力常数。

初始条件h₀=79 248 m, v₀=7 333.792 8 m/s, γ₀=-0.018 57 rad, ψ₀=0 rad, θ₀=0 rad, ϕ₀=0 rad。终端条件为h_f=24 384 m, v_f=762 m/s, γ_f=-0.087 27 rad。优化方法均取101个节点对本算例进行求解, OpenMP采用4线程并行, 为确保结果的可靠性和准确性, 最终的优化耗时结果是基于10次优化的平均耗时。仿真结果如表 3所示, 状态变量随时间变化曲线如图 5所示, 控制变量随时间变化曲线如图 6所示。

表3

不同优化方法得到的优化结果

图5

状态变量随时间变化曲线

图6

控制变量随时间变化曲线

从优化结果可以看出, 2种方法的目标函数都是相同的, 均为-0.501 63 rad, 但本文方法的优化耗时明显少于对比方法, 求解时间为1 562 ms, 相比于对比方法减少了30.8%。

通过以上2个算例可以看出, 结合超对偶数方法和OpenMP的综合计算方法, 可以在保证计算精度的同时, 较大程度地提升轨迹优化问题的求解效率。超对偶数方法和OpenMP并行计算方法的引入对轨迹优化问题的求解效率有着显著影响, 结合2种方法能够较大限度地发挥计算资源的优势, 缩短计算时间, 提高计算效率。

5 结论

本文深入分析了飞行器轨迹优化问题, 结合超对偶数方法与OpenMP并行计算方法, 显著提高了轨迹优化效率。具体而言, 本文通过引入超对偶数方法准确识别了NLP二阶偏导数矩阵的稀疏型, 确定其中非零元素的位置, 有效地解决了二阶偏导数的冗余计算问题, 而引入OpenMP则充分利用了其多核处理器的优势, 通过并行计算NLP的一阶/二阶偏导数的非零元素以及动力学系统在各离散点处的值, 大幅缩短了计算时间。仿真结果表明, 采用这种综合方法不仅保证了计算精度, 而且提升了计算效率。对于轨道转移问题, 优化时间相比于GPOPS-Ⅱ大幅减少了约85.8%, 相比于对比方法减少了35.7%。对于超声速飞行器再入问题, 优化时间与对比方法相比减少了30.8%。相比于对比方法, 本文方法在处理不同轨迹优化问题时均有显著的效率提升。

未来的研究可以探索更高效的并行计算策略, 进一步提升计算速度。同时, 还可以将这一方法扩展到其他类型的优化问题中, 以验证其通用性和适应性。

References

ZHANG Can, HU Dongdong, YE Lei, et al. Overview of foreign hypersonic vehicle technology development in 2017[J]. Tactical Missile Technology, 2018(1): 47–50 (in Chinese) [Google Scholar]
BETTS J T. Practical methods for optimal control and estimation using nonlinear programming[M]. Philadelphia: Society for Industrial and Applied Mathematics, 2010. [Google Scholar]
ZHAO Jisong, ZHANG Jianhong, LI Shuang. Rapid and high accuracy approach for hypersonic glide vehicle reentry trajectory optimization[J]. Journal of Astronautics, 2019, 40(9): 1034–1043 (in Chinese) [Google Scholar]
REN Pengfei, WANG Hongbo, ZHOU Guofeng. Reentry trajectory optimization of hypersonic vehicle based on adaptive pseudo-spectrum method[J]. Journal of Beijing University of Aeronautics and Astronautics, 2019, 45(11): 2257–2265 (in Chinese) [Google Scholar]
GILL P E, MURRAY W, SAUNDERS M A. SNOPT: an SQP algorithm for large-scale constrained optimization[J]. SIAM Journal on Optimization, 2002, 12(4): 979–1006 [Article] [Google Scholar]
BIEGLER L T, ZAVALA V M. Large-scale nonlinear programming using IPOPT: an integrating framework for enterprise-wide dynamic optimization[J]. Computers & Chemical Engineering, 2009, 33(3): 575–582 [Google Scholar]
AGAMAWI Y M, RAO A V. Comparison of derivative estimation methods in optimal control using direct collocation[J]. AIAA Journal, 2020, 58(1): 341–354 [Article] [Google Scholar]
PATTERSON M A, RAO A V. Exploiting sparsity in direct collocation pseudospectral methods for solving optimal control problems[J]. Journal of Spacecraft and Rockets, 2012, 49(2): 364–377 [Google Scholar]
ZHAO Jisong. Exploiting sparsity in local collocation methods for solving trajectory optimization problems[J]. Journal of Astronautics, 2017, 38(12): 1263–1272 (in Chinese) [Google Scholar]
FIKE J, ALONSO J. The development of hyper-dual numbers for exact second-derivative calculations[C]//49th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 2011: 886 [Google Scholar]
CSENDES T. Developments in Reliable Computing[M]. Dordrecht: Kluwer Academic Publishers, 1999, 77–104 [Google Scholar]
AGAMAWI Y M, RAO A V. CGPOPS: a C++software for solving multiple-phase optimal control problems using adaptive Gaussian quadrature collocation and sparse nonlinear programming[J]. ACM Transactions on Mathematical Software, 2020, 46(3): 1–38 [Google Scholar]
PATTERSON M A, RAO A V. GPOPS-Ⅱ: a Matlab software for solving multiple-phase optimal control problems using hp-adaptive Gaussian quadrature collocation methods and sparse nonlinear programming[J]. ACM Trans on Mathematical Software, 2014, 41(1): 1–37 [Google Scholar]
DU Zhihui. Parallel programming techniques in high-performance computing: MPI parallel programming[M]. Beijing: Tsinghua University Press, 2001 (in Chinese) [Google Scholar]
BETTS J T, HUFFMAN W P. Trajectory optimization on a parallel processor[J]. Journal of Guidance, Control, and Dynamics, 1991, 14(2): 431–439 [Article] [Google Scholar]
ANTONY T, GRANT M J. Rapid indirect trajectory optimization on highly parallel computing architectures[J]. Journal of Spacecraft and Rockets, 2017, 54(5): 1081–1091 [Article] [Google Scholar]
ZHONG Xinyu, WANG Xun, ZHANG Xuejun. The parallel computing research of trajectory simulation based on MPI[J]. Aerospace Control, 2015, 33(3): 63–67 (in Chinese) [Google Scholar]
MENG Zhenhua. Research on high performance parallel computing architecture based on FPGA+DSP[D]. Beijing: The First Academy of China Aerospace Science and Technology Corporation, 2018 (in Chinese) [Google Scholar]
ZHU Boling. Research and program development of general trajectory optimization method based on collocation method[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2023 (in Chinese) [Google Scholar]

All Tables

表1

轨道转移问题NLP二阶偏导数矩阵稀疏型比较

不同优化方法的优化结果

不同优化方法得到的优化结果

All Figures

	图1 NLP一阶偏导数矩阵示意图
In the text

	图2 非零元素并行计算流程示意图
In the text

	图3 状态变量随时间变化曲线
In the text

	图4 控制量随时间变化曲线
In the text

	图5 状态变量随时间变化曲线
In the text

	图6 控制变量随时间变化曲线
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.

[1] ZHANG Can, HU Dongdong, YE Lei, et al. Overview of foreign hypersonic vehicle technology development in 2017[J]. Tactical Missile Technology, 2018(1): 47–50 (in Chinese) [Google Scholar]

[2] BETTS J T. Practical methods for optimal control and estimation using nonlinear programming[M]. Philadelphia: Society for Industrial and Applied Mathematics, 2010. [Google Scholar]

[3] ZHAO Jisong, ZHANG Jianhong, LI Shuang. Rapid and high accuracy approach for hypersonic glide vehicle reentry trajectory optimization[J]. Journal of Astronautics, 2019, 40(9): 1034–1043 (in Chinese) [Google Scholar]

[4] REN Pengfei, WANG Hongbo, ZHOU Guofeng. Reentry trajectory optimization of hypersonic vehicle based on adaptive pseudo-spectrum method[J]. Journal of Beijing University of Aeronautics and Astronautics, 2019, 45(11): 2257–2265 (in Chinese) [Google Scholar]

[5] GILL P E, MURRAY W, SAUNDERS M A. SNOPT: an SQP algorithm for large-scale constrained optimization[J]. SIAM Journal on Optimization, 2002, 12(4): 979–1006 [Article] [Google Scholar]

[6] BIEGLER L T, ZAVALA V M. Large-scale nonlinear programming using IPOPT: an integrating framework for enterprise-wide dynamic optimization[J]. Computers & Chemical Engineering, 2009, 33(3): 575–582 [Google Scholar]

[7] AGAMAWI Y M, RAO A V. Comparison of derivative estimation methods in optimal control using direct collocation[J]. AIAA Journal, 2020, 58(1): 341–354 [Article] [Google Scholar]

[8] PATTERSON M A, RAO A V. Exploiting sparsity in direct collocation pseudospectral methods for solving optimal control problems[J]. Journal of Spacecraft and Rockets, 2012, 49(2): 364–377 [Google Scholar]

[9] ZHAO Jisong. Exploiting sparsity in local collocation methods for solving trajectory optimization problems[J]. Journal of Astronautics, 2017, 38(12): 1263–1272 (in Chinese) [Google Scholar]

[10] FIKE J, ALONSO J. The development of hyper-dual numbers for exact second-derivative calculations[C]//49th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 2011: 886 [Google Scholar]

[11] CSENDES T. Developments in Reliable Computing[M]. Dordrecht: Kluwer Academic Publishers, 1999, 77–104 [Google Scholar]

[12] AGAMAWI Y M, RAO A V. CGPOPS: a C++software for solving multiple-phase optimal control problems using adaptive Gaussian quadrature collocation and sparse nonlinear programming[J]. ACM Transactions on Mathematical Software, 2020, 46(3): 1–38 [Google Scholar]

[13] PATTERSON M A, RAO A V. GPOPS-Ⅱ: a Matlab software for solving multiple-phase optimal control problems using hp-adaptive Gaussian quadrature collocation methods and sparse nonlinear programming[J]. ACM Trans on Mathematical Software, 2014, 41(1): 1–37 [Google Scholar]

[14] DU Zhihui. Parallel programming techniques in high-performance computing: MPI parallel programming[M]. Beijing: Tsinghua University Press, 2001 (in Chinese) [Google Scholar]

[15] BETTS J T, HUFFMAN W P. Trajectory optimization on a parallel processor[J]. Journal of Guidance, Control, and Dynamics, 1991, 14(2): 431–439 [Article] [Google Scholar]

[16] ANTONY T, GRANT M J. Rapid indirect trajectory optimization on highly parallel computing architectures[J]. Journal of Spacecraft and Rockets, 2017, 54(5): 1081–1091 [Article] [Google Scholar]

[17] ZHONG Xinyu, WANG Xun, ZHANG Xuejun. The parallel computing research of trajectory simulation based on MPI[J]. Aerospace Control, 2015, 33(3): 63–67 (in Chinese) [Google Scholar]

[18] MENG Zhenhua. Research on high performance parallel computing architecture based on FPGA+DSP[D]. Beijing: The First Academy of China Aerospace Science and Technology Corporation, 2018 (in Chinese) [Google Scholar]

[19] ZHU Boling. Research and program development of general trajectory optimization method based on collocation method[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2023 (in Chinese) [Google Scholar]