## 66  Min Energy Orbit Transfer

### 66.1  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)

### 66.2  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);
```

### 66.3  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 Optimization Inc. Development license  999001. Valid to 2011-02-05
=====================================================================================
Problem: ---  1: Min Energy Transfer            f_k       0.523205792272900430
sum(|constr|)      0.000002310623001201
f(x_k) + sum(|constr|)      0.523208102895901580
f(x_0)      0.569999999999999170

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

FuncEv  235 GradEv 1975 ConstrEv  234 ConJacEv 1975 Iter  206 MinorIter  231
CPU time: 9.203125 sec. Elapsed time: 9.203000 sec.
```
```===== * * * =================================================================== * * *
TOMLAB - Tomlab Optimization Inc. Development license  999001. Valid to 2011-02-05
=====================================================================================
Problem: ---  1: Min Energy Transfer            f_k       0.523085806417755370
sum(|constr|)      0.000000004935871044
f(x_k) + sum(|constr|)      0.523085811353626420
f(x_0)      0.522037502371971220

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

FuncEv  113 GradEv 1732 ConstrEv  112 ConJacEv 1732 Iter   79 MinorIter  109
CPU time: 19.593750 sec. Elapsed time: 20.375000 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');

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
```