GPOPS-II教程(2): 可复用火箭再入大气层最优轨迹规划问题

问题描述

考虑一类可复用火箭再入大气层最优轨迹规划问题，其动力学方程为

{ r ˙ = v sin ⁡ γ , θ ˙ = v cos ⁡ γ sin ⁡ ψ r cos ⁡ ϕ , ϕ ˙ = v cos ⁡ γ cos ⁡ ψ r , v ˙ = − F d m − F g sin ⁡ γ , γ ˙ = F l cos ⁡ σ m v − ( F g v − v r ) cos ⁡ γ , ψ ˙ = F l sin ⁡ σ m v cos ⁡ γ + v cos ⁡ r sin ⁡ ψ 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​​−Fg​sinγ,γ˙​=mvFl​cosσ​−(vFg​​−rv​)cosγ,ψ˙​=mvcosγFl​sinσ​+rvcosrsinψtanϕ​,​

边界条件为

{ r ( 0 ) = 79248 + R e  m  ,   r ( t f ) = 79248 + R e  m , θ ( 0 ) = 0  deg  ,   θ ( t f ) = Free , ϕ ( 0 ) = 0  deg  ,   ϕ ( t f ) = Free , v ( 0 ) = 7803  m/s  ,   v ( t f ) = 762  m/s , γ ( 0 ) = − 1  deg  ,   γ ( t f ) = − 5  deg , ψ ( 0 ) = 90  deg  ,   ψ ( t f ) = 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 = − ϕ ( t f ) . 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个步骤，分别是：

  1. 初始参数设置；
  2. 边界条件设置；
  3. 初值猜测；
  4. 设置GPOPS求解器参数；
  5. 求解；
  6. 画图。

那么，一步一步地来写代码吧。

1. 初始参数设置

%% 01.初始参数设置
%-------------------------------------------------------------------------%
%----------------------- 设置问题的求解边界 ------------------------------%
%-------------------------------------------------------------------------%
cft2m = 0.3048;
cft2km = cft2m/1000;
cslug2kg = 14.5939029;
%-------------------------------------%
%             Problem Setup           %
%-------------------------------------%
auxdata.Re = 20902900*cft2m;              % Equatorial Radius of Earth (m)
auxdata.S  = 2690*cft2m^2;                % Vehicle Reference Area (m^2)
auxdata.cl(1) = -0.2070;                  % Parameters for lift coefficient
auxdata.cl(2) = 1.6756;       
auxdata.cd(1) = 0.0785;       
auxdata.cd(2) = -0.3529;       
auxdata.cd(3) = 2.0400;
auxdata.b(1)  = 0.07854;      
auxdata.b(2)  = -0.061592;    
auxdata.b(3)  = 0.00621408;
auxdata.H     = 23800*cft2m;              % Density Scale Height (m)
auxdata.al(1) = -0.20704;    
auxdata.al(2) = 0.029244;
auxdata.rho0  = 0.002378*cslug2kg/cft2m^3;% Sea Level Atmospheric Density (slug/ft^3)
auxdata.mu    = 1.4076539e16*cft2m^3;     % Earth Gravitational Parameter (ft^^3/s^2) 
auxdata.mass  = 6309.433*cslug2kg;      

% 初始条件
t0 = 0;
alt0 = 260000*cft2m;
rad0 = alt0+auxdata.Re;
lon0 = 0;
lat0 = 0;
speed0 = 25600*cft2m;
fpa0   = -1*pi/180;
azi0   = 90*pi/180;

% 终端条件
altf = 80000*cft2m;
radf = altf+auxdata.Re;
speedf = 2500*cft2m;
fpaf   = -5*pi/180;
azif   = -90*pi/180;

%----------------------------------------------------%
% 时间、状态和控制量的上界和下界
%----------------------------------------------------%
tfMin = 0;            tfMax = 3000;
radMin = auxdata.Re;  radMax = rad0;
lonMin = -pi;         lonMax = -lonMin;
latMin = -70*pi/180;  latMax = -latMin;
speedMin = 10;        speedMax = 45000;
fpaMin = -80*pi/180;  fpaMax =  80*pi/180;
aziMin = -180*pi/180; aziMax =  180*pi/180;
aoaMin = -90*pi/180;  aoaMax = -aoaMin;
bankMin = -90*pi/180; bankMax =   1*pi/180;

2. 边界条件设置

%-------------------------------------------------------------------------%
%------------------------ 将求解边界设置于问题中 -------------------------%
%-------------------------------------------------------------------------%
bounds.phase.initialtime.lower = t0;
bounds.phase.initialtime.upper = t0;
bounds.phase.finaltime.lower = tfMin;
bounds.phase.finaltime.upper = tfMax;
bounds.phase.initialstate.lower = [rad0, lon0, lat0, speed0, fpa0, azi0];
bounds.phase.initialstate.upper = [rad0, lon0, lat0, speed0, fpa0, azi0];
bounds.phase.state.lower = [radMin, lonMin, latMin, speedMin, fpaMin, aziMin];
bounds.phase.state.upper = [radMax, lonMax, latMax, speedMax, fpaMax, aziMax];
bounds.phase.finalstate.lower = [radf, lonMin, latMin, speedf, fpaf, aziMin];
bounds.phase.finalstate.upper = [radf, lonMax, latMax, speedf, fpaf, aziMax];
bounds.phase.control.lower = [aoaMin, bankMin];
bounds.phase.control.upper = [aoaMax, bankMax];

3. 初值猜测

%-------------------------------------------------------------------------%
%------------------------------- 初值猜想 --------------------------------%
%-------------------------------------------------------------------------%
tGuess = [0; 1000];
radGuess = [rad0; radf];
lonGuess = [lon0; lon0+10*pi/180];
latGuess = [lat0; lat0+10*pi/180];
speedGuess = [speed0; speedf];
fpaGuess = [fpa0; fpaf];
aziGuess = [azi0; azif];
aoaGuess = [0; 0];
bankGuess = [0; 0];

guess.phase.state   = [radGuess, lonGuess, latGuess, speedGuess, fpaGuess, aziGuess];
guess.phase.control = [aoaGuess, bankGuess];
guess.phase.time    = tGuess;

4. 设置GPOPS求解器参数

%-------------------------------------------------------------------------%
%---------------------------- 设置求解器参数 -----------------------------%        
%-------------------------------------------------------------------------%
meshphase.colpoints = 4*ones(1,10);
meshphase.fraction = 0.1*ones(1,10);

setup.name = 'Reusable-Launch-Vehicle-Entry-Problem';
setup.functions.continuous = @rlvEntryContinuous;
setup.functions.endpoint   = @rlvEntryEndpoint;
setup.auxdata = auxdata;
setup.mesh.phase = meshphase;
setup.bounds = bounds;
setup.guess = guess;
setup.nlp.solver = 'ipopt';
setup.derivatives.supplier = 'sparseCD';
setup.derivatives.derivativelevel = 'second';
setup.scales.method = 'automatic-bounds';
setup.mesh.method = 'hp1';
setup.mesh.tolerance = 1e-6;
setup.mesh.colpointsmin = 4;
setup.mesh.colpointsmax = 16;

5. 求解

%-------------------------------------------------------------------------%
%----------------------- 使用 GPOPS2 求解最优控制问题 --------------------%
%-------------------------------------------------------------------------%
output = gpops2(setup);
solution = output.result.solution;
toc;

time = solution.phase(1).time;
altitude  = (solution.phase(1).state(:,1)-auxdata.Re)/1000;
longitude = solution.phase(1).state(:,2)*180/pi;
latitude  = solution.phase(1).state(:,3)*180/pi;
speed     = solution.phase(1).state(:,4)/1000;
fpa       = solution.phase(1).state(:,5)*180/pi;
azimuth   = solution.phase(1).state(:,6)*180/pi;
aoa       = solution.phase(1).control(:,1)*180/pi;
bank      = solution.phase(1).control(:,2)*180/pi;

6. 画图

figure('Color',[1,1,1])
pp = plot(time,altitude,'-o', 'markersize', 7, 'linewidth', 1.5);
xl = xlabel('Time (s)');
yl = ylabel('Altitude (km)');
title('Altitude');
set(xl,'FontSize',18);
set(yl,'FontSize',18);
set(gca,'FontSize',16,'FontName','Times New Roman');
set(pp,'LineWidth',1.25);
print -dpng rlvAltitude.png

figure('Color',[1,1,1])
plot(longitude,latitude,'-o', 'markersize', 7, 'linewidth', 1.5);
xl = xlabel('Longitude (deg)');
yl = ylabel('Latitude (deg)');
title('Longitude and Latitude');
set(xl,'FontSize',18);
set(yl,'FontSize',18);
set(gca,'FontSize',16,'FontName','Times New Roman');
set(pp,'LineWidth',1.25);
print -dpng rlvLonLat.png

figure('Color',[1,1,1])
plot(time,speed,'-o', 'markersize', 7, 'linewidth', 1.5);
xl = xlabel('Time (s)');
yl = ylabel('Speed (km/s)');
title('Speed')
set(xl,'FontSize',18);
set(yl,'FontSize',18);
set(gca,'FontSize',16,'FontName','Times New Roman');
set(pp,'LineWidth',1.25);
print -dpng rlvSpeed.png

figure('Color',[1,1,1])
plot(time,fpa,'-o', 'markersize', 7, 'linewidth', 1.5);
yl = xlabel('Time (s)');
xl = ylabel('Flight Path Angle (deg)');
title('Flight Path Angle')
set(xl,'FontSize',18);
set(yl,'FontSize',18);
set(gca,'FontSize',16,'FontName','Times New Roman');
set(pp,'LineWidth',1.25);
print -dpng rlvFlightPathAngle.png

figure('Color',[1,1,1])
plot(time,azimuth,'-o', 'markersize', 7, 'linewidth', 1.5);
yl = xlabel('Time (s)');
xl = ylabel('Azimuth Angle (deg)');
set(xl,'FontSize',18);
set(yl,'FontSize',18);
set(gca,'FontSize',16,'FontName','Times New Roman');
set(pp,'LineWidth',1.25);
print -dpng rlvAzimuthAngle.png

figure('Color',[1,1,1])
plot(time,aoa,'-o', 'markersize', 7, 'linewidth', 1.5);
yl = xlabel('Time (s)');
xl = ylabel('Angle of Attack (deg)');
set(xl,'FontSize',18);
set(yl,'FontSize',18);
set(gca,'FontSize',16,'YTick',[16.5 17 17.5],'FontName','Times New Roman');
set(pp,'LineWidth',1.25);
print -dpng rlvAngleofAttack.png

figure('Color',[1,1,1])
plot(time,bank,'-o', 'markersize', 7, 'linewidth', 1.5);
yl = xlabel('Time (s)');
xl = ylabel('Bank Angle (deg)');
set(xl,'FontSize',18);
set(yl,'FontSize',18);
set(gca,'FontSize',16,'FontName','Times New Roman');
set(pp,'LineWidth',1.25);
print -dpng rlvBankAngle.png

continuous function

写这部分代码的时候分成三步来写。

第一步，把所有要用的变量全部导入进来。

function phaseout = rlvEntryContinuous(input)

rad = input.phase.state(:,1);
lon = input.phase.state(:,2);
lat = input.phase.state(:,3);
speed = input.phase.state(:,4);
fpa = input.phase.state(:,5);
azimuth = input.phase.state(:,6);
aoa = input.phase.control(:,1);
bank = input.phase.control(:,2);

cd0 = input.auxdata.cd(1);
cd1 = input.auxdata.cd(2);
cd2 = input.auxdata.cd(3);
cl0 = input.auxdata.cl(1);
cl1 = input.auxdata.cl(2);
mu  = input.auxdata.mu;
rho0 = input.auxdata.rho0;
H = input.auxdata.H;
S = input.auxdata.S;
mass = input.auxdata.mass;
altitude = rad - input.auxdata.Re;

第二步，计算在求解动力学方程时会用到的变量。

CD = cd0+cd1*aoa+cd2*aoa.^2;

rho = rho0*exp(-altitude/H);
CL = cl0+cl1*aoa;
gravity = mu./rad.^2;
dynamic_pressure = 0.5*rho.*speed.^2;
D = dynamic_pressure.*S.*CD./mass;
L = dynamic_pressure.*S.*CL./mass;
slon = sin(lon);
clon = cos(lon);
slat = sin(lat);
clat = cos(lat);
tlat = tan(lat);
sfpa = sin(fpa);
cfpa = cos(fpa);
sazi = sin(azimuth);
cazi = cos(azimuth);
cbank = cos(bank);
sbank = sin(bank);

第三步，根据动力学方程写出代码。

重新复习一下动力学方程，为

{ r ˙ = v sin ⁡ γ , θ ˙ = v cos ⁡ γ sin ⁡ ψ r cos ⁡ ϕ , ϕ ˙ = v cos ⁡ γ cos ⁡ ψ r , v ˙ = − F d m − F g sin ⁡ γ , γ ˙ = F l cos ⁡ σ m v − ( F g v − v r ) cos ⁡ γ , ψ ˙ = F l sin ⁡ σ m v cos ⁡ γ + v cos ⁡ r sin ⁡ ψ 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​​−Fg​sinγ,γ˙​=mvFl​cosσ​−(vFg​​−rv​)cosγ,ψ˙​=mvcosγFl​sinσ​+rvcosrsinψtanϕ​,​

根据动力学方程，就可以对应地写出，代码如下。

raddot   = speed.*sfpa;
londot   = speed.*cfpa.*sazi./(rad.*clat);
latdot   = speed.*cfpa.*cazi./rad;
speeddot = -D-gravity.*sfpa;
fpadot   = (L.*cbank-cfpa.*(gravity-speed.^2./rad))./speed;
azidot   = (L.*sbank./cfpa + speed.^2.*cfpa.*sazi.*tlat./rad)./speed;

phaseout.dynamics  = [raddot, londot, latdot, speeddot, fpadot, azidot];
end

endpoint function

endpoint function按照性能指标函数写即可，性能指标为

J = − ϕ ( t f ) . J = -\phi(t_f). J=−ϕ(tf​).

那么，代码如下。

% ----------------------------------------------------------------------- %
% ------------------------ BEGIN: rlvEntryEndpoint.m -------------------- %
% ----------------------------------------------------------------------- %
function output = rlvEntryEndpoint(input)
latf = input.phase.finalstate(3);
output.objective = -latf;
end
% ----------------------------------------------------------------------- %
% ------------------------- END: rlvEntryEndpoint.m --------------------- %
% ----------------------------------------------------------------------- %

完整代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 功能描述：可复用火箭再入大气层最优轨迹规划问题
% 文件名解释：mainReentry.m 中，main 代表 主函数，
%             Re-entry 代表 再入航天器
% 作者：Lei Lie
% 时间：2024/06/24
% 版本：1.0
% - 完成了代码框架的初始搭建
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
clear;
close all;
tic;
%% 01.初始参数设置
%-------------------------------------------------------------------------%
%----------------------- 设置问题的求解边界 ------------------------------%
%-------------------------------------------------------------------------%
cft2m = 0.3048;
cft2km = cft2m/1000;
cslug2kg = 14.5939029;
%-------------------------------------%
%             Problem Setup           %
%-------------------------------------%
auxdata.Re = 20902900*cft2m;              % Equatorial Radius of Earth (m)
auxdata.S  = 2690*cft2m^2;                % Vehicle Reference Area (m^2)
auxdata.cl(1) = -0.2070;                  % Parameters for lift coefficient
auxdata.cl(2) = 1.6756;       
auxdata.cd(1) = 0.0785;       
auxdata.cd(2) = -0.3529;       
auxdata.cd(3) = 2.0400;
auxdata.b(1)  = 0.07854;      
auxdata.b(2)  = -0.061592;    
auxdata.b(3)  = 0.00621408;
auxdata.H     = 23800*cft2m;              % Density Scale Height (m)
auxdata.al(1) = -0.20704;    
auxdata.al(2) = 0.029244;
auxdata.rho0  = 0.002378*cslug2kg/cft2m^3;% Sea Level Atmospheric Density (slug/ft^3)
auxdata.mu    = 1.4076539e16*cft2m^3;     % Earth Gravitational Parameter (ft^^3/s^2) 
auxdata.mass  = 6309.433*cslug2kg;      

% 初始条件
t0 = 0;
alt0 = 260000*cft2m;
rad0 = alt0+auxdata.Re;
lon0 = 0;
lat0 = 0;
speed0 = 25600*cft2m;
fpa0   = -1*pi/180;
azi0   = 90*pi/180;

% 终端条件
altf = 80000*cft2m;
radf = altf+auxdata.Re;
speedf = 2500*cft2m;
fpaf   = -5*pi/180;
azif   = -90*pi/180;

%----------------------------------------------------%
% 时间、状态和控制量的上界和下界
%----------------------------------------------------%
tfMin = 0;            tfMax = 3000;
radMin = auxdata.Re;  radMax = rad0;
lonMin = -pi;         lonMax = -lonMin;
latMin = -70*pi/180;  latMax = -latMin;
speedMin = 10;        speedMax = 45000;
fpaMin = -80*pi/180;  fpaMax =  80*pi/180;
aziMin = -180*pi/180; aziMax =  180*pi/180;
aoaMin = -90*pi/180;  aoaMax = -aoaMin;
bankMin = -90*pi/180; bankMax =   1*pi/180;

%% 02.边界条件设置
%-------------------------------------------------------------------------%
%------------------------ 将求解边界设置于问题中 -------------------------%
%-------------------------------------------------------------------------%
bounds.phase.initialtime.lower = t0;
bounds.phase.initialtime.upper = t0;
bounds.phase.finaltime.lower = tfMin;
bounds.phase.finaltime.upper = tfMax;
bounds.phase.initialstate.lower = [rad0, lon0, lat0, speed0, fpa0, azi0];
bounds.phase.initialstate.upper = [rad0, lon0, lat0, speed0, fpa0, azi0];
bounds.phase.state.lower = [radMin, lonMin, latMin, speedMin, fpaMin, aziMin];
bounds.phase.state.upper = [radMax, lonMax, latMax, speedMax, fpaMax, aziMax];
bounds.phase.finalstate.lower = [radf, lonMin, latMin, speedf, fpaf, aziMin];
bounds.phase.finalstate.upper = [radf, lonMax, latMax, speedf, fpaf, aziMax];
bounds.phase.control.lower = [aoaMin, bankMin];
bounds.phase.control.upper = [aoaMax, bankMax];

%% 03.初值猜测
%-------------------------------------------------------------------------%
%------------------------------- 初值猜想 --------------------------------%
%-------------------------------------------------------------------------%
tGuess = [0; 1000];
radGuess = [rad0; radf];
lonGuess = [lon0; lon0+10*pi/180];
latGuess = [lat0; lat0+10*pi/180];
speedGuess = [speed0; speedf];
fpaGuess = [fpa0; fpaf];
aziGuess = [azi0; azif];
aoaGuess = [0; 0];
bankGuess = [0; 0];

guess.phase.state   = [radGuess, lonGuess, latGuess, speedGuess, fpaGuess, aziGuess];
guess.phase.control = [aoaGuess, bankGuess];
guess.phase.time    = tGuess;


%% 04.设置GPOPS求解器参数
%-------------------------------------------------------------------------%
%---------------------------- 设置求解器参数 -----------------------------%        
%-------------------------------------------------------------------------%
meshphase.colpoints = 4*ones(1,10);
meshphase.fraction = 0.1*ones(1,10);

setup.name = 'Reusable-Launch-Vehicle-Entry-Problem';
setup.functions.continuous = @rlvEntryContinuous;
setup.functions.endpoint   = @rlvEntryEndpoint;
setup.auxdata = auxdata;
setup.mesh.phase = meshphase;
setup.bounds = bounds;
setup.guess = guess;
setup.nlp.solver = 'ipopt';
setup.derivatives.supplier = 'sparseCD';
setup.derivatives.derivativelevel = 'second';
setup.scales.method = 'automatic-bounds';
setup.mesh.method = 'hp1';
setup.mesh.tolerance = 1e-6;
setup.mesh.colpointsmin = 4;
setup.mesh.colpointsmax = 16;


%% 05.求解
%-------------------------------------------------------------------------%
%----------------------- 使用 GPOPS2 求解最优控制问题 --------------------%
%-------------------------------------------------------------------------%
output = gpops2(setup);
solution = output.result.solution;
toc;

time = solution.phase(1).time;
altitude  = (solution.phase(1).state(:,1)-auxdata.Re)/1000;
longitude = solution.phase(1).state(:,2)*180/pi;
latitude  = solution.phase(1).state(:,3)*180/pi;
speed     = solution.phase(1).state(:,4)/1000;
fpa       = solution.phase(1).state(:,5)*180/pi;
azimuth   = solution.phase(1).state(:,6)*180/pi;
aoa       = solution.phase(1).control(:,1)*180/pi;
bank      = solution.phase(1).control(:,2)*180/pi;

%% 06.画图

figure('Color',[1,1,1])
pp = plot(time,altitude,'-o', 'markersize', 7, 'linewidth', 1.5);
xl = xlabel('Time (s)');
yl = ylabel('Altitude (km)');
title('Altitude');
set(xl,'FontSize',18);
set(yl,'FontSize',18);
set(gca,'FontSize',16,'FontName','Times New Roman');
set(pp,'LineWidth',1.25);
print -dpng rlvAltitude.png

figure('Color',[1,1,1])
plot(longitude,latitude,'-o', 'markersize', 7, 'linewidth', 1.5);
xl = xlabel('Longitude (deg)');
yl = ylabel('Latitude (deg)');
title('Longitude and Latitude');
set(xl,'FontSize',18);
set(yl,'FontSize',18);
set(gca,'FontSize',16,'FontName','Times New Roman');
set(pp,'LineWidth',1.25);
print -dpng rlvLonLat.png

figure('Color',[1,1,1])
plot(time,speed,'-o', 'markersize', 7, 'linewidth', 1.5);
xl = xlabel('Time (s)');
yl = ylabel('Speed (km/s)');
title('Speed')
set(xl,'FontSize',18);
set(yl,'FontSize',18);
set(gca,'FontSize',16,'FontName','Times New Roman');
set(pp,'LineWidth',1.25);
print -dpng rlvSpeed.png

figure('Color',[1,1,1])
plot(time,fpa,'-o', 'markersize', 7, 'linewidth', 1.5);
yl = xlabel('Time (s)');
xl = ylabel('Flight Path Angle (deg)');
title('Flight Path Angle')
set(xl,'FontSize',18);
set(yl,'FontSize',18);
set(gca,'FontSize',16,'FontName','Times New Roman');
set(pp,'LineWidth',1.25);
print -dpng rlvFlightPathAngle.png

figure('Color',[1,1,1])
plot(time,azimuth,'-o', 'markersize', 7, 'linewidth', 1.5);
yl = xlabel('Time (s)');
xl = ylabel('Azimuth Angle (deg)');
set(xl,'FontSize',18);
set(yl,'FontSize',18);
set(gca,'FontSize',16,'FontName','Times New Roman');
set(pp,'LineWidth',1.25);
print -dpng rlvAzimuthAngle.png

figure('Color',[1,1,1])
plot(time,aoa,'-o', 'markersize', 7, 'linewidth', 1.5);
yl = xlabel('Time (s)');
xl = ylabel('Angle of Attack (deg)');
set(xl,'FontSize',18);
set(yl,'FontSize',18);
set(gca,'FontSize',16,'YTick',[16.5 17 17.5],'FontName','Times New Roman');
set(pp,'LineWidth',1.25);
print -dpng rlvAngleofAttack.png

figure('Color',[1,1,1])
plot(time,bank,'-o', 'markersize', 7, 'linewidth', 1.5);
yl = xlabel('Time (s)');
xl = ylabel('Bank Angle (deg)');
set(xl,'FontSize',18);
set(yl,'FontSize',18);
set(gca,'FontSize',16,'FontName','Times New Roman');
set(pp,'LineWidth',1.25);
print -dpng rlvBankAngle.png

%% 函数模块部分
% ----------------------------------------------------------------------- %
% ----------------------- BEGIN: rlvEntryContinuous.m ------------------- %
% ----------------------------------------------------------------------- %
function phaseout = rlvEntryContinuous(input)

rad = input.phase.state(:,1);
lon = input.phase.state(:,2);
lat = input.phase.state(:,3);
speed = input.phase.state(:,4);
fpa = input.phase.state(:,5);
azimuth = input.phase.state(:,6);
aoa = input.phase.control(:,1);
bank = input.phase.control(:,2);

cd0 = input.auxdata.cd(1);
cd1 = input.auxdata.cd(2);
cd2 = input.auxdata.cd(3);
cl0 = input.auxdata.cl(1);
cl1 = input.auxdata.cl(2);
mu  = input.auxdata.mu;
rho0 = input.auxdata.rho0;
H = input.auxdata.H;
S = input.auxdata.S;
mass = input.auxdata.mass;
altitude = rad - input.auxdata.Re;

CD = cd0+cd1*aoa+cd2*aoa.^2;

rho = rho0*exp(-altitude/H);
CL = cl0+cl1*aoa;
gravity = mu./rad.^2;
dynamic_pressure = 0.5*rho.*speed.^2;
D = dynamic_pressure.*S.*CD./mass;
L = dynamic_pressure.*S.*CL./mass;
slon = sin(lon);
clon = cos(lon);
slat = sin(lat);
clat = cos(lat);
tlat = tan(lat);
sfpa = sin(fpa);
cfpa = cos(fpa);
sazi = sin(azimuth);
cazi = cos(azimuth);
cbank = cos(bank);
sbank = sin(bank);

raddot   = speed.*sfpa;
londot   = speed.*cfpa.*sazi./(rad.*clat);
latdot   = speed.*cfpa.*cazi./rad;
speeddot = -D-gravity.*sfpa;
fpadot   = (L.*cbank-cfpa.*(gravity-speed.^2./rad))./speed;
azidot   = (L.*sbank./cfpa + speed.^2.*cfpa.*sazi.*tlat./rad)./speed;

phaseout.dynamics  = [raddot, londot, latdot, speeddot, fpadot, azidot];
end
% ----------------------------------------------------------------------- %
% ------------------------ END: rlvEntryContinuous.m -------------------- %
% ----------------------------------------------------------------------- %

% ----------------------------------------------------------------------- %
% ------------------------ BEGIN: rlvEntryEndpoint.m -------------------- %
% ----------------------------------------------------------------------- %
function output = rlvEntryEndpoint(input)

latf = input.phase.finalstate(3);
output.objective = -latf;
end
% ----------------------------------------------------------------------- %
% ------------------------- END: rlvEntryEndpoint.m --------------------- %
% ----------------------------------------------------------------------- %

代码仿真结果

高度：
在这里插入图片描述

速度：

在这里插入图片描述

经纬度：

在这里插入图片描述

飞行路径角度：

在这里插入图片描述

攻角：

在这里插入图片描述

倾斜角：

在这里插入图片描述

最后

关于可复用火箭再入大气层的最优轨迹规划问题，没有写得很详细每个参数为什么这么取的原因，是我觉得现在大家都大概率不需要解决这个问题。所以诸如攻角、仰角、迎角之类的变量不去解释为什么要这么写。

只是希望通过这样一系列的例子告诉大家GPOPS-II应该怎么使用的思路。

  • 动力学方程应该怎么写？
  • 性能指标应该怎么写？
  • 约束应该怎么写？
  • 猜测应该怎么写？
  • 画图怎么画？

要被解决的问题是无限的，只有掌握了方法论，才能一法通时万法通。

欢迎通过邮箱联系我：lordofdapanji@foxmail.com

来信请注明你的身份，否则恕不回信。

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.mfbz.cn/a/744500.html

如若内容造成侵权/违法违规/事实不符，请联系我们进行投诉反馈qq邮箱809451989@qq.com，一经查实，立即删除！

相关文章

解决chrome浏览器总是将对站点的http访问改为https的问题

解决chrome浏览器总是将对站点的http访问改为https的问题

问题：vue项目本地运行出来的地址是http开头的，但在chrome浏览器中无法访问，在Edge浏览器就可以，发现是chrome总是自动将http协议升级为https。 已试过的有效的方法： 一、无痕模式下访问 无痕模式下访问不会将http自…
阅读更多...
推送电子邮箱与其他营销手段如何有效结合？

推送电子邮箱与其他营销手段如何有效结合？

推送电子邮箱的效果如何？怎么优化邮件推送的策略？ 将推送电子邮箱与其他营销手段有效结合，可以显著提升营销效果和用户体验。AokSend将探讨如何将推送电子邮箱与社交媒体营销、内容营销、搜索引擎优化（SEO）等手段相结…
阅读更多...
Python 实现Excel转TXT，或TXT文本导入Excel

Python 实现Excel转TXT，或TXT文本导入Excel

Excel是一种具有强大的数据处理和图表制作功能的电子表格文件，而TXT则是一种简单通用、易于编辑的纯文本文件。将Excel转换为TXT可以帮助我们将复杂的数据表格以文本的形式保存，方便其他程序读取和处理。而将TXT转换为Excel则可以将文本文件中的数据导入…
阅读更多...
鸿蒙应用开发 - 软件安装 - DevEco

鸿蒙应用开发 - 软件安装 - DevEco

第一步 前往下载点下载安装包下载中心 | 华为开发者联盟-HarmonyOS开发者官网，共建鸿蒙生态https://developer.huawei.com/consumer/cn/download/ 根据自身需求下载对应安装包 第二步 点击打开安装包,配置安装路径 我个人选择放E盘,避免占用c盘空间 第三步 …
阅读更多...
Linux集群自动化维护-Ansible

Linux集群自动化维护-Ansible

1.1Ansible概述 自动化运维：批量管理，批量分发，批量执行，维护。。是python写的 批量管理工具： Ansible（无客户端）：无客户端，基于ssh进行管理与维护 Saltstack &#…
阅读更多...
RabbitMQ中lazyqueue队列

RabbitMQ中lazyqueue队列

lazyqueue队列非常强悍 springboot注解方式开启 // 使用注解的方式lazy.queue队列模式 非常GoodRabbitListener(queuesToDeclare Queue(name "lazy.queue",durable "true",arguments Argument(name "x-queue-mode",value "lazy&…
阅读更多...
RocketMQ源码学习笔记：Broker启动流程

RocketMQ源码学习笔记：Broker启动流程

这是本人学习的总结，主要学习资料如下 马士兵教育rocketMq官方文档 目录 1、Broker启动流程2、一些重要的类2.1、MappedFile2.2、MessgeStore2.3、MessageStore的加载启动流程 3、技术亮点3.1、 内存映射3.1.1、简介3.1.2、源码 1、Broker启动流程 Broker启动流程…
阅读更多...
upload-labs第14关

upload-labs第14关

upload-labs第14关 第十四关一、源代码分析代码审计 二、绕过分析a. 制作图片码首先需要一个照片，然后其次需要一个eval.php。 b.上传图片码上传成功 c.结合文件包含漏洞进行访问访问：http://192.168.1.110/upload-labs-master/include.php?filehttp://…
阅读更多...
【Spine学习16】之 人物面部绑定

【Spine学习16】之 人物面部绑定

1、创建头部骨骼 一根头骨 以头骨为父结点创建一个面部控制器face-holder 2、创建头发和face面部控制结点的变换约束 左右头发的约束指向为face结点 3、设定后发的变换约束，约束指向为face结点，反方向移动 设置参数为-100 同理，耳朵也依…
阅读更多...
6.25C高级

6.25C高级

终端输入两个数&#xff0c;判断两数是否相等&#xff0c;如果不相等&#xff0c;判断大小关系 #!/bin/bash if [ $1 -eq $2 ] thenecho $1$2 elif [ $1 -gt $2 ] thenecho "$1>$2" elseecho "$1<$2" fi 2.已知网址www.hqyj.com&#xff0c;使用e…
阅读更多...
Java-day01--基础知识

Java-day01--基础知识

1、计算机基础知识&#xff1a; 计算机主要是由硬件和软件组成&#xff0c;软件指的是特定顺序的计算机指令&#xff0c;硬件主要可以看成是系统软件和应用软件等。 目前java主流分成三种&#xff1a;javase&#xff08;标准版&#xff09;、javame&#xff08;小型版&#x…
阅读更多...
TCP: 传输控制协议

TCP: 传输控制协议

TCP: 传输控制协议 TCP的服务TCP 的首部小结 本系列文章旨在巩固网络编程理论知识&#xff0c;后续将结合实际开展深入理解的文章。 TCP的服务 T C P和U D P都使用相同的网络层&#xff08;I P&#xff09;&#xff0c;T C P却向应用层提供与U D P完全不同的服务。 T C P提供一…
阅读更多...
数据结构-----【链表:刷题】

数据结构-----【链表:刷题】

-------------------------------------------基础题参照leetcode---------------------------------------------------------------------------------------------------------- 【2】两数相加 /*** Definition for singly-linked list.* struct ListNode {* int val;…
阅读更多...
Midjourney 出图效果总比别人差？看看是不是这3点没做好！

Midjourney 出图效果总比别人差？看看是不是这3点没做好！

前言 大家好&#xff0c;我是AI绘画咪酱~ 很多小伙伴都会在「优设 AI 绘画交流群」以及文章评论区中非常积极地和我探讨 Midjourney 的使用问题&#xff0c;我在帮助大家的同时自己也学到了不少新的技巧。今天就写一篇文章解释 3 个比较常见的疑问&#xff0c;希望对大家有帮…
阅读更多...
vue项目无后台版本打包上传到服务器

vue项目无后台版本打包上传到服务器

打包项目 也可以在文件目录下npm run build 生成dist文件夹 将dist文件夹里的所有文件拷贝到站点的根目录&#xff0c;这里使用宝塔面板进行操作 前提你得先创建站点&#xff0c;域名绑定等操作
阅读更多...
IIC学习笔记（立创STMF4开发板）

IIC学习笔记（立创STMF4开发板）

目录 #I2C涉及相关知识 #I2C相关介绍 欢迎指正&#xff0c;希望对你&#xff0c;有所帮助&#xff01;&#xff01;&#xff01; 个人学习笔记&#xff0c;参考文献&#xff0c;链接最后&#xff01;&#xff01;&#xff01; #I2C涉及相关知识 SDA串行数据线&#xff1a; Ser…
阅读更多...
win11误删生物识别设备

win11误删生物识别设备

不小心给删掉了&#xff0c;我当时删的时候没有打勾&#xff0c;解决方法如下&#xff1a; ************ 很简单&#xff0c; 1.点一下 设置 然后找到下面的界面&#xff0c; 2.找到圈出来的功能 3.再找 “可选更新”的功能点进去 4.一下界面可能有一个能勾选&#xff0c;也可…
阅读更多...
怎么样才能让老旧的和颜色受损、丢失的照片重新上色呢？

怎么样才能让老旧的和颜色受损、丢失的照片重新上色呢？

怎么样才能让老旧的和颜色受损、丢失的照片重新上色呢&#xff1f;大家有时候在家中打扫卫生的时候&#xff0c;偶然发现了自己爸爸妈妈以前拍的照片&#xff0c;但是照片颜色已经受损的很严重了&#xff0c;几乎就是黑白的颜色&#xff0c;很难看清楚爸爸妈妈年轻时候的样子&a…
阅读更多...
泵设备的监测控制和智慧运维

泵设备的监测控制和智慧运维

泵是一种输送流体或使流体增压的机械。它通过各种工作原理&#xff08;如离心、柱塞等&#xff09;将机械能转换为流体的动能或压力能&#xff0c;从而实现液体的输送、提升、循环等操作。 泵的一些具体应用场景&#xff1a; 1.智能水务&#xff1a;在城市供水管网中&#xff…
阅读更多...
Lesson 39 Don‘t drop it!

Lesson 39 Don‘t drop it!

Lesson 39 Don’t drop it! 词汇 front n. 前面 搭配&#xff1a;in front of … 在……前面&#xff08;外部&#xff09; 　　　in the front of … 在……前面&#xff08;内部&#xff09; 例句&#xff1a;Bobby坐在Sam的前面。 　　　Bobby is sitting in front of Sam…
阅读更多...
最新文章

Copyright @ 2022~2023

友情链接： 我的编程经验分享 一条长河网 尧图网站开发 尧图建网站