Optimization of Rocket Launch

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

[Les Beckham]

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.

Answer 1

[James Tursa]

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]

Thanks for the pointers. The equations I'm using are:

Unfortunately my gravity turn doesn't seem to be working properly. After the kickover manoeuvre, the gravity turn equation (dγ/dt) doesn't pull the rocket over enough to achieve horizontal flight. I'll need to figure out how I can do that now.

[James Tursa]

Thanks. Can you post the entire problem statement also? I need to understand the context of these equations and be sure what the other variables shown are in this context. E.g., is v really velocity (a 2-element or 3-element vector) as written, or is it speed (a scalar) as you have stated.

[Hamish Brown]

There isn't a set problem statement as such, but I will try and summarise what I'm doing as best I can.

At this stage, I'm looking at simulating the flight of a 'stageless' rocket which consumes its structure as its fuel. I've based some initial conditions off of the Miura 1 launcher (thrust, mass etc) and an end payload of 100kg. At the end of flight, the rocket's mass will be equal to the payload. I'm wanting the rocket to reach an altitude of around 350km (LEO) with a end horizontal velocity of about 7.5 km/s. The velocity (speed) is indeed a scalar, but I split it into horizontal and vertical components. I'll attach the code I have. I've modelled the Earth as a flat surface.

[James Tursa]

"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]

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.