Verbesserungen am Skript für die Finite Burn Trajectory Optimization mit MATLAB

Ich habe ein MATLAB-Skript geschrieben, um den Zweipunkt-Grenzgleichungslöser bvp4c zu verwenden, um eine Trajektorienoptimierung des Calculus of Variations des Problems einer einzelnen endlichen Verbrennung mit einer freien Küste zu implementieren, um den optimalen Zündpunkt auszuwählen. Ich hätte gerne Beiträge zu Fehlern, die ich gemacht habe, oder Möglichkeiten, den Code zu verbessern. Was es tut, ist die optimale endliche Verbrennung von einer Umlaufbahn von 185 x 185 km auf eine Umlaufbahn von 185 x 1000 km mit einer optionalen Neigungsänderung. Es stellt ein Mehrpunkt-Grenzwertproblem auf, um ein Burn-Coast-Problem unter Verwendung der Primer-Vektorgleichung zu lösen.

Der Code basiert auf MATLAB-Code in den Anhängen von Optimal Control with Aerospace Applications von Longuski, Guzman und Prussing. Die Auslaufzeit und Brennzeiten wurden als freie Parameter implementiert und das Grenzwertproblem verwendet eine normalisierte Zeit (Tau), sodass [0,1] die Auslaufzeit und [1,2] die Brenndauer ist.

Die ODE implementiert die Bewegungsgleichungen und die Grundvektorgleichung für die Zentralkraftbewegung in kartesischen 3D-Koordinaten. Ich glaube, ich habe diese richtig verstanden, könnte aber jemanden gebrauchen, der meine Arbeit noch einmal überprüft. Die Multiplikation mit den Küsten- oder Brennzeiten in den Rückgabewert der ODE ist notwendig, um die Ableitung umzuwandeln$\frac{d}{dt}$zu$\frac{d}{d \tau}$.

Die Randbedingungen sind etwas, zu dem ich eine Frage habe, ob ich es richtig gemacht habe. Anfangslage, Geschwindigkeit und Masse sind triviale Bedingungen. Es gibt auch 14 Bedingungen, die allen Variablen auf der linken und rechten Seite zwischen der Küste und dem Brand (Tau = 1) Kontinuität auferlegen.

Die verbleibenden Bedingungen sind die 5 Beschränkungen für die Endbedingungen der Umlaufbahn und dann die 3 Transversalitätsbedingungen, um die wahre Anomalie zu berücksichtigen, dass der Befestigungspunkt frei ist, und die Zeiten der Küste und der Verbrennung, die frei sind. Für die Endbahnbedingungen verwende ich den Drehimpulsvektor und die ersten beiden Komponenten des Exzentrizitätsvektors sind gleich der Zielbahn. Was ich für Transveralitätsbedingungen verwende, ist das, wenn der Hamiltonian zerlegt wird$H = H_0 + T S$wobei T die Schubkraft und S die Schaltfunktion ist, dann stelle ich ein$H_0(t_2) = 0$,$H_0(t_0) = 0$und dann mit$S = ( \lvert p_v \rvert - 1)$einstellen$\lvert p_v(t_1) \rvert = 1$. Mir fehlt etwas Selbstvertrauen, dass ich das ganz richtig verstanden habe.

Die Ergebnisse scheinen ziemlich anständig zu sein, aber ich bin ein wenig überrascht, dass ich es nicht sehe$\lvert p_v(t_2) \rvert = 1$. Ich sehe auch Konvergenzprobleme bei Problemen mit einer Neigung über etwa 25 Grad (aber das geht in 300 Sekunden / 3000 dV-Verbrennungen).

Ich habe die Einheiten für Position und Geschwindigkeit noch nicht normalisiert, also sind sie noch drin$m$und$m/s$.

Für mäßigere Verbrennungen scheint es ziemlich gut zu funktionieren und liefert Brennzeiten innerhalb von etwa einer Sekunde des impulsiven Verbrennungsmodells und sieht so aus, als ob es die Verbrennung mit einer negativen Küste von ungefähr der Hälfte der Brennzeit korrekt anführt.

Ich würde mich jedoch über jede Hilfe freuen, um Fehler zu finden, die ich gemacht habe, oder Möglichkeiten, diesen Code zu verbessern, um ihn robuster zu machen.

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