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

news2024/10/7 8:29:00

问题描述

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

{ 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.coloradmin.cn/o/1862229.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

malloc和new的本质区别

malloc和new的本质区别

目录 一、结论 二、示例 1.实现类T 2.用malloc分配类T的内存空间 3.用new分配类T的内存空间 一、结论 malloc 和 new 都是用于在运行时动态分配内存的机制。但它们之间存在一些本质的区别，主要是在使用方面，现在我们直接说结论，然后在通过…
阅读更多...
日光模拟器对显示器光干涉影响

日光模拟器对显示器光干涉影响

太阳光模拟器应用领域 抬头显示器TFT日光照射仿真太阳光模拟器LED显示器阳光耐候老化测试仪器中的光照实验抬头显示器TFT日光照射仿真太阳光模拟器全光谱阳光太阳辐射环境模拟系统中的光LCD和OLED显示器强光试验太阳光模拟器日光模拟器的光谱匹配度测量方法新型LED太阳模拟器光…
阅读更多...
Android音频系统

Android音频系统

最近在做UAC的项目，大概就是接收内核UAC的事件，也就是声音相关事件。然后就是pcm_read和AudioTrackr->write之间互传。感觉略微有点奇怪，所以简单总结一下。 1 UAC的简要流程 open_netlink_socket 打开内核窗口，类似于ioctl。…
阅读更多...
Verilog的逻辑系统及数据类型（一）:四值逻辑系统

Verilog的逻辑系统及数据类型（一）:四值逻辑系统

目录 1. Verilog采用的四值逻辑系统2.主要数据类型2.1 net（线网）2.2 寄存器类 （register)2.3 Verilog中net和register声明语法2.3.1 net声明2.3.2 寄存器声明 2.4 选择正确的数据类型2.5 选择数据类型时常犯的错误2.5.1 信号类型确定方法总结…
阅读更多...
使用 Spring Boot 3.x 与图形学技术，添加电子印章防伪特征

使用 Spring Boot 3.x 与图形学技术，添加电子印章防伪特征

使用 Spring Boot 3.x 与图形学技术，添加电子印章防伪特征 在电子办公和无纸化办公日益普及的今天，电子印章的使用越来越广泛。然而，如何确保电子印章的安全性和防伪能力成为了一个亟待解决的问题。本文将通过 Spring Boot 3.x 和图形学技术，深入探讨如何为电子印章添加防…
阅读更多...
速卖通自养号测评：安全高效的推广手段

速卖通自养号测评：安全高效的推广手段

在速卖通平台上，卖家们常常寻求各种方法来提升商品的曝光、转化率和店铺权重。其中，自养号测评作为一种低成本、高回报的推广方式，备受关注。然而，若操作不当，也可能带来风险。以下是如何安全有效地进行自养号测评的指…
阅读更多...
label studio数据标注平台的自动化标注使用

label studio数据标注平台的自动化标注使用

（作者：陈玓玏） 开源项目，欢迎star哦，https://github.com/tencentmusic/cube-studio 做图文音项目过程中，我们通常会需要进行数据标注。label studio是一个比较好上手的标注平台，可以直接搜…
阅读更多...
MAB规范（3）：Chapter6 Glossary 术语表

MAB规范（3）：Chapter6 Glossary 术语表

第6章 - 术语表 此章不做过多的批注，都是些简单的术语解释。
阅读更多...
【算法学习】判断点在多边形内外的算法以及确定内外两点连线与边界的交点

【算法学习】判断点在多边形内外的算法以及确定内外两点连线与边界的交点

1.前言： 在GIS开发中，经常会遇到确定一个坐标点是否在一块区域的内部这一问题。 如果这个问题不是一个单纯的数学问题，例如：在判断DEM、二维图像像素点、3D点云点等含有自身特征信息的这些点是否在一个区域范围内部的时候&#x…
阅读更多...
Java三层框架的解析

Java三层框架的解析

