PROPT Orbit Raising Maximum Radius

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

Problem description

Maximize:

$J=r(t_{f})$

subject to the dynamic constraints

${\frac {d_{r}}{dt}}=u$ ${\frac {du}{dt}}={\frac {v^{2}}{r}}-{\frac {mmu}{r^{2}}}+T*{\frac {w1}{m}}$ ${\frac {dv}{dt}}=-u*{\frac {v}{r}}+T*{\frac {w2}{m}}$ ${\frac {dm}{dt}}=-{\frac {T}{g0*ISP}}$

the boundary conditions

$r(t_{0})=1$ $u(t_{0})=1$ $u(t_{f})=0$ $v(t_{0})=({\frac {mmu}{r(t_{0})}})^{0.5}$ $v(t_{f})=({\frac {mmu}{r(t_{f})}})^{0.5}$ $m(t_{0})=1$

and the path constraint

$w_{1}^{2}+w_{2}^{2}=1$

The control pitch angel is not being used in this formulation. Instead two control variables (w1,w2) are used to for the thrust direction. A path constraint ensures that (w1,w2) is a unit vector.

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

Problem setup

t_F    = 3.32;
mmu    = 1;
m_0   = 1;   r_0     = 1;            u_0         = 0;
u_f   = 0;   v_0     = sqrt(mmu/r_0); rmin        = 0.9;
rmax  = 5;   umin    = -5;           umax        = 5;
vmin  = -5;  vmax    = 5;            mmax        = m_0;
mmin  = 0.1;

T  = 0.1405;
Ve = 1.8758;

toms t
p1 = tomPhase('p1', t, 0, t_F, 40);
setPhase(p1);

tomStates r u v m
tomControls w1 w2

% Initial guess
x0 = {icollocate({
    r == r_0
    u == u_0 + (u_f-u_0)*t/t_F
    v == v_0
    m == m_0})
    collocate({w1 == 0; w2 == 1})};

% Boundary constraints
cbnd = {initial({
    r == r_0
    u == u_0
    v == v_0
    m == m_0
    })
    final({u == u_f
    v - sqrt(mmu/r) == 0})};

% Box constraints
cbox = {
    rmin <= icollocate(r) <= rmax
    umin <= icollocate(u) <= umax
    vmin <= icollocate(v) <= vmax
    mmin <= icollocate(m) <= mmax
    -1 <= collocate(w1) <= 1
    -1 <= collocate(w2) <= 1};

% ODEs and path constraints
ceq = collocate({
    dot(r) == u
    dot(u) == v^2/r-mmu/r^2+T*w1/m
    dot(v) == -u*v/r+T*w2/m
    dot(m) == -T/Ve
    w1^2+w2^2 == 1});

% Objective
objective = -final(r);

Solve the problem

options = struct;
options.name = 'Orbit Raising Problem Max Radius';
options.solver = 'snopt';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);

Problem type appears to be: lpcon
Time for symbolic processing: 0.22212 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license  999007. Valid to 2011-12-31
=====================================================================================
Problem: ---  1: Orbit Raising Problem Max Radius  f_k      -1.518744202740535800
                                          sum(|constr|)      0.000017265830060524
                                 f(x_k) + sum(|constr|)     -1.518726936910475200
                                                 f(x_0)     -0.999999999999991120

Solver: snopt.  EXIT=0.  INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv    1 ConstrEv   59 ConJacEv   59 Iter   36 MinorIter  346
CPU time: 0.202801 sec. Elapsed time: 0.226000 sec. 

Plot result

subplot(2,1,1)
ezplot([r u v m]);
legend('r','u','v','m');
title('Orbit Raising Problem Max Radius state variables');

subplot(2,1,2)
ezplot([w1 w2 atan2(w1,w2)])
legend('w_1', 'w_2', '\phi');
title('Orbit Raising Problem Max Radius control');