PROPT Hang Glider Control

From TomWiki
Jump to navigationJump to search

Notice.png

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

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

Problem Formulation

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


subject to:



cL is the control variable.

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

Problem setup

toms t
toms t_f

for n=[10 80]


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

    tomStates x dx y dy
    tomControls cL

    % Initial guess
    % Note: The guess for t_f must appear in the list before
    % expression involving t.
    if n == 10
        x0 = {t_f == 105
            icollocate({
            dx == 13.23; x  == dx*t
            dy == -1.288; y  == 1000+dy*t
            })
            collocate(cL==1.4)};
    else
        x0 = {t_f == 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 <= t_f <= 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});

    % Oscilation penalty, to avoid overfitting
    penalty = 0.03*integrate(dot(dy).^2+dot(dx).^2);

    % Objective
    objective = -final(x)+penalty;

Solve the problem

    options = struct;
    options.name = 'Hang Glider';
    solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
Problem type appears to be: con
Time for symbolic processing: 0.65513 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license  999007. Valid to 2011-12-31
=====================================================================================
Problem: ---  1: Hang Glider                    f_k   -1280.926971802052500000
                                       sum(|constr|)      0.000000003742813828
                              f(x_k) + sum(|constr|)  -1280.926971798309800000
                                              f(x_0)  -1389.149999999999600000

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

FuncEv   62 GradEv   60 ConstrEv   60 ConJacEv   60 Iter   43 MinorIter  197
CPU time: 0.093601 sec. Elapsed time: 0.088000 sec. 

Problem type appears to be: con
Time for symbolic processing: 0.65458 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license  999007. Valid to 2011-12-31
=====================================================================================
Problem: ---  1: Hang Glider                    f_k   -1247.682043489516000000
                                       sum(|constr|)      0.000000072989492234
                              f(x_k) + sum(|constr|)  -1247.682043416526500000
                                              f(x_0)  -1280.808409831168300000

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

FuncEv   86 GradEv   84 ConstrEv   84 ConJacEv   84 Iter   63 MinorIter 1062
CPU time: 1.513210 sec. Elapsed time: 1.520000 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(t_f,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');

HangGlider 01.png

HangGlider 02.png