Optimization of Rocket Launch

111 views (last 30 days)
Hamish Brown
Hamish Brown on 12 Oct 2022
Commented: Hamish Brown on 13 Oct 2022
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));
Temp(i+1) = -131.21 + 0.00299*y(i);
p(i+1) = 2.488*((Temp(i)+273.1)/216.6)^-11.388;
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;
gam(i+1) = gam(i) - ((g(i)/v(i)-v(i)/(R_E+y (i)))*cos(gam(i)))*dt;
%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;
%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));
  1 Comment
Les Beckham
Les Beckham on 12 Oct 2022
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.

Sign in to comment.

Answers (1)

James Tursa
James Tursa on 13 Oct 2022
Edited: James Tursa on 13 Oct 2022
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
Hamish Brown
Hamish Brown on 13 Oct 2022
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.

Sign in to comment.


Find more on Earth and Planetary Science in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by