PROPT Min Energy Orbit Transfer

This page is part of the PROPT Manual. See PROPT Manual.

Problem description

Minimum energy orbit transfer of a spacecraft with limited variable thrust.

From the paper: Low thrust, high-accuracy trajectory optimization, I. Michael Ross, Qi Gong and Pooya Sekhavat.

Programmers: Li Jian (Beihang University)

% Copyright (c) 2009-2009 by Tomlab Optimization Inc.

Problem setup

% Array with consecutive number of collocation points
narr = [40 80];

toms t;
toms t_f;
t0 = 0;
tfmax = 57;

for n=narr

    %p = tomPhase('p', t, 0, t_f-t0, n, [], 'cheb');
    p = tomPhase('p', t, 0, t_f-t0, n);
    setPhase(p)

    tomStates r theta vr vt
    tomControls ur ut

    % Parameters
    r0 = 1; theta0 = 0; vr0 = 0; vt0 = 1;
    r_f = 4;             vrf = 0; vtf = 0.5;
    umax = 0.01;

    % Initial state
    xi=[r0; theta0; vr0; vt0];

    % Initial guess
    if n==narr(1)
        x0 = {t_f == 57;
            icollocate({r == xi(1); theta == xi(2); vr == xi(3); vt == xi(4)})
            collocate({ur == 0; ut == umax})};
    else
        x0 = {t_f == tfopt;
            icollocate({r == xopt1; theta == xopt2; vr == xopt3; vt == xopt4});
            collocate({ur == uopt1; ut == uopt2})};
    end

    % Box constraints
    cbox = {10 <= t_f<= tfmax;
        1 <= collocate(r) <= 4;
        0 <= collocate(vr) <= 0.5;
        0 <= collocate(vt) <= 1;
        -umax <= collocate(ur) <= umax;
        -umax <= collocate(ut) <= umax};

    % Boundary constraints
    cbnd = {initial({r == r0; theta == theta0; vr == vr0; vt == vt0})
        final({r == r_f; vr == 0; vt == vtf})};

    % ODEs and path constraints
    d_r      = vr;
    dtheta  = vt./r;
    dvr     = vt.*vt./r - 1.0./r./r + ur;
    dvt     = -vr.*vt./r + ut;

    ceq = collocate({
        dot(r) == d_r;
        dot(theta) == dtheta;
        dot(vr) == dvr;
        dot(vt) == dvt;
        0<=(ur.*ur+ut.*ut).^0.5<=umax});

    % Objective
    objective = integrate((ur.^2+ut.^2).^0.5);

Solve the problem

    options = struct;
    options.type = 'con';
    options.name = 'Min Energy Transfer';
    Prob = sym2prob(objective, {cbox, cbnd, ceq}, x0, options);
    Prob.KNITRO.options.ALG = 1;
    Prob.KNITRO.options.HONORBNDS = 0;
    Result = tomRun('knitro', Prob, 1);
    solution = getSolution(Result);

    xopt1 = subs(r,solution);
    xopt2 = subs(theta,solution);
    xopt3 = subs(vr,solution);
    xopt4 = subs(vt,solution);
    uopt1 = subs(ur,solution);
    uopt2 = subs(ut,solution);
    tfopt = subs(t_f,solution);

===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license  999007. Valid to 2011-12-31
=====================================================================================
Problem: ---  1: Min Energy Transfer            f_k       0.523146271130297880
                                       sum(|constr|)      0.000000048726915650
                              f(x_k) + sum(|constr|)      0.523146319857213490
                                              f(x_0)      0.569999999999998950

Solver: KNITRO.  EXIT=0.  INFORM=0.
Interior/Direct NLP KNITRO
Locally optimal solution found

FuncEv  873 GradEv 6168 ConstrEv  872 ConJacEv 6168 Iter  605 MinorIter  852
CPU time: 8.283653 sec. Elapsed time: 8.325000 sec. 

===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license  999007. Valid to 2011-12-31
=====================================================================================
Problem: ---  1: Min Energy Transfer            f_k       0.523087180930804970
                                       sum(|constr|)      0.000000000489698152
                              f(x_k) + sum(|constr|)      0.523087181420503140
                                              f(x_0)      0.521985852703766070

Solver: KNITRO.  EXIT=0.  INFORM=0.
Interior/Direct NLP KNITRO
Locally optimal solution found

FuncEv  453 GradEv 3679 ConstrEv  452 ConJacEv 3679 Iter  262 MinorIter  434
CPU time: 15.568900 sec. Elapsed time: 15.573000 sec. 

end

% Get final solution
t = subs(icollocate(t),solution);
r = subs(icollocate(r),solution);
theta = subs(icollocate(theta),solution);
vr = subs(icollocate(vr),solution);
vt = subs(icollocate(vt),solution);
ur = subs(icollocate(ur),solution);
ut = subs(icollocate(ut),solution);

t1 = 0:0.5:solution.t_f;
r_inter  = interp1p(t,r,t1);
theta_inter  = interp1p(t,theta,t1);
ur_inter = interp1p(t,ur,t1);
ut_inter = interp1p(t,ut,t1);

% Plot final solution
figure(1)
subplot(2,2,1)
plot(t,r,'*-');
legend('r');
title('radius');

subplot(2,2,2)
plot(t,theta,'*-');
legend('theta');
title('flight angle');

subplot(2,2,3)
plot(t,vr,'r*-',t,vt,'b*-');
legend('vr','vt');
title('velocity');

subplot(2,2,4)
plot(t,ur,'r+-',t,ut,'b+-');
legend('ur','ut');
title('Spacecraft controls');

figure(2)
polar(theta_inter,r_inter,'*-')
grid on
axis equal