I am trying to write a code to optimize the wing planform for a high altitude UAV. I have 18 design variables that determine the shape of a 2 panel wing. These are listed in the attached code, with (i) designating panel 1 and (ii) designating panel 2. Airfoil data is supplied as a list of constants.
These feed into a lifting line code (from http://www.dept.aoe.vt.edu/~marchman/software/) that will determine the lift and induced drag coefficients of the wing (I know that this will not take account of the wing sweep or dihedral). I will later add to this to include a weight estimate of the wing that should use sweep as a variable. This forms a fitness function for a multiobjective optimization. I would like to set bounds and constraints on some of the variables (not all of them) but I am having difficulty in applying the example given by Matlab to this context - gamultiobj does not seem to accept a constraint function. Examples of these would be setting a minimum wing area and ensuring that Lift = Weight = (1/2 * density * velocity^2 * Wing Area * Lift Coefficient) or (1/2 * density * velocity^2 * Wing Area * Lift Coefficient)-Weight >= 0.
I would also appreciate any other feedback on how to write a good optimization routine as I am new to Matlab and coding in general.
So to summarize my questions:
1) How do I specify equality, inequality and bounds in a multiobjective optimization? 2) How could I improve my optimization routine?
Any help would be greatly appreciated.
Script calling gamultiobj:
FitnessFunction = @LLT_test;
nvar = 18;
% Bounds
LB = [];
UB = [];
A = [];
b = []; Aeq = [];
beq = [];
options = gaoptimset('PopulationSize',100);
[x, fval] = gamultiobj(FitnessFunction, nvar, A,b,Aeq,beq,LB,UB, options);
Fitness Function:
function y = LLT_test(x)
writetofile=0;
% No. of Design Points in Lifting Line Algorithm
n=30;
% Wing Geometry
% x(1) = Root Chord (i)
% x(2) = Root Chord (ii)
% x(3) = Tip Chord (i)
% x(4) = Tip Chord (ii)
% x(5) = Panel Length (i)
% x(6) = Panel Length (ii)
% x(7) = y-coord of root (i)
% x(8) = y-coord of root (ii)
% x(9) = y-coord of tip (i)
% x(10) = y-coord of tip (ii)
% x(11) = Area of panel (i)
% x(12) = Area of panel (ii)
% x(13) = Half chord sweep (i)
% x(14) = Half chord sweep (ii)
% x(15) = Dihedral (i)
% x(16) = Dihedral (ii)
% x(17) = root incidence
% x(18) = tip twist (relative to root)
% y(1) = Wingspan
% y(2) = Aspect Ratio
% y(3) = Lift Coefficient
% y(4) = Induced Drag Coefficient
% y(5) = Effective Wingspan
% y(6) = Effective Aspect Ratio
% y(7) = Wing Area
% y(8) = Effective Wing Area
% Wing Span and Effective Wing Span
y(1) = x(5) * cos(x(15)) + x(6) * cos(x(16));
if x(15)==0 && x(16)==0
else
end
% Taper Ratio
Lambda = (x(5) / (y(5)/2)) * (1 + x(3)/x(1)) + (x(6)/(y(5)/2) * (x(3) + x(4)) / x(3)) - 1;
% Wing Area and Effective Wing Area
y(7) = x(11) * cos(x(15)) + x(12) * cos(x(16));
y(8) = x(11) + x(12);
% Aspect Ratio and effective Aspect Ratio
y(2) = y(1)^2/y(7);
y(6) = y(5)^2 / y(8);
% Aerofoil Properties
a=6.9; % Lift Curve Slope [1/rad]
al0= -6.56; % Zero Lift Angle of Attack [deg]
% Lifting Line Algorithm
% Reference: Marchman, J.F. "MONOPLEQN" [Matlab Code] Available from:
odd=1:2:2*n-1;odd=odd';
phi=pi/(2*n):pi/(2*n):pi/2;phi=phi';
al=(x(17)+x(18)*cos(phi)-al0)*pi/180;
ch=x(1)-(x(1)-x(4))*cos(phi);
bn=ch.*al.*sin(phi)*a/(4*y(1));
for j=1:n
A(:,j)=sin(odd(j)*phi).*(sin(phi)+odd(j)*ch*a/(4*y(1)));
end
x=inv(A)*bn;
y(3) = pi*y(2)*x(1);
y(4)=odd'*(x.*x)*pi/y(2);
EL_Cdi=y(3)^2/(pi*y(2));
m=10;
lc=zeros(m,1);
phi2=(0:(pi/2)/(m-1):pi/2)';
odd2=(1:2:2*n-1)';
lc=x'*sin(odd2*phi2');
lc=4*y(1)*lc'./(x(1)-(x(1)-x(4)).*cos(phi2));
disp(sprintf('Aspect Ratio =%2d',y(2)));
disp(sprintf('Taper Ratio =%2d',Lambda));
disp(sprintf('Root Incidence =%2d',x(17)));
disp(sprintf('Washout =%2d\n',-x(18)));
for i=1:length(phi2)
disp(sprintf(' %2d %8.5f %8.5f %8.5f',floor(phi2(i)*180/pi),lc(i),cos(phi2(i)),lc(i)/y(3)));
end
disp([' Cl =',num2str(y(3))]);
disp([' Cdi =',num2str(y(4))]);
disp(['Elliptical Loading Cdi =',num2str(EL_Cdi)]);
if writetofile fprintf(fid,'Aspect Ratio =%2d\n',y(2));
fprintf(fid,'Taper Ratio =%2d\n',Lambda);
fprintf(fid,'Root Incidence =%2d\n',x(17));
fprintf(fid,'Washout =%2d\n\n',-x(18));
for i=1:length(phi2)
fprintf(fid,' %2d %8.5f %8.5f %8.5f\n',floor(phi2(i)*180/pi),lc(i),cos(phi2(i)),lc(i)/y(3));
end
fprintf(fid,' Cl = %10.7f\n',y(3));
fprintf(fid,' Cdi = %10.7f\n',y(4));
fprintf(fid,'Elliptical Loading Cdi = %10.7f\n',EL_Cdi);
fclose(fid);
disp(['Created filename',blanks(1),filename]);
end
