Hang Glider Control

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

Problem Formulation
Problem setup
Solve the problem
Extract optimal states and controls from solution
Plot result

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