问题描述
考虑一类可复用火箭再入大气层最优轨迹规划问题,其动力学方程为
{ r˙=vsinγ,θ˙=vcosγsinψrcosϕ,ϕ˙=vcosγcosψr,v˙=−Fdm−Fgsinγ,γ˙=Flcosσmv−(Fgv−vr)cosγ,ψ˙=Flsinσmvcosγ+vcosrsinψtanϕr, \left \{ \begin{aligned} &\dot r = v \sin \gamma, \\ &\dot \theta = \frac{v \cos \gamma \sin \psi}{r \cos \phi}, \\ &\dot \phi = \frac{v \cos \gamma \cos \psi}{r}, \\ &\dot v = - \frac{F_d}{m}-F_g \sin \gamma, \\ &\dot \gamma = \frac{F_l \cos \sigma}{m v} - (\frac{F_g}{v} - \frac{v}{r}) \cos \gamma, \\ &\dot \psi = \frac{F_l \sin \sigma}{m v \cos \gamma} + \frac{v \cos r \sin \psi \tan \phi}{r}, \end{aligned} \right . ⎩ ⎨ ⎧r˙=vsinγ,θ˙=rcosϕvcosγsinψ,ϕ˙=rvcosγcosψ,v˙=−mFd−Fgsinγ,γ˙=mvFlcosσ−(vFg−rv)cosγ,ψ˙=mvcosγFlsinσ+rvcosrsinψtanϕ,
边界条件为
{ r(0)=79248+Re m , r(tf)=79248+Re m,θ(0)=0 deg , θ(tf)=Free,ϕ(0)=0 deg , ϕ(tf)=Free,v(0)=7803 m/s , v(tf)=762 m/s,γ(0)=−1 deg , γ(tf)=−5 deg,ψ(0)=90 deg , ψ(tf)=Free. \left \{ \begin{array}{lcl} r(0) = 79248 + R_e \ \text{m}\ &,\ &r(t_f) = 79248 + R_e \ \text{m}, \\ \theta(0) = 0 \ \text{deg}\ &,\ & \theta(t_f) = \text{Free}, \\ \phi(0) = 0 \ \text{deg}\ &,\ & \phi(t_f) = \text{Free}, \\ v(0) = 7803 \ \text{m/s}\ &,\ & v(t_f) = 762 \ \text{m/s}, \\ \gamma(0) = -1 \ \text{deg}\ &,\ & \gamma(t_f) = -5 \ \text{deg}, \\ \psi(0) = 90 \ \text{deg}\ &,\ & \psi(t_f) = \text{Free}. \end{array} \right. ⎩ ⎨ ⎧r(0)=79248+Re m θ(0)=0 deg ϕ(0)=0 deg v(0)=7803 m/s γ(0)=−1 deg ψ(0)=90 deg , , , , , , r(tf)=79248+Re m,θ(tf)=Free,ϕ(tf)=Free,v(tf)=762 m/s,γ(tf)=−5 deg,ψ(tf)=Free.
性能指标为
J=−ϕ(tf). J = -\phi(t_f). J=−ϕ(tf).
参考文献: [1] Betts J T. Practical methods for optimal control and estimation using nonlinear programming[M]. Society for Industrial and Applied Mathematics, 247-252, 2010.
GPOPS代码
main function
虽然这个最优控制问题很复杂,不过不要着急,心里要有一个顺序,按照顺序一步一步写下去就行。
按照我的习惯,main function
一般分成6个步骤,分别是:
- 初始参数设置;
- 边界条件设置;
- 初值猜测;
- 设置GPOPS求解器参数;
- 求解;
- 画图。
那么,一步一步地来写代码吧。
1. 初始参数设置
%% 01.初始参数设置
%-------------------------------------------------------------------------%
%----------------------- 设置问题的求解边界 ------------------------------%
%-------------------------------------------------------------------------%
cft2m = 0.3048;
cft2km = cft2m/1000;
cslug2kg = 14.5939029;
%-------------------------------------%
% Problem Setup %
%-------------------------------------%
auxdata.Re