Code covered by the BSD License  

Highlights from
Two-dimensional, Low-thrust Earth-to-Mars Trajectory Analysis with MATLAB

Two-dimensional, Low-thrust Earth-to-Mars Trajectory Analysis with MATLAB

by

 

27 Jun 2013 (Updated )

Determines optimal, two-dimensional low-thrust Earth-to-Mars interplanetary trajectories.

ilt_opt.m
% ilt_opt.m           May 6, 2014

% two-dimensional low-thrust Earth-to-Mars trajectory optimization

% minimize transfer time or maximize final mass
% with piecewise-linear steering

% time and state variables

%  t    = simulation time
%  y(1) = radial distance (r)
%  y(2) = radial component of velocity (u)
%  y(3) = tangential component of velocity (v)
%  y(4) = polar angle (theta)

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

global xmu aunit rkcoef tfactor nsegments thrmag mass0 iopt tof

global acc_thrust mdot mp_dot alpha_wrk xinitial xfinal

% conversion factors - radians to/from degrees

rtd = 180.0 / pi;

dtr = pi / 180.0;

% astronomical unit (kilometers)

aunit = 149597870.691;

% gravitational constant of the sun (au^3/day^2)

smu = 0.2959122082855912D-03;

% gravitational constant of the sun (kilometer^3/second^2)

xmu = smu * aunit^3 / 86400.0^2;

% non-dimensional to dimensional time conversion factor

tfactor = sqrt(1.0 / smu);

% initialize rkf78 algorithm

rkcoef = 1;

clc; home;

fprintf('\n\nprogram ilt_opt - Earth-to-Mars low-thrust trajectory\n');

% request name of simulation definition data file

[filename, pathname] = uigetfile('*.dat', 'Please select the input file to read');

% read data file

[fid, iopt, nsegments, mass0, thrmag, xisp, time_g, time_lb, ...
    time_ub, throttle_g, throttle_lb, throttle_ub, alpha_g, ...
    alpha_lb, alpha_ub, jd_guess] = read_ilt(filename);

% compute non-dimensional thrust acceleration ratio

acc_thrust = (0.001 * thrmag / mass0) / (xmu / aunit^2);

% compute normalized propellant flow rate

mdot = (thrmag / (9.80665 * xisp)) * 86400.0;

mp_dot = (mdot / mass0) * tfactor;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initial and final times and states
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% define state vector at initial time

xinitial(1) = 1.0;

xinitial(2) = 0.0;

xinitial(3) = 1.0;

xinitial(4) = 0.0;

% final conditions

xfinal(1) = 1.52368;

xfinal(2) = 0.0;

xfinal(3) = sqrt(1.0 / xfinal(1));

if (iopt == 1)
    
    % initial guess for non-dimensional transfer time
    
    xg(1) = time_g / tfactor;
    
else
    
    % initial guess for throttle setting
    
    xg(1) = throttle_g;
    
    tof = time_g / tfactor;
    
end

% initial guess for steering angles (radians)

xg(2:nsegments + 1) = alpha_g;

% transpose initial guess

xg = xg';

if (iopt == 1)
    
    % upper and lower bounds for non-dimensional transfer time
    
    xlb(1) = time_lb / tfactor;
    
    xub(1) = time_ub / tfactor;
    
else
    
    % upper and lower bounds for throttle setting
    
    xlb(1) = throttle_lb;
    
    xub(1) = throttle_ub;
    
end

% upper and lower bounds for steering angles (radians)

xlb(2:nsegments + 1) = alpha_lb;

xub(2:nsegments + 1) = alpha_ub;

% transpose bounds

xlb = xlb';

xub = xub';

% define lower and upper bounds on objective function (transfer time)

flow(1) = 0.0;

fupp(1) = +Inf;

% define bounds on final state vector equality constraints

flow(2) = xfinal(1);
fupp(2) = xfinal(1);

flow(3) = xfinal(2);
fupp(3) = xfinal(2);

flow(4) = xfinal(3);
fupp(4) = xfinal(3);

flow = flow';

