Mejoras en el script para la optimización de trayectoria de quemado finito con MATLAB

Escribí un script de MATLAB para usar el solucionador de ecuaciones de límite de dos puntos bvp4c para implementar una optimización de trayectoria de Cálculo de variaciones del problema de una sola quemadura finita con una costa libre para elegir el punto de ignición óptimo. Me gustaría cualquier aporte sobre los errores que he cometido o formas de mejorar el código. Lo que hace es resolver el quemado finito óptimo desde una órbita de 185x185 km hasta una órbita de 185x1000 km, con un cambio de inclinación opcional. Establece un problema de valor límite multipunto para resolver un problema de costa quemada utilizando la ecuación del vector de imprimación.

El código se basa en el código de MATLAB en los apéndices de Optimal Control with Aerospace Applications de Longuski, Guzman y Prussing. El tiempo de inercia y los tiempos de quema se han implementado como parámetros libres y el problema del valor límite utiliza el tiempo normalizado (tau) de modo que [0,1] es el período de inercia y [1,2] es el período de quema.

La ODE implementa las ecuaciones de movimiento y la ecuación del vector primario para el movimiento de la fuerza central en coordenadas cartesianas 3D. Creo que los he entendido bien, pero me vendría bien que alguien revisara mi trabajo dos veces. La multiplicación por la costa o tiempos de quema en el valor de retorno de la ODE es necesaria para convertir la derivada de$\frac{d}{dt}$a$\frac{d}{d \tau}$.

Las condiciones de contorno son algo sobre lo que tengo una duda sobre si lo he hecho correctamente. La posición inicial, la velocidad y la masa son condiciones triviales. También hay 14 condiciones que imponen continuidad en todas las variables del lado izquierdo y del lado derecho entre la costa y la quema (tau = 1).

Las condiciones restantes son las 5 restricciones en las condiciones orbitales terminales, y luego las 3 condiciones de transversalidad para dar cuenta de la verdadera anomalía de que el punto de conexión esté libre y los tiempos de la costa y la quema estén libres. Para las condiciones orbitales terminales, estoy usando el vector de momento angular y los dos primeros componentes del vector de excentricidad son iguales a la órbita objetivo. Lo que estoy usando para las condiciones de transversalidad es que si el hamiltoniano se divide en$H = H_0 + T S$donde T es el empuje y S la función de conmutación, entonces configuro$H_0(t_2) = 0$,$H_0(t_0) = 0$y luego con$S = ( \lvert p_v \rvert - 1)$colocar$\lvert p_v(t_1) \rvert = 1$. Me falta algo de confianza en que lo haya entendido del todo bien.

Los resultados parecen ser bastante decentes, pero estoy un poco sorprendido de que no veo$\lvert p_v(t_2) \rvert = 1$. También veo problemas de convergencia para cualquier problema que involucre una inclinación superior a unos 25 grados (pero eso se convierte en quemaduras de 300 segundos / 3000 dV).

Todavía no he normalizado las unidades de posición y velocidad, por lo que todavía están en$m$y$m/s$.

Para quemados más moderados, parece funcionar bastante bien y brinda tiempos de quemado dentro de aproximadamente un segundo del modelo de quemado impulsivo y parece que lidera el quemado correctamente con una costa negativa de aproximadamente la mitad del tiempo de quemado.

Sin embargo, agradecería cualquier ayuda para encontrar errores que he cometido o formas de mejorar este código para hacerlo más sólido.

close all; clear all; clc;

global r0 v0 m0 rT vT g0
global MU Thrust Mdot

MU = 3.9860044189e+14;  % earth
Thrust = 232.7 * 1000;  % N; LR-91 232.7 kN
m0 = 32.74 * 1000;      % kg; 32.74t - 1g start accel
isp = 316;              % sec; LR-91
g0 = 9.80665;           % m/s; standard gravity
ve = isp * g0;          % m/s
a0 = Thrust / m0        % m/s^2; initial accel
tau = ve / a0           % s; time to burn rocket to zero
Mdot = Thrust / ve;     % kg/sec

rearth = 6.371e+6;      % m; earth radius
r185 = rearth + 0.185e+6;
r1000 = rearth + 1.000e+6;

% initial position (185x185)
r0 = [ r185, 0, 0 ];
v185 = sqrt(MU/r185);
v0 = [ 0, v185, 0 ];

% target orbital parameters (185x1000)
rT = [ r185, 0, 0 ];
smaT = ( r185 + r1000 ) / 2;
inc = 10;  % degrees
vTm = sqrt(MU * (2/r185 - 1/smaT) );
vT = [ 0, vTm * cosd(inc), vTm * sind(inc)];

vBurn = vT - v0;
dV = norm(vBurn)

% list initial conds
yinit = [r0 v0 vBurn/norm(vBurn) 0 0 0 m0 ];  % initial state and costate
tb_guess = tau * (1 - exp(-dV/ve) )
tc_guess = - tb_guess / 2

% sanity checks on the impulsive burn model
mf_impulsive = m0 - tb_guess * Mdot
af_impulsive = Thrust / mf_impulsive / g0

% parameters
parameters = [ tc_guess, tb_guess ];

% number of timeslices
Nt = 41;
tau = [
  linspace(0,1,Nt)'
  linspace(1,2,Nt)'
];

% initial guess
solinit = bvpinit(tau, yinit, parameters);

% bump up the mesh by a bit
option=bvpset('Nmax',50000);

