- 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
65 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
5 Comments
James Tursa
on 13 Oct 2022
Edited: James Tursa
on 13 Oct 2022
"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?
See Also
Categories
Find more on Visualization 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!