Optimization of Rocket Launch

121 次查看（过去 30 天）
Hamish Brown
Hamish Brown 2022-10-12
评论： Hamish Brown 2022-10-13
I've written a matlab code that simulates the first 500 seconds of a stageless rocket's flight. At the end, it reaches horizontal flight at an altitude of around 350km and velocity of about 7.5km/s. I used Euler integration in a 'for loop' to solve the equations of motion. Now, I'd like to optimise the problem by seeing which trajectory maximises the payload into a 350km orbit with a final velocity of about 7.5km/s. The optimisation variables could be the mass and ISP of the fuel, thrust, the time at which the pitch over manoeuvre is started at and the angle at which the pitchover occurs. I've never done an optimisation problem in MATlab before so am unsure of where to start. Here is the solver part of the code. Thanks for any help.
mass_payload = 150 ; %mass of payload (kg)
mass_fuel_struc = 2500; %mass of the rocket structure including fuel (kg)
m0(i) = mass_payload + mass_fuel_struc; %total rocket mass (kg)
R_E = 6.3781*10^6; %radius of the Earth (m)
M_E = 5.9722*10^24; %mass of the Eath (kg)
G = 6.674*10^-11; %universal gravitational constant (m^3 kg^-1 s^-2)
thrust = 30000 ; %initial engine thrust (N)
isp = 292; %ISP of fuel (s)
orbit = 350000; %target orbit altitude (m)
sim_time = 600; time = 0:dt:sim_time; %total simulation time
length_sim = length(time);
mdot = mass_fuel_struc/burntime; %mass fuel consumption rate (kg/s)
for i = 1:length_sim-1
%density calculator
if y<=11000
Temp(i+1) = 15.04-0.00649*y(i);
p(i+1) = 101.29*((Temp(i)+273.1)/288.08)^5.258;
elseif y <=25000
Temp(i+1) = -56.46;
p(i+1) = 22.65*exp(1.73-0.000157*y(i));
else
Temp(i+1) = -131.21 + 0.00299*y(i);
p(i+1) = 2.488*((Temp(i)+273.1)/216.6)^-11.388;
end
rho(i+1) = p(i)/(0.2869*(Temp(i)+273.1));
% gravity model
g(i+1) = (G*M_E)/(R_E+y(i))^2;
%drag force D
D(i+1) = 0.5*rho(i)*v(i)^2*Area*cd;
%speed magnitude v
v(i+1) = v(i) + (T/m0(i) - D(i)/m0(i) - g(i)*sin(gam(i)))*dt;
%path angle gam (from 90 to 0 degrees)
if i<=25/dt
gam(i+1) = gam(i);
elseif 25/dt< i <= 40/dt
gam(i+1) = gam(i) - 0.0085*gam(i)*dt;
else
gam(i+1) = gam(i) - ((g(i)/v(i)-v(i)/(R_E+y (i)))*cos(gam(i)))*dt;
end
%x distance downrange
x(i+1) = x(i) + (R_E/(R_E+y(i))*v(i)*cos(gam(i)))*dt;
%altitude y
y(i+1) = y(i) + (v(i)*sin(gam(i)))*dt;
VG(i) = g(i)*sin(gam(i));
VD(i) = D(i)/m0(i);
% rocket mass
if i <= burntime/dt
m0(i+1) = m0(i) - mdot*dt;
else %rocket has run out of fuel so mass is now constant
m0(i+1) = m0(round(burntime/dt));
T = 0;
end
%dynamic pressure
q(i+1) = (rho(i) * v(i)^2)/2;
%Vertical Velocity
Vy(i+1) = v(i)*sin(gam(i));
%Horizontal Velocity
Vx(i+1) = v(i)*cos(gam(i));
end
  1 个评论
Les Beckham
Les Beckham 2022-10-12
You will have to get your simulation working before you can do any optimization. The posted code does not work.
A few additional thoughts.
  1. I don't see anything in this code that controls the final altitude of the trajectory to achieve the desired 350 km so you are probably going to have to run a separate optimization on your path angle algorithm for each combination of the other parameters to achieve the right final altitude.
  2. Generally you wouldn't have the ability to freely adjust the quantity of fuel, so you might want to eliminate from your optimized parameters, or at least put some limits on how much it can adjust.
  3. The same is true of thrust and ISP.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

James Tursa
James Tursa 2022-10-13
编辑：James Tursa 2022-10-13
You have proposed work that IMO is way beyond what is achievable with your code. Constrained optimization of this sort is not trivial. Trying to maximize payload subject to final altitude and velocity constraints will heavily depend on flight path, which your program isn't even set up to adjust for optimization as far as I can tell. My suggestion is to back off your goals and try something much simpler. Like get your program working properly for always hitting the correct altitude at horizontal flight, i.e. constrain the end point. Then start to work on the optimization stuff later.
It looks like you are integrating x, y, speed, and flight path angle directly. This is a bit unusual ... maybe you could post the differential equations you are using for this so we can compare it with your code.
Finally, this line doesn't do what you expect:
elseif 25/dt< i <= 40/dt
It should be this instead
elseif i <= 40/dt
  5 个评论
显示 3更早的评论
James Tursa
James Tursa 2022-10-13
编辑：James Tursa 2022-10-13
"I've modelled the Earth as a flat surface"
Not sure about this. Why would there be a RE/(RE+h) factor in the downrange equation if you modelled the surface as flat? It seems to me that the RE/(RE+h) factor is to acccount for a downrange measurement along the curved surface of a planet with radius RE at arbitrary height h. That is, the ratio of a circle circumference at the surface divided by a circle circumference at altitude h is RE/(RE+h), and downrange x is measured along the surface at radius RE. At least that is what it looks like is going on in this equation. Can you elaborate on this?
Hamish Brown
Hamish Brown 2022-10-13
yes, you're right. The equations of motion make the assumptions that the Earth is spherical and non-rotating. As such, the coriolis effect is ignored and centrifugal forces.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Geoscience 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by