# Hang Glider Control

Benchmarking Optimization Software with COPS Elizabeth D. Dolan and Jorge J. More ARGONNE NATIONAL LABORATORY

## Problem Formulation

Find u(t) over t in [0; tF ] to maximize

subject to:

cL is the control variable.

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


## Problem setup

toms t
toms tf

for n=[10 80]

    p = tomPhase('p', t, 0, tf, n, [], 'cheb');
setPhase(p);

tomStates x dx y dy
tomControls cL

% Initial guess
% Note: The guess for tf must appear in the list before
% expression involving t.
if n == 10
x0 = {tf == 105
icollocate({
dx == 13.23; x  == dx*t
dy == -1.288; y  == 1000+dy*t
})
collocate(cL==1.4)};
else
x0 = {tf == tf_opt
icollocate({
dx == dx_opt; x  == x_opt
dy == dy_opt; y  == y_opt
})
collocate(cL == cL_opt)};
end

% Box constraints
cbox = {
0.1 <= tf <= 200
0   <= icollocate(x)
0   <= icollocate(dx)
0   <= collocate(cL) <= 1.4};

% Boundary constraints
cbnd = {initial({x  == 0; dx == 13.23
y  == 1000; dy == -1.288})
final({dx == 13.23; y  == 900; dy == -1.288})};

% Various constants and expressions
m = 100;      g = 9.81;
uc = 2.5;     r0 = 100;
c0  = 0.034;  c1  = 0.069662;
S   = 14;     rho = 1.13;

r = (x/r0-2.5).^2;
u = uc*(1-r).*exp(-r);
w = dy-u;
v = sqrt(dx.^2+w.^2);
sinneta = w./v;
cosneta = dx./v;
D = 1/2*(c0+c1*cL.^2).*rho.*S.*v.^2;
L = 1/2*cL.*rho.*S.*v.^2;

% ODEs and path constraints
ceq = collocate({
dot(x)  == dx
dot(dx) == 1/m*(-L.*sinneta-D.*cosneta)
dot(y)  == dy
dot(dy) == 1/m*(L.*cosneta-D.*sinneta)-g
dx.^2+w.^2 >= 0.01});

% Objective
objective = -final(x);


## Solve the problem

    options = struct;
options.name = 'Hang Glider';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);

Problem type appears to be: lpcon
===== * * * =================================================================== * * *
TOMLAB - Tomlab Optimization Inc. Development license  999001. Valid to 2010-02-05
=====================================================================================
Problem: ---  1: Hang Glider                    f_k   -1281.388593956428200000
sum(|constr|)      0.000000000068965670
f(x_k) + sum(|constr|)  -1281.388593956359300000
f(x_0)  -1389.149999999998700000

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

FuncEv    1 ConstrEv   69 ConJacEv   69 Iter   48 MinorIter  225
CPU time: 0.265625 sec. Elapsed time: 0.266000 sec.

Problem type appears to be: lpcon
===== * * * =================================================================== * * *
TOMLAB - Tomlab Optimization Inc. Development license  999001. Valid to 2010-02-05
=====================================================================================
Problem: ---  1: Hang Glider                    f_k   -1305.253702077172000000
sum(|constr|)      0.000000043627467339
f(x_k) + sum(|constr|)  -1305.253702033544400000
f(x_0)  -1281.388593956420700000

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

FuncEv    1 ConstrEv   80 ConJacEv   80 Iter   67 MinorIter 1182
CPU time: 4.125000 sec. Elapsed time: 4.172000 sec.


## Extract optimal states and controls from solution

    x_opt  = subs(x,solution);
dx_opt = subs(dx,solution);
y_opt  = subs(y,solution);
dy_opt = subs(dy,solution);
cL_opt = subs(cL,solution);
tf_opt = subs(tf,solution);

end


## Plot result

figure(1)
ezplot(x,y);
xlabel('Hang Glider x');
ylabel('Hang Glider y');
title('Hang Glider trajectory.');

figure(2)
subplot(2,1,1)
ezplot([dx; dy]);
legend('vx','vy');
title('Hang Glider speeds dxdt and dydt');

subplot(2,1,2)
ezplot(cL);
legend('cL');
title('Hang Glider lift coefficient');