Code covered by the BSD License  

Highlights from
Aerospace Trajectory Optimization Using Direct Transcription and Collocation

Aerospace Trajectory Optimization Using Direct Transcription and Collocation

by

David Eagle (view profile)

 

06 Nov 2012 (Updated )

Demonstrates the solution of an aerospace trajectory optimization problem.

dtofunc_trap (x)
function [f, g] = dtofunc_trap (x)

% objective function and equality constraints

% trapezoidal collocation method

% inputs

%  x = current nlp variable values

% outputs

%  f = vector of equality constraints and
%      objective function evaluated at x

% Orbital Mechanics with MATLAB

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

global nc_defect ndiffeq nnodes nlp_state ncv tarray

% compute state vector defect equality constraints

for k = 1:1:nnodes - 1

    % time elements

    tk = tarray(k);

    tkp1 = tarray(k + 1);

    % state vector elements

    if (k == 1)
        
       % first node
       
       nks = 0;
       
       nkp1s = ndiffeq;
       
    else
        
       % reset to previous node
       
       nks = (k - 1) * ndiffeq;
       
       nkp1s = nks + ndiffeq;
       
    end
   
    for i = 1:1:ndiffeq
        
        xk(i) = x(nks + i);

        xkp1(i) = x(nkp1s + i);
        
    end

    % control variable elements

    if (k == 1)
        
       % first node
       
       nkc = nlp_state;
       
       nkp1c = nkc + ncv;
       
    else
        
       % reset to previous node
       
       nkc = nlp_state + (k - 1) * ncv;
       
       nkp1c = nkc + ncv;
       
    end

    for i = 1:1:ncv
        
        uk(i) = x(nkc + i);

        ukp1(i) = x(nkp1c + i);
        
    end

    % compute state vector defects for current node

    reswrk = defect_trap(tk, tkp1, xk, xkp1, uk, ukp1);

    % load defect array for this node

    for i = 1:1:ndiffeq
        
        resid(nks + i) = reswrk(i);
        
    end
end

% set active defect constraints
% (offset by 1)

for i = 1:1:nc_defect
    
    f(i + 1) = resid(i);
    
end

% current final state vector

xfinal(1) = x(nlp_state - 2);

xfinal(2) = x(nlp_state - 1);

xfinal(3) = x(nlp_state);

% objective function (maximize final radius)

f(1) = -xfinal(1);

% compute auxillary equality constraints (final boundary conditions)

f(nc_defect + 2) = xfinal(2);

f(nc_defect + 3) = xfinal(1) * xfinal(3)^2 - 1.0;

% transpose

f = f';

% no derivatives

g = [];

Contact us