Differential Equation Parameter Optimization with Free Time Parameter

I am trying to optimize 3 parameters in a system of differential equations using fminunc. The three parameters are theta0, thetaf and tcutoff, which model the initial thrust angle, final thrust angle and engine cut-off time respectively, of a rocket that is launched from the Moon's surface and ends up in a circular 100km orbit around the Moon when its engine is turned off. Below are the functions and main file that perform the optimization and running of the simulation output:

The first file is main.m, which performs the parameter optimization and runs the simulation once the parameters have been found by fminunc:

clear
clc
close all
format long g

%Near optimal solution parameters found using Mathematica's FindRoot function
theta0 = 1.15055; %Initial thrust angle (radians)
thetaf = -0.867385; %Final thrust angle (radians)
tcutoff = 430.086; %Engine cutoff time (s)

%Parameter optimization
parameters_init = [0, 0, 450] %Initial parameter guess values
parameters = fminunc(@cost_function, parameters_init)

%%Running simulation with parameter outputs
G = 6.672e-11; %Gravitational constant
M = 7.3476e22; %Mass of Moon (kg)
R = 1.737e6; %Radius of Moon (m)
g = 9.81; %Gravitational acceleration on Earth (m/s^2);
mprop = 2353; %Propellant mass (kg)
mstruct = 2150; %Structural mass (kg)
T = 16000; %Total engine thrust (N)
Isp = 311; %Engine specific impulse (s)
mdot = T/(g*Isp); %Engine mass flow rate (kg/s)
tburn = mprop/mdot; %Engine burn time (s)
m0 = mprop + mstruct; %Initial launcher mass (kg)

x0 = 0; %Initial x-position (m)
y0 = R; %Initial y-position (m)
vx0 = 0; %Initial x-velocity (m/s)
vy0 = 0; %Initial y-velocity (m/s)

rorbit = R+100000; %Desired orbital altitude above Moon's surface (m)
vorbit = sqrt((G*M)/rorbit); %Orbital velocity required to maintain circular orbit at altitude r = rorbit

tmax = tburn%tcutoff; %Simulation running time (s)
tspan = [0:tmax]; %Use : notation to produce finer integration steps
xinit = [x0; y0; vx0; vy0; m0]; %Initial conditions

[t,x] = ode23t(@(t,x) equations_of_motion(t, x, G, g, M, Isp, T, tburn, parameters), tspan, xinit);
vtotal = sqrt(x(:,3).^2 + x(:,4).^2); %Orbital velocity
vradial = (x(:,1).*x(:,3) + x(:,2).*x(:,4))./sqrt(x(:,1).^2 + x(:,2).^2); %Radial velocity

%Plotting
plot(x(:,1), x(:,2))
xlabel('x-position');
ylabel('y-position');
pbaspect([1 1 1])
grid on
figure
plot(t, vtotal)
xlabel('time (s)');
ylabel('Orbital velocity (m/s)');
grid on
figure
plot(1:length(t),vradial)
xlabel('time (s)');
ylabel('Radial velocity (m/s)');
grid on

The second file is cost_function.m, which contains the call to ode45 and calculates the objective function to be minimized by fminunc:

function J = cost_function(parameters)
%function [outputs] = function(inputs)

%x = [x1, x2, x3, x4, x5] = state vector = output from ODE solver
%x1 = x-position
%x2 = y-position
%x3 = x-velocity
%x4 = y-velocity
%x5 = mass flow rate

%x(end, 1) = x-position at t = tfinal
%x(end, 2) = y-position at t = tfinal
%x(end, 3) = x-velocity at t = tfinal
%x(end, 4) = y-velocity at t = tfinal

G = 6.672e-11; %Gravitational constant
M = 7.3476e22; %Mass of Moon (kg)
R = 1.737e6; %Radius of Moon (m)
g = 9.81; %Gravitational acceleration on Earth (m/s^2);
mprop = 2353; %Propellant mass (kg)
mstruct = 2150; %Structural mass (kg)
T = 16000; %Total engine thrust (N)
Isp = 311; %Engine specific impulse (s)
mdot = T/(g*Isp); %Engine mass flow rate (kg/s)
tburn = mprop/mdot; %Engine burn time (s)
m0 = mprop + mstruct; %Initial launcher mass (kg)

x0 = 0; %Initial x-position (m)
y0 = R; %Initial y-position (m)
vx0 = 0; %Initial x-velocity (m/s)
vy0 = 0; %Initial y-velocity (m/s)

rorbit = R+100000; %Desired orbital altitude above Moon's surface (m)
vorbit = sqrt((G*M)/rorbit); %Orbital velocity required to maintain circular orbit at altitude r = rorbit

