Optimization

381Appendix B

yyy = func_val(X);

LAMBDA = LAMBDA + 2*RK*EQCONSTR;

BETA = BETA + 2*RK*(max([ICONSTR; -BETA./(2*RK)]));

fnew = func_val(X);

checkconstr1 = max([ineqconstr_val(X);checkconstr]);

disp([i X FVALUE norm([checkconstr1 eqconstr_val(X)])]);

if abs(fnew-fprev) < epsilon1

break

end

end

fprintf('___________________________________________________\n')

plot(X(1), (1+X(2)^2-X(1)-0.1*sin(3*pi*X(1))),'r*')

hold on

end

xlabel('f1')

ylabel('f2')

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% MATLAB code func_val.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% computes augmented Lagrangian function value

%

function y = func_val(x)

global LAMBDA RK BETA EQCONSTR ICONSTR FVALUE W1 W2

y = (1+x(2)^2-x(1)-0.1*sin(3*pi*x(1)));

g = ineqconstr_val(x);

h = eqconstr_val(x);

EQCONSTR = h;

ICONSTR = g;

FVALUE = y;

y = y + LAMBDA*EQCONSTR' + RK*EQCONSTR*EQCONSTR' + sum(BETA.*

max([ICONSTR; -BETA./(2*RK)])) + sum(RK*(max([ICONSTR; -BETA./

(2*RK)])).^2);

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% MATLAB code ineqconstr_val.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% computes value of inequality constraint

function g = ineqconstr_val(x)

global LAMBDA RK BETA EQCONSTR ICONSTR FVALUE W1 W2

g(1) = x(1)-1;

g(2) = -x(1);

g(3) = x(2)-2;

g(4) = -x(2)-2;

g(5) = x(1)- W1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

382 Appendix B

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% File name pso.m

% Particle Swarm Optimization algorithm with dynamic weights

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

% lb -> lower bound of variables

% ub -> upper bound of variables

% x -> position of individual

% v -> velocity of individual

% rand -> random number from 0 to 1

% fitness -> fitness of individual

% pbest -> best fitness achieved by individual

% gbest -> best fitness of group

%

clear all

clc

format long

global W1 W2

pop = 200;

phi_1 = 0.5;

phi_2 = 0.5;

nmax = 120;

weight = linspace(1,0.4,nmax);

lb = [0 -2];

ub = [1 2];

W1 = 0;

W2 = 1;

for i = 1:length(lb)

for j = 1:pop

x(i,j) = lb(i) + (ub(i)-lb(i))*rand;

v(i,j) = 0;

end

end

for i = 1:pop

fitness(i) = func1(x(:,i));

pbest(i) = fitness(i);

px(i,:) = x(:,i);

end

[gbest, location] = min(fitness);

gx = x(:,location);

for i = 1:nmax

W1 = abs(sin(2*pi*i/150));

W2 = 1-W1;

for j = 1:pop