引言：欢迎各位点击收看本篇博客，在历经很多的艰辛，我也是成功由小白浅浅进入了入门行列，也是收货到很多的知识，每次看黑马的JavaWeb课程视频，才使一个小菜鸡见识到了Java前后端是如何进行交互访问的&#x…
阅读更多...
20240626 每日AI必读资讯

20240626 每日AI必读资讯

🌍警告！OpenAI宣布全面封锁中国API接入！ - 7月9号开始封锁不支持的国家API - 如果在OpenAI不允许的国家使用其 API 将面临封杀 🔗 警告！OpenAI 宣布全面封锁中国 API 接入-CSDN博客 🎵索尼、环球音乐、华…
阅读更多...
29-Matplotlib数学表达式

29-Matplotlib数学表达式

Matplotlib数学表达式 Matplotlib 中的文本字符串都可以使用 Text Markup（一种文本标记语言）显现出来，具体的使用方法是将文本标记符放在一对美元符号$内，语法格式如下： #数学表达式 plt.title(r$\alpha > \beta$…
阅读更多...
【c语言】二级指针

【c语言】二级指针

1，定义 本质还是从指针的角度去理解，只不过存的指针的值 2，使用方法
阅读更多...
小程序的基本使用

小程序的基本使用

【 0 】前言 【 0 】 这个就是js代码的存放地方 app.json // pages/banner/banner.js Page({/*** 页面的初始数据*/data: {},/*** 生命周期函数--监听页面加载*/onLoad(options) {},/*** 生命周期函数--监听页面初次渲染完成*/onReady() {},/*** 生命周期函数--监听页面显示…
阅读更多...
互联网应用主流框架整合之Spring Boot运维体系

互联网应用主流框架整合之Spring Boot运维体系

先准备个简单的系统，配置和代码如下 # 服务器配置 server:# 服务器端口port: 8001# Spring Boot 配置 spring:# MVC 配置mvc:# Servlet 配置servlet:# Servlet 的访问路径path: /sbd# 应用程序配置application:# 应用程序名称name: SpringBootDeployment# 配置数据…
阅读更多...
云动态摘要 2024-06-25

云动态摘要 2024-06-25

给您带来云厂商的最新动态，最新产品资讯和最新优惠更新。 最新产品更新 Web应用防火墙 - 验证码支持微信小程序接入 阿里云 2024-06-25 支持客户从微信小程序场景下接入，提供人机识别的安全防护。 工业数字模型驱动引擎 - iDME控制台换新升级 华为云…
阅读更多...
cropperjs 裁剪/框选图片

cropperjs 裁剪/框选图片

1.效果 2.使用组件 <!-- 父级 --><Cropper ref"cropperRef" :imgUrl"url" searchImg"searchImg"></Cropper>3.封装组件 <template><el-dialog :title"title" :visible.sync"dialogVisible" wi…
阅读更多...
亚马逊云科技官方活动：一个月拿下助理架构师SAA＋云从业者考试认证（送半价折扣券）

亚马逊云科技官方活动：一个月拿下助理架构师SAA＋云从业者考试认证（送半价折扣券）

为了帮助大家考取AWS SAA和AWS云从业者认证&#xff0c;小李哥争取到了大量考试半价50%折扣券&#xff0c;使用折扣券考试最多可省75刀(545元人民币)。 领取折扣券需要加入云师兄必过班群&#xff0c;在群中免费领取。目前必过班群招募到了超过200名小伙伴&#xff0c;名额有限…
阅读更多...
前端自动化

前端自动化

前端自动化的内容 自动化代码检查自动化测试自动化构建自动化部署自动化文档 前端自动化的最佳实践
阅读更多...
（六）使用统计学方法进行变量有效性测试（43道选择题）

（六）使用统计学方法进行变量有效性测试（43道选择题）

本文整理了使用统计学方法进行变量有效性测试相关的练习题&#xff0c;共43道&#xff0c;适用于想巩固理论基础的同学。来源&#xff1a;如荷学数据科学题库&#xff08;CDA二级-第7章&#xff09;。 1&#xff09; 2&#xff09; 3&#xff09; 4&#xff09; 5&#xff09;…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023