Усовершенствования скрипта для оптимизации траектории конечного прожига с помощью MATLAB

Я написал сценарий MATLAB, чтобы использовать решатель двухточечного граничного уравнения bvp4c для реализации оптимизации траектории с помощью вариационного исчисления проблемы одиночного конечного горения со свободным выбегом для выбора оптимальной точки воспламенения. Я хотел бы получить информацию об ошибках, которые я сделал, или о способах улучшения кода. Что он делает, так это определяет оптимальное конечное горение от орбиты 185x185 км до орбиты 185x1000 км с необязательным изменением наклонения. Он устанавливает многоточечную краевую задачу для решения задачи ожога-побережья с использованием векторного уравнения праймера.

Код основан на коде MATLAB в приложениях к книге « Оптимальное управление с аэрокосмическими приложениями » Лонгуски, Гусмана и Пруссинга. Время выбега и время горения были реализованы как свободные параметры, а в краевой задаче используется нормализованное время (тау), так что [0,1] — это период выбега, а [1,2] — период горения.

ODE реализует уравнения движения и векторное уравнение праймера для движения центральной силы в трехмерных декартовых координатах. Думаю, я все правильно понял, но мог бы кто-нибудь перепроверить мою работу. Умножение на время выбега или время горения в возвращаемом значении ODE необходимо для преобразования производной от$\frac{d}{dt}$к$\frac{d}{d \tau}$.

Граничные условия - это то, о чем у меня есть вопрос, правильно ли я все сделал. Начальное положение, скорость и масса — тривиальные условия. Есть также 14 условий, которые обеспечивают непрерывность всех переменных слева и справа между побережьем и горением (tau = 1).

Остальные условия - это 5 ограничений на конечные орбитальные условия, а затем 3 условия трансверсальности, учитывающие истинную аномалию свободного места точки присоединения, а также свободное время выбега и ожога. Для конечных орбитальных условий я использую вектор углового момента, а первые две компоненты вектора эксцентриситета равны целевой орбите. Что я использую для условий трансверальности, так это то, что если гамильтониан разбит на$H = H_0 + T S$где Т — тяга, а S — функция переключения, то устанавливаю$H_0(t_2) = 0$,$H_0(t_0) = 0$а потом с$S = ( \lvert p_v \rvert - 1)$поставил$\lvert p_v(t_1) \rvert = 1$. Мне не хватает уверенности в том, что я все понял правильно.

Результаты кажутся довольно приличными, но я немного удивлен, что не вижу$\lvert p_v(t_2) \rvert = 1$. Я также вижу проблемы с конвергенцией для любых проблем, связанных с наклоном выше примерно 25 градусов (но это приводит к сжиганию 300 секунд / 3000 dV).

Я еще не нормализовал единицы измерения положения и скорости, поэтому они все еще$m$и$m/s$.

Для более умеренных прожигов он, кажется, работает довольно хорошо и дает время прожига примерно в секунду по сравнению с моделью импульсного прожига и выглядит так, как будто он правильно ведет прожиг с отрицательным выбегом примерно в половину времени горения.

Я был бы признателен за любую помощь в поиске ошибок, которые я сделал, или способов улучшить этот код, чтобы сделать его более надежным.

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