PROPT Hang Glider Control
From TomWiki
Jump to navigationJump to search
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');