- 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.
- 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.
- The same is true of thrust and ISP.

# Optimization of Rocket Launch

111 views (last 30 days)

Show older comments

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 Comment

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.

### Answers (1)

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

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!