Hang Glider Control

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

Contents

Problem Formulation

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

$$ J = x $$

subject to:

$$ \frac{d^{2}x}{dt^{2}} = \frac{1}{m}*(-L*sin(neta)-D*cos(neta)) $$

$$ \frac{d^{2}y}{dt^{2}} = \frac{1}{m}*(L*cos(neta)-D*sin(neta)) - g $$

$$ sin(neta) = \frac{w}{v} $$

$$ cos(neta) = \frac{\frac{dx}{dt}}{v} $$

$$ v = \sqrt{(\frac{dx}{dt})^2+w^2} $$

$$ w = \frac{dy}{dt}-u $$

$$ u = u_0*(1-r)*exp^{-r} $$

$$ r = (\frac{x}{r_0}-2.5)^2 $$

$$ u_0 = 2.5 $$

$$ r_0 = 100 $$

$$ D = \frac{1}{2}*(c_0+c_1**c_L^2)*rho*S*v^2 $$

$$ L = \frac{1}{2}*c_L*rho*S*v^2 $$

$$ c_0 = 0.034 $$

$$ c_1 = 0.069662 $$

$$ S = 14 $$

$$ rho = 1.13 $$

$$ 0 <= c_L <= c_{max} $$

$$ x >= 0 $$

$$ \frac{dx}{dt} >= 0 $$

$$ c_{max} = 1.4 $$

$$ m = 100 $$

$$ g = 9.81 $$

$$ [x_0 \ y_0] = [0 \ 1000] $$

$$ [y_{t_f}] = 900 $$

$$ [\frac{dx}{dt}_0 \ \frac{dy}{dt}_0] = [13.23 \ -1.288] $$

$$ [\frac{dx}{dt}_{t_f} \ \frac{dy}{dt}_{t_f}] = [13.23 \ -1.288] $$

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