fupp = fupp';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve low-thrust shooting problem
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

snscreen('on');

[x, f, inform, xmul, fmul] = snopt(xg, xlb, xub, flow, fupp, 'ilt_shoot');

% print results

fprintf('\n\nprogram ilt_opt - Earth-to-Mars low-thrust trajectory\n');

if (iopt == 1)
    
    fprintf('\nminimize transfer time\n');
    
else
    
    fprintf('\nmaximize final mass\n');
    
end

fprintf ('\nnumber of segments  %2i \n', nsegments);

fprintf('\ninitial mass            %12.6f kilograms\n', mass0);

fprintf('\nthrust magnitude        %12.6f newtons\n', thrmag);

fprintf('\nspecific impulse        %12.6f seconds\n', xisp);

fprintf('\ninitial state vector \n');

fprintf('\n  radius                %12.6f (AU)\n', xinitial(1));
fprintf('\n  radial velocity       %12.6f (AU/day)\n', xinitial(2));
fprintf('\n  transverse velocity   %12.6f (AU/day)\n', xinitial(3));

fprintf('\n\nfinal state vector \n');

fprintf('\n  radius                %12.6f (AU)\n', f(2));
fprintf('\n  radial velocity       %12.6f (AU/day)\n', f(3));
fprintf('\n  transverse velocity   %12.6f (AU/day)\n\n', f(4));

fprintf('\ntransfer time           %12.6f (non-dimensional)\n', tof);
fprintf('\ntransfer time           %12.6f days \n', tfactor * tof);

if (iopt == 1)
    
    % minimize transfer time (fixed throttle setting)
    
    fprintf('\nthrottle setting        %12.6f \n', throttle_g);
    
else
    
    % maximize final mass (optimal throttle setting)
    
    fprintf('\nthrottle setting        %12.6f \n', x(1));
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create control, state and trajectory graphics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute duration of each time interval

if (iopt == 1)
    
    % minimize transfer time option
    
    deltat = x(1) / nsegments;
    
else
    
    % maximize final mass option
    
    deltat = tof / nsegments;
    
end

% specify number of differential equations

neq = 6;

% truncation error tolerance

tetol = 1.0e-10;

% initialize initial time

ti = -deltat;

% set initial conditions

yi(1) = xinitial(1);

yi(2) = xinitial(2);

yi(3) = xinitial(3);

yi(4) = xinitial(4);

yi(5) = 0.0;

yi(6) = 0.0;

% load initial graphics data

npts = 1;

x1(npts) = 0.0;

x2(npts) = yi(1);

x3(npts) = yi(2);

x4(npts) = yi(3);

x5(npts) = x(2);

x6(npts) = yi(4);

% step size guess (non-dimensional time)

h = (1200.0 / 86400.0) / tfactor;

% integrate for all segments

for i = 1:1:nsegments
    
    % current steering angle
    
    alpha_wrk = x(i + 1);
    
    % increment initial and final times
    
    ti = ti + deltat;
    
    tf = ti + deltat;
    
    % integrate from current ti to tf
    
    yfinal = rkf78('ilt_eqm', neq, ti, tf, h, tetol, yi);
    
    % create graphics data
    
    npts = npts + 1;
    
    x1(npts) = tf;
    
    x2(npts) = yfinal(1);
    
    x3(npts) = yfinal(2);
    
    x4(npts) = yfinal(3);
    
    x5(npts) = alpha_wrk;
    
    x6(npts) = yfinal(4);
    
    % reset integration vector
    
    yi = yfinal;
    
    % check for end of simulation
    
    if (tf >= tof)
        
        break;
        
    end
    
end

fprintf('\npropellant mass         %12.6f kilograms\n', mass0 * yfinal(5));

fprintf('\nfinal mass              %12.6f kilograms\n\n', ...
    mass0 * (1.0 - yfinal(5)));

% determine departure and arrival calendar dates

transfer_time = tfactor * tof;

transfer_angle = rtd * x6(end);

[jd_depart, tsynodic] = departure(jd_guess, transfer_time, transfer_angle);