% solve
sol = bvp4c(@Burn_odes, @Burn_bcs, solinit, option);

% extract times
tc = sol.parameters(1)
tb = sol.parameters(2)

% pull out values for times
Z = deval(sol, tau);

% convert taus to times
for i=1:length(tau)
  if tau(i) <= 1
    time(i) = tau(i) * tc;
  else
    time(i) = tc + ( tau(i) - 1 ) * tb;
  end
end

% extract solution vars
x_sol = Z(1,:);
y_sol = Z(2,:);
z_sol = Z(3,:);
r_sol = sqrt( x_sol.^2 + y_sol.^2 + z_sol.^2 );
vx_sol = Z(4,:);
vy_sol = Z(5,:);
vz_sol = Z(6,:);
v_sol = sqrt( vx_sol.^2 + vy_sol.^2 + vz_sol.^2 );
pvx_sol = Z(7,:);
pvy_sol = Z(8,:);
pvz_sol = Z(9,:);
pv_sol = sqrt( pvx_sol.^2 + pvy_sol.^2 + pvz_sol.^2 );
m_sol = Z(13,:);

for i = 1:length(tau)
  r = Z(1:3, i);
  v = Z(4:6, i);
  h(i) = norm(cross(r,v));
end

% plots

figure;
subplot(3,2,1);
plot(time,m_sol/1000);
ylabel('mass, tons', 'fontsize', 14);
xlabel('Time, sec', 'fontsize', 14);
hold on;
grid on;

subplot(3,2,2);
plot(time,r_sol/1000);
ylabel('radius, km', 'fontsize', 14);
xlabel('Time, sec', 'fontsize', 14);
hold on;
grid on;

subplot(3,2,3);
plot(time,v_sol/1000);
ylabel('velocity, km/s', 'fontsize', 14);
xlabel('Time, sec', 'fontsize', 14);
hold on;
grid on;

subplot(3,2,4);
plot(time,h);
ylabel('h', 'fontsize', 14);
xlabel('Time, sec', 'fontsize', 14);
hold on;
grid on;

subplot(3,2,5);
plot(time,pv_sol);
ylabel('pv magnitude', 'fontsize', 14);
xlabel('Time, sec', 'fontsize', 14);
hold on;
grid on;

figure;
plot(x_sol/1000, y_sol/1000);
grid on;
axis equal
xlabel('x, km', 'fontsize', 14);
ylabel('y, km', 'fontsize', 14);


%
% Boundary conditions
%

function PSI = Burn_bcs(yleft, yright, parameters)

Y0 = yleft(:,1);
Y1 = yleft(:,2);  % == yright(:,1)
Y2 = yright(:,2);

global r0 v0 m0 rT vT MU g0 Thrust

ri = Y0(1:3);
vi = Y0(4:6);
pvi = Y0(7:9);
pri = Y0(10:12);
mi = Y0(13);
ri3 = norm(ri)^3;

rf = Y2(1:3);
vf = Y2(4:6);
pvf = Y2(7:9);
prf = Y2(10:12);
mf = Y2(13);
rf3 = norm(rf)^3;

r1 = Y1(1:3);
v1 = Y1(4:6);
pv1 = Y1(7:9);
pr1 = Y1(10:12);
r13 = norm(r1)^3;

hT = cross(rT, vT);
hf = cross(rf, vf);

eT = - (rT/norm(rT) + cross(hT,vT)/MU);
ef = - (rf/norm(rf) + cross(hf,vf)/MU);

emiss = ef - eT';

uf = pvf / norm(pvf);

H0ti = dot(pri, vi) - MU * dot(pvi, ri) / ri3;
H0t1 = dot(pr1, v1) - MU * dot(pv1, r1) / r13;
H0tf = dot(prf, vf) - MU * dot(pvf, rf) / rf3;
Htf = H0tf - Thrust * ( norm(pvf) - 1 );

PSI = [
    Y0(1:3) - r0'                       % 3 - initial position
    Y0(4:6) - v0'                       % 3 - initial velocity
    Y0(13) - m0                         % 1 - initial mass
    norm(Y1(7:9)) - 1                   % 1 - norm(pv1) = 1 (so H(t1) = 0 on both sides of t1)
    hf - hT'                            % 3 - terminal constraint on angular momentum
    emiss(1:2)                          % 2 - terminal constraint on ecc vector
    H0ti                                % 1 - H0(t0) = 0
    H0tf                                % 1 - H0(tf) = 0
    yleft(:,2) - yright(:,1)            % 14 - values at t1 are continuous
    ];

end

%
% Equations of motion
%

function dX_dtau = Burn_odes(tau, X, k, parameters)

global MU Thrust Mdot

if k == 1
  thr = 0;
  md = 0;
  tinterval = parameters(1);  % tc
else
  thr = Thrust;
  md = Mdot;
  tinterval = parameters(2);  % tb
end

r = X(1:3);
v = X(4:6);
pv = X(7:9);
pr = X(10:12);
m = X(13);

u = pv / norm(pv);

Fm = thr / m;
r3 = norm(r)^3;
r2 = dot(r,r);
r5 = norm(r)^5;

rdot  = v;
vdot  = - MU * r / r3 + Fm * u;
pvdot = pr;
prdot = - MU / r5 * ( 3 * r' * r - r2 * eye(3,3) ) * pv;

dX_dtau = tinterval*[rdot', vdot', pvdot', prdot', -md ];

end