tmax = tburn;%tcutoff; %Simulation running time (s)
tspan = [0:tmax]; %Use : notation to produce finer integration steps
xinit = [x0; y0; vx0; vy0; m0]; %Initial conditions

[t,x] = ode45(@(t,x) equations_of_motion(t, x, G, g, M, Isp, T, tburn, parameters), tspan, xinit);
vtotal = sqrt(x(:,3).^2 + x(:,4).^2); %Orbital velocity
vradial = (x(:,1).*x(:,3) + x(:,2).*x(:,4))./sqrt(x(:,1).^2 + x(:,2).^2); %Radial velocity

%Objective/cost function to minimize
%1) Spacecraft must be at r = rorbit when t = tcutoff
%2) Spacecraft must be at v = vorbit when t = tcutoff
%3) Spacecraft must have vradial = 0 when t = tcutoff
J = (rorbit - sqrt(x(end,1).^2+x(end,2).^2)).^2 + (vorbit - sqrt(x(end,3).^2+x(end,4).^2)).^2 + (0 - vradial(end)).^2
end

And the third file, equations_of_motion.m, contains the system of ODEs and the parameter values (theta0, thetaf, tcutoff) to be optimized:

function xdot = equations_of_motion(t, x, G, g, M, Isp, T, tburn, parameters)
%function [ output_args ] = f( input_args )
%t = current time
%x = current state vector (x1, x2, x3, x4, x5)
%g = gravitational acceleration

      %System of equations
      %x1 = x-position, x'(t) = vx(t) = x3
      %x2 = y-position, y'(t) = vy(t) = x4
      %x3 = x-velocity, vx'(t) = ax(t)
      %x4 = y-velocity, vy'(t) = ay(t)
      %x5 = launcher mass = m(t)
      %Note that variable mass x(5)(t) is used as an acceleration term in
      %xdot(3)(t) and xdot(4)(t)

      theta0 = parameters(1);
      thetaf = parameters(2);
      tcutoff = parameters(3);

      %Linear-tangent steering law
      LTS = atan(tan(theta0) - (tan(theta0) - tan(thetaf)).*(t./tburn));

      if t <= tcutoff
      xdot = [
              x(3); %xdot1: x-velocity
              x(4); %xdot2: y-velocity
              -(G*M*x(1))/((x(1)^2+x(2)^2)^(3/2))+(T/x(5))*cos(LTS); %xdot3: x-acceleration
              -(G*M*x(2))/((x(1)^2+x(2)^2)^(3/2))+(T/x(5))*sin(LTS); %xdot4: y-acceleration
              -T/(g*Isp) %xdot5: engine mass flow rate
              ];
      else
          xdot = [
              x(3); %xdot1: x-velocity
              x(4); %xdot2: y-velocity
              -(G*M*x(1))/((x(1)^2+x(2)^2)^(3/2)); %xdot3: x-acceleration
              -(G*M*x(2))/((x(1)^2+x(2)^2)^(3/2)); %xdot4: y-acceleration
              0 %xdot5: engine mass flow rate
              ];
      end
  end

I had previously written the simulation in Mathematica, where the optimal parameter values were theta0 = 1.15055 radians, thetaf = -0.867385 radians and tcutoff = 430.086 seconds, which I've added to the top of the main.m file for reference. It seems that fminunc works beautifully at finding the correct values for theta0 and thetaf. However, the value for tcutoff does not change from the initial guess value that I give using parameters_init in main.m, and I'm sure this is due to me incorrectly implementing the tcutoff parameter using an if statement in the equations_of_motion.m file. As such, I was wondering how I could better implement the third parameter, tcutoff, so that it is able to be optimized to the expected optimal value of tcutoff = 430.086 seconds, instead of not being changed at all. Any help would be much appreciated! Below are the comparison graphs from Mathematica's output, and the current output from MATLAB:

Mathematica:

Cartesian xy position when engine turns off:

Orbital velocity when engine turns off:

Radial velocity when engine turns off:

MATLAB:

Cartesian xy position when engine turns off:

Orbital velocity when engine turns off:

Radial velocity when engine turns off:

Answer 1

[Alan Weiss]

There is too much here for me to be able to diagnose anything in short order. However, I wonder whether you might just need to take larger finite difference sizes and maybe central finite differences as suggested in the documentation.

Also, for minimizing a nonlinear sum of squares, it is sometimes more beneficial to use lsqnonlin instead of fminunc, making sure to modify the objective function to return the three components that you currently square and sum; I mean, for lsqnonlin, don't do the summing and squaring in the objective.

Alan Weiss

MATLAB mathematical toolbox documentation