[cd_depart, ut_depart] = jd2str(jd_depart);

fprintf('\ndeparture calendar date      ');

disp(cd_depart);

jd_arrive = jd_depart + transfer_time;

[cd_arrive, ut_arrive] = jd2str(jd_arrive);

fprintf('\narrival calendar date        ');

disp(cd_arrive);

fprintf('\ntransfer angle          %12.6f degrees\n', transfer_angle);

fprintf('\nsynodic period          %12.6f days', tsynodic);
fprintf('\n                        %12.6f years\n\n', tsynodic / 365.25);

%%%%%%%%%%%%%%%%%%%%%%%%%
% create graphic displays
%%%%%%%%%%%%%%%%%%%%%%%%%

t = (jd_depart - 2451545.0) / 36525.0;

t2 = t * t;

% celestial longitude of the Earth at departure (radians)

theta_earth =  mod(dtr * (100.466449 + 35999.3728519 * t - 0.00000568 * t2), 2.0 * pi);

figure(1);

x(end+1) = x(end);

stairs(tfactor * x1(1:end), rtd * x(2:end));

title('Two-dimensional Earth-to-Mars Low-thrust Trajectory', 'FontSize', 16);

xlabel('mission elapsed time (days)', 'FontSize', 12);

ylabel('steering angle (degrees)', 'FontSize', 12);

grid;

print -depsc -tiff -r300 control.eps

figure(2);

plot(tfactor * x1, x2);

title('Two-dimensional Earth-to-Mars Low-thrust Trajectory', 'FontSize', 16);

xlabel('mission elapsed time (days)', 'FontSize', 12);

ylabel('radial distance (au)', 'FontSize', 12);

grid;

print -depsc -tiff -r300 rdistance.eps

figure(3);

plot(tfactor * x1, x3);

title('Two-dimensional Earth-to-Mars Low-thrust Trajectory', 'FontSize', 16);

xlabel('mission elapsed time (days)', 'FontSize', 12);

ylabel('radial velocity (au/day)', 'FontSize', 12);

grid;

print -depsc -tiff -r300 vradial.eps

figure(4);

plot(tfactor * x1, x4);

title('Two-dimensional Earth-to-Mars Low-thrust Trajectory', 'FontSize', 16);

xlabel('mission elapsed time (days)', 'FontSize', 12);

ylabel('transverse velocity (au/day)', 'FontSize', 12);

grid;

print -depsc -tiff -r300 vtransverse.eps

figure(5);

plot(tfactor * x1, rtd * (x6 + theta_earth));

title('Two-dimensional Earth-to-Mars Low-thrust Trajectory', 'FontSize', 16);

xlabel('mission elapsed time (days)', 'FontSize', 12);

ylabel('polar angle (degrees)', 'FontSize', 12);

grid;

print -depsc -tiff -r300 polar_angle.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create trajectory display
%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(6);

hold on;

axis([-1.7 1.7 -1.6 1.6]);

plot (0, 0, 'y*');

% create array of polar angles (radians)

t = 0: pi / 50.0: 2.0 * pi;

% plot Earth orbit

plot(xinitial(1) * sin(t), xinitial(1) * cos(t), 'Color', 'b');

% plot destination planet orbit

plot(f(2) * sin(t), f(2) * cos(t), 'Color', 'r');

% plot and label transfer orbit

plot(x2.* cos(x6 + theta_earth), x2.* sin(x6 + theta_earth), 'Color', 'k');

plot(x2(1) * cos(x6(1) + theta_earth), x2(1) * sin(x6(1) + theta_earth), 'ko');

plot(x2(end) * cos(x6(end) + theta_earth), x2(end) * sin(x6(end) + theta_earth), 'ko');

% label plot and axes

xlabel('X coordinate (AU)', 'FontSize', 12);

ylabel('Y coordinate (AU)', 'FontSize', 12);

title('Two-dimensional Earth-to-Mars Low-thrust Trajectory', 'FontSize', 16);

grid;

axis equal;

zoom on;

print -depsc -tiff -r300 trajectories.eps

Contact us