PROPT Orbit Raising Maximum Radius: Difference between revisions

From TomWiki
Jump to navigationJump to search
No edit summary
No edit summary
Line 6: Line 6:


<math> J = r(t_f) </math>
<math> J = r(t_f) </math>


subject to the dynamic constraints
subject to the dynamic constraints


<math> \frac{d_r}{dt} = u </math>
<math> \frac{d_r}{dt} = u </math>
<math> \frac{du}{dt} = \frac{v^2}{r}-\frac{mmu}{r^2}+T*\frac{w1}{m} </math>
<math> \frac{du}{dt} = \frac{v^2}{r}-\frac{mmu}{r^2}+T*\frac{w1}{m} </math>
<math> \frac{dv}{dt} = -u*\frac{v}{r}+T*\frac{w2}{m} </math>
<math> \frac{dv}{dt} = -u*\frac{v}{r}+T*\frac{w2}{m} </math>
<math> \frac{dm}{dt} = -\frac{T}{g0*ISP} </math>
<math> \frac{dm}{dt} = -\frac{T}{g0*ISP} </math>


the boundary conditions
the boundary conditions


<math> r(t_0) = 1 </math>
<math> r(t_0) = 1 </math>
<math> u(t_0) = 1 </math>
<math> u(t_0) = 1 </math>
<math> u(t_f) = 0 </math>
<math> u(t_f) = 0 </math>
<math> v(t_0) = (\frac{mmu}{r(t_0)})^{0.5} </math>
<math> v(t_0) = (\frac{mmu}{r(t_0)})^{0.5} </math>
<math> v(t_f) = (\frac{mmu}{r(t_f)})^{0.5} </math>
<math> v(t_f) = (\frac{mmu}{r(t_f)})^{0.5} </math>
<math> m(t_0) = 1 </math>
<math> m(t_0) = 1 </math>


and the path constraint
and the path constraint


<math> w_1^2+w_2^2 = 1 </math>
<math> w_1^2+w_2^2 = 1 </math>


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

Revision as of 08:10, 9 November 2011

Notice.png

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

Problem description

Maximize:


subject to the dynamic constraints


the boundary conditions


and the path constraint


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');

OrbitRaisingMaxRadius 01.png