v(:,j) = weight(i)*v(:,j) + phi_1*rand*(px(j,:)'-x(:,j)) +

phi_2*rand*(gx-x(:,j));

x(:,j) = x(:,j) + v(:,j);

for k = 1:length(x(:,j))

if x(k,j) < lb(k) || x(k,j) > ub(k)

x(k,j) = lb(k) + (ub(k)-lb(k))*rand;

383Appendix B

end

end

fitness(j) = func1(x(:,j));

if fitness(j) < pbest(j)

pbest(j) = fitness(j);

px(j,:) = x(:,j);

end

F1(j) = x(1,j);

F2(j) = (1+x(2,j)^2-x(1,j)-0.1*sin(3*pi*x(1,j)));

end

[gbest, location] = min(pbest);

gx = x(:,location);

[gx' gbest];

plot(F1,F2,'r*')

pause(0.1)

end

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% File name func1.m

% Enter the function to be optimized

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

function [y] = func1(x)

global W1 W2

y = W1*x(1) + W2*(1+x(2)^2-x(1)-0.1*sin(3*pi*x(1)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% MATLAB code sqp.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

% n_of_var -> number of design variables

% x = [-1.5 1.5] -> starting value of x

% epsilon1 -> constant used for terminating the algorithm

% delx -> required for gradient computation

% falpha_prev -> function value at first/previous iteration

% deriv -> gradient vector

% quadprog -> MATLAB function to solve quadratic programming

% LAMBDA -> Lagrange multipliers

clear all

clc

warning off

n_of_var = 5;

n_of_eqcons = 1;

n_of_iqcons = 12;

scale_factor = 1;

384 Appendix B

global LAMBDA RK BETA EQCONSTR ICONSTR FVALUE W1 W2

LAMBDA = zeros(1,n_of_eqcons);

BETA = zeros(1,n_of_iqcons);

X = [0.5 25 31 0.5 0.5];

RK = 1;

A = eye(length(X));

epsilon1 = 1e-2;

delx = 1e-3;

for kk = 640:20:1630

W1 = (kk-1)/100;

checkconstr = zeros(1,n_of_iqcons);

fprintf(' No. x-vector f(x) |Cons.| \n')

fprintf('_________________________________________________\n')

checkconstr = zeros(1,n_of_iqcons);

for i = 1:30

deriv_f = grad_vec_f(X,delx,n_of_var,scale_factor);

sec_deriv_f = hessian(X,delx);

deriv_eqcon = grad_vec_eqcon(X,delx,n_of_eqcons);

deriv_ineqcon = grad_vec_ineqcon(X,delx,n_of_iqcons);

options = optimset('Display','off');

x = quadprog(sec_deriv_f,deriv_f,deriv_ineqcon,-ineqconstr_

val(X),deriv_eqcon,-eqconstr_val(X),[],[],X,options);

fprev = func_val(X);

X = X+x';

yyy = func_val(X);

LAMBDA = LAMBDA + 2*RK*EQCONSTR;

BETA = BETA + 2*RK*(max([ICONSTR; -BETA./(2*RK)]));

fnew = func_val(X);

checkconstr1 = max([ineqconstr_val(X);checkconstr]);

disp([i X FVALUE norm([checkconstr1 eqconstr_val(X)])]);

if abs(fnew-fprev) < epsilon1

break

end

end

fprintf('_________________________________________________\n')

[ xcp, area, volume] = dynamics(X);

plot(xcp, area,'k*','LineWidth',1.5)

hold on

end

xlabel('X_{cp}')

ylabel('A')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% MATLAB code dynamics.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xcp, area, volume] = dynamics(X)

385Appendix B

beta1 = [0.075759311956522

0.001381832173914

0.005825624927536

0.087880851017943

0.079788071083505

0.103099119565217

0.000091411652174

-0.000083773739130

0.024378102714516

0.052872440763746];

cm = -0.278527718840579+beta1'*[X(1);X(2);X(3);X(4);X(5);

X(1)^2;X(2)^2;X(3)^2;X(4)^2;X(5)^2];

beta2 = [0.150130422826087

-0.003965447971014

0.014019523043478

0.103639584627328

0.108001867667355

0.047906228260870

0.000186851275362

-0.000219940956522

-0.000485274327121

0.016298096388315];

cn = -0.314286112399353+beta2’*[X(1);X(2);X(3);X(4);X(5);

X(1)^2;X(2)^2;X(3)^2;X(4)^2;X(5)^2];

xcp = cm/cn;

r1 = X(1)*cos(deg2rad(X(2)));

r2 = X(1)*cos(deg2rad(X(2))) + X(4)*tan(deg2rad(X(2)));

r3 = r2 + X(5)*tan(deg2rad(X(3)));

area = 2*pi*X(1)*X(1)*(1-sin(deg2rad(X(2)))) +

pi*(r1+r2)*sqrt((r2-r1)^2 + X(4)^2) +

pi*(r3+r2)*sqrt((r3-r2)^2 + X(5)^2) + pi*r3^2;

cap = (pi*r1*(1-sin(deg2rad(X(2))))/6) * (3*r1*r1 +

(r1*(1-sin(deg2rad(X(2))))^2));

volume = cap + 0.3333*pi*X(4)*(r2^2 + r1^2 + r2*r1) +

0.3333*pi*X(5)*(r3^2 + r2^2 + r2*r3);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Chapter 9

Code Name Details

sqp.m (main program) Sequential quadratic programming (SQP) method (for

multidisciplinary design optimization [MDO] application)

discipline1.m (function) Output from rst discipline