轨迹优化入门必备 — GPOPS-II进阶

1. 简介

在上一篇文章中，介绍了GPOPS-II的基本用法和入门示例

轨迹优化入门必备 — GPOPS-II

这篇文章将更进一步，介绍GPOPS-II的进阶使用，侧重多阶段轨迹优化、初始轨迹与初始网格设置的工程经验

2. 多阶段轨迹优化

在航天任务、再入任务、以及飞行器的分级发射过程中，最优控制问题往往天然具有多阶段特性

此处以Launch Vehicle Ascent Problem(火箭发射问题)为例

该问题包含四个阶段：

• 第一阶段：助推器工作
• 第二阶段：助推器分离
• 第三阶段：一级火箭分离
• 第四阶段：二级入轨

目标函数：$J=-m\left(t_f\right)$

运动学方程：$\begin{array}{l}\dot{\mathbf{r}}\left(t\right)=\mathbf{v}\left(t\right)\\ \dot{\mathbf{v}}\left(t\right)=\frac{T\mathbf{u}+\mathbf{D}}{m}+\mathbf{g}\\ \dot{m}\left(t\right)=-\frac{T}{g_0{I}_{sp}}\end{array}$

路径约束：${\Vert \mathbf{u}\left(t\right)\Vert }_2=1$

阶段连接约束：$\begin{array}{l}\mathbf{r}\left(t_f^{\left(i\right)}\right)=\mathbf{r}\left(t_0^{\left(i+1\right)}\right)\\ \mathbf{v}\left(t_f^{\left(i\right)}\right)=\mathbf{v}\left(t_0^{\left(i+1\right)}\right)\\ m\left(t_f^{\left(i\right)}\right)=m\left(t_0^{\left(i+1\right)}\right)\end{array},i=1,2,3$

边界条件：$\begin{array}{l}\mathbf{r}\left(0\right)={\mathbf{r}}_0,\mathbf{v}\left(0\right)={\mathbf{v}}_0\\ a\left(t_f\right)=a_f,e\left(t_f\right)=e_f,i\left(t_f\right)=i_f,\Omega \left(t_f\right)={\Omega }_f,\omega \left(t_f\right)={\omega }_f\end{array}$

在 GPOPS-II 中，每个阶段都可以拥有自己的时间范围、状态变量、控制变量以及动力学方程

阶段之间的状态需要通过事件约束(eventgroup)实现连续性

2.1 主函数：launchMain.m

阶段时间固定：通过将下界finaltime.lower和上界finaltime.upper设定为同一固定值

终端时间自由：通过设置不同的上下界，让优化器自动决定阶段时间

归一化处理：为了提高收敛性，GPOPS-II通常要求对问题进行缩放

clear all; clc

% 物理参数
earthRadius         = 6378145;
gravParam           = 3.986012e14;
initialMass         = 301454;
earthRotRate        = 7.29211585e-5;
seaLevelDensity     = 1.225;
densityScaleHeight  = 7200;
g0                  = 9.80665;

% 归一化尺度
scales.length       = earthRadius;
scales.speed        = sqrt(gravParam/scales.length);
scales.time         = scales.length/scales.speed;
scales.acceleration = scales.speed/scales.time;
scales.mass         = initialMass;
scales.force        = scales.mass*scales.acceleration;
scales.area         = scales.length^2;
scales.volume       = scales.area.*scales.length;
scales.density      = scales.mass/scales.volume;
scales.gravparam    = scales.acceleration*scales.length^2;

omega               = earthRotRate*scales.time;
auxdata.omegaMatrix = omega*[0 -1 0;1 0 0;0 0 0];
auxdata.mu          = gravParam/scales.gravparam;
auxdata.cd          = 0.5;
auxdata.sa          = 4*pi/scales.area;
auxdata.rho0        = seaLevelDensity/scales.density;
auxdata.H           = densityScaleHeight/scales.length;
auxdata.Re          = earthRadius/scales.length;
auxdata.g0          = g0/scales.acceleration;

lat0                = 28.5*pi/180;         
x0                  = auxdata.Re*cos(lat0);
y0                  = 0;                   
z0                  = auxdata.Re*sin(lat0);
r0                  = [x0 y0 z0];