How do I specify equality, inequality and bounds in a multiobjective optimization and how could I improve this optimization?

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

    y(5) = y(1);

else

    y(5) = x(5) + x(6);

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:

% http://www.dept.aoe.vt.edu/~marchman/software/

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

Answer 1

[Alan Weiss]

If you use fgoalattain instead of gamultiobj then you can use all constraint types. The goal attainment formulation is, perhaps, a bit obscure, but it is a way to perform multiobjective optimization.

However, depending on the code that calculates the objective function, you might be in the situation where small changes in parameters do not necessarily give a change in the objective function. If so, use the suggestions in Optimizing a Simulation or ODE.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

Answer 2

[HACRS4]

Thank you, I will try using fgoalattain.