%%%%%%%%% Nine Problems are solved here %%%
%%%%%%%%
%% Q1 Minmization of unconstrained fitness function using
GA MATLAB
% By Dr Endalew Ayenew
% Objective function
% array x and returns the objective function value
function f= EvalObj(x)
f= x (1) ^2+x (2) ^2+x (3) ^2; % objective function
end
% Main optimization function using the GA
rng default % Initialize the random number generator
options= optimoptions('ga','HybridFcn',@fmincon);
options.PlotFcn= 'gaplotbestf';
options.PopulationSize= 30;
options.MaxGenerations= 100;
nvar= 3;
LB= [0 0 0];
UB= [15 15 15];
[x,fval]= ga(@EvalObj,nvar,[],[],[],[],LB,UB,[],options)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%
%% Q2 Minmization of constrained fitness function using GA
% by Dr Endalew Ayenew
% objective function
function y=myFitness(x)
y= 100*(x(1)^2-x(2))^2+(1-x(1))^2
end
% Constraint function
function [C, Ceq]= myConstraints(x)
C= [1.5+x(1)*x(2)+x(1)-x(2);10-x(1)*x(2)]; %linear equality
constraints:
Ceq=[];% No nonlinear equality constraints:
end
% main code for minmization of the constrained fitness function
by GA
ObjFcn=@myFitness;
nVars=2;
LB=[0 0];
UB=[1 13];
ConsFun=@myConstraints;
Optionsoptimoptions('ga','MutationFcn',@mutationadaptfeasibl
e);
[x,fval]=ga(ObjFcn, nVars, [],[],[],[],LB,UB,ConsFun,options);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%
%% Q3 Optimization of multiobjective functions using GA by
Dr Endalew Ayenew
% multiobjective optimization using GA
function y =simplemultiobjective(x)
% Initialize for two objectives
y = zeros(2,1);
% Compute first objective
for i = 1:2
y(1) = y(1) - 10*exp((x(i)^2 + x(i+1)^2));
end
% Compute second objective
for i = 1:3
y(2) = y(2) + abs(x(i))^0.8 + 5*cos(x(i)^3);
end
% main function for optimization of multiobjective functions by
GA
FitnessFunction = @simplemultiobjective; % Function handle
to the fitness function
nVars = 3; % Number of decision variables
lb = [-5 -5 -5]; % Lower bound
ub = [5 5 5]; % Upper bound
A = []; % No linear inequality constraints
b = []; % No linear inequality constraints
Aeq = []; % No linear equality constraints
beq = []; % No linear equality constraints
options = optimoptions(@gamultiobj,'PlotFcn',@gaplotpareto);
[x,Fval,exitFlag,Output] = gamultiobj(FitnessFunction,nVars,A,
...
b,Aeq,beq,lb,ub,options)
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
%% Q4 pid controller parameters tuning using
Optimization toolbox GA solver
%% By Dr Endalew Ayenew
%objective function
function [J]=pid_optim(x)
% plant
n=[1];
d =[1 10 20];
plant=tf(n,d);
%controller
kp=x(1)
ki=x(2)
kd=x(3)
s=tf('s');
cont=kp+ki/s+kd*s;
step(feedback(plant*cont,1));
dt=0.1; t=0:dt:1;
% cost function formation=ITAE
e=1-step(feedback(plant*cont,1),t);
J=sum(t'.*abs(e)*dt)
%% goto optimization toolbox and select ga in solver and set all
other parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%
%% Q5 Unconstrained optimization problem is solved by using
PSO MATLAB code
% array x and returns the objective function value
% Dr Endalew Ayenew
function f= EvalObj(x)
f= x(1)^2+x(2)^2+x(3)^2; % objective function
end
% PSO main program
rng default % Initialize the random seed
options= optimoptions('particleswarm','HybridFcn',@fmincon);
options.PlotFcn= 'pswplotbestf';
options.SwarmSize= 50;
options.MaxIterations= 100;
nvar= 3; % Number of decision variables
LB= [-5 -5 -5]; % Lower bounds
UB= [1000 200 100]; % Upper bounds
[x,fval]= particleswarm(@EvalObj,nvar,LB,UB,options)
%%’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’’%
%%%%
%% Q6 minimization of constrained objective function by PSO
%% By Dr Endalew Ayenew
function f = confun(x)
c0= [];
% minimization of objective function
of= 10*(x(1)-1)^2+20*(x(2)-2)^2+30*(x(3)-3)^2;
% constraints (all constraints must be converted into,= 0
type)
% if there is no constraints then comments all c0 lines below
c0(1) = x(1)+x(2)+x(3)-5; % <=0 type constaint
c0(2) = x(1)^2+2*x(2)-x(3); % <=0 type constaint
%define penalty for each constraint
w= [500 20]; % penalty weight on each constraint vaiolation
p= 0;
for i = 1:length(c0)
if c0(i)>0
p= p+w(i)*abs(c0(i)); % penalty
end
end
f = of + p; % objective function + penalty
% main PSO program for optimization of constrained problem
tic
clear all
close all
LB= [0 0 0]; % Lower bound
UB= [10 10 10];% upper bound
VLB= -abs(UB-LB); VUB= abs(UB-LB); % Lower/Upper
bound of velocity
% PSO parameters values
m = 3; % no. of variables or swarms
n= 100; % population size
wmax = 0.9; % max inertia
wmin = 0.4; % min inertia
c1 = 2; c2 = 2; % acceleration factors
% PSO main program
maxite = 1000; % max no. of iteration
maxrun = 1; % max no. of runs or populations
for run = 1:maxrun
% PSO intialization
for i = 1:n % max no. of runs or populations
for j = 1:m % no. of variables
x0(i,j) = LB(j)+rand*(UB(j)-LB(j));
end
end
x = x0; % initial population
v = 0.1*x0; % initial velocity
for i = 1:n
f0(i)= confun(x0(i,:));
end
[fmin0,index0]=min(f0);
pbest =x0; % initial personel best
gbest =x0(index0,:); % initial global best
% PSO Algorithm stated
ite = 1;
tolerance =1;
while ite <= maxite && tolerance> 10^-12
w = wmax-(wmax-wmin)*ite/maxite; % update inertia
% PSO velocity updates
for i = 1:n
for j = 1:m
v(i,j) = w*v(i,j)+c1*rand()*(pbest(i,j)-x(i,j))
+c2*rand()*(gbest(j)-x(i,j));
v(i,j)= min(max(VLB(j),v(i,j)),VUB(j)); % velocity
clamping
end
end
% PSO position update
for i = 1:n
for j = 1:m
x(i,j) = x(i,j)+v(i,j);
x(i,j)= min(max(LB(j),x(i,j)),UB(j)); % position
clamping
end
end
% evaluating fitness
for i = 1:n
f(i) = confun(x(i,:));
end
% updateing pbest and fitness
for i = 1:n
if f(i)<f0(i)
pbest(i,:)= x(i,:);
f0(i)= f(i);
end
end
[fmin,index]=min(f0); % Finding out the best particle
ffmin(ite, run)=fmin; % storing best fitness
ffite(run)=ite; % storing iteration count
% updating gbest and best fitness
if fmin<fmin0
gbest = pbest(index,:);
fmin0 =fmin;
end
% calculating tolerance
if ite>100
tolerance = abs(ffmin(ite-100,run)-fmin0);
end
%displaying relative results
if ite == 1
fprintf('Iteration Best particle Objective function\n');
end
fprintf('%8g %8g %8.4f\n', ite, index, fmin0);
ite = ite+1;
end
% PSO algorithm end
gbest;
fvalue=10*(gbest(1)-1)^2+20*(gbest(2)-2)^2+30*(gbest(3)-
3)^2;
fff(run)=fvalue;
rgbest(run,:)=gbest;
fprintf('-----------------------------------------\n');
end
% PSO main program is ended
fprintf('\n');
fprintf('*********************************************
*\n');
fprintf('Final result\n');
[bestfun, bestrun]=min(fff)
best_variables=rgbest(bestrun,:)
fprintf('*********************************************
*\n');
toc
% PSO Convergence characteritics
plot(ffmin(1:ffite(bestrun),bestrun),'-k');
xlabel('Fitness function value');
title('PSO Convergence characteritics')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
%% Q7 Maximization of constrained objective function by
PSO
% by Dr Endalew Ayenew
function f = confun(x)
c0= [];
% Q6 is changed to maximization
% Maximization of objective function
of= 10*(x(1)-1)^2+20*(x(2)-2)^2+30*(x(3)-3)^2;
% constraints (all constraints must be converted into ,= 0
type)
% if there is no constraints then comments all c0 lines below
c0(1) = x(1)+x(2)+x(3)-5; % <=0 type constaint
c0(2) = x(1)^2+2*x(2)-x(3); % <=0 type constaint
%define penalty for each constraint
w= [500 20]; % penalty weight on each constraint vaiolation
p= 0;
for i = 1:length(c0)
if c0(i)>0
p= p+w(i)*abs(c0(i)); % penalty
end
end
f = of + p; % objective function + penalty
end
% main PSO program for maximization of constrained problem
tic
clear all
close all
LB= [0 0 0]; % Lower bound
UB= [10 10 10];% upper bound
VLB= -abs(UB-LB); VUB= abs(UB-LB); % Lower/Upper
bound of velocity
% PSO parameters values
m = 3; % no. of variables or swarms
n= 100; % population size
wmax = 0.9; % max inertia
wmin = 0.4; % min inertia
c1 = 2; c2 = 2; % acceleration factors
% PSO main program
maxite = 1000; % max no. of iteration
maxrun = 1; % max no. of runs or populations
for run = 1:maxrun
% PSO intialization
for i = 1:n % max no. of runs or populations
for j = 1:m % no. of variables
x0(i,j) = LB(j)+rand*(UB(j)-LB(j));
end
end
x = x0; % initial population
v = 0.1*x0; % initial velocity
for i = 1:n
f0(i)= confun(x0(i,:));
end
[fmin0,index0]=min(f0);
pbest =x0; % initial personel best
gbest =x0(index0,:); % initial global best
% PSO Algorithm stated
ite = 1;
tolerance =1;
while ite <= maxite && tolerance> 10^-12
w = wmax-(wmax-wmin)*ite/maxite; % update inertia
% PSO velocity updates
for i = 1:n
for j = 1:m
v(i,j) = w*v(i,j)+c1*rand()*(pbest(i,j)-x(i,j))
+c2*rand()*(gbest(j)-x(i,j));
v(i,j)= min(max(VLB(j),v(i,j)),VUB(j)); % velocity
clamping
end
end
% PSO position update
for i = 1:n
for j = 1:m
x(i,j) = x(i,j)+v(i,j);
x(i,j)= min(max(LB(j),x(i,j)),UB(j)); % position
clamping
end
end
% evaluating fitness
for i = 1:n
f(i) = confun(x(i,:));
end
% updateing pbest and fitness
for i = 1:n
if f(i)>f0(i) % < for minimization is changed to > for
maximization
pbest(i,:)= x(i,:);
f0(i)= f(i);
end
end
[fmin,index]=min(f0); % Finding out the best particle
ffmin(ite, run)=fmin; % storing best fitness
ffite(run)=ite; % storing iteration count
% updating gbest and best fitness
if fmin>fmin0 % < for minimization is changed to > for
maximization
gbest = pbest(index,:);
fmin0 =fmin;
end
% calculating tolerance
if ite>100
tolerance = abs(ffmin(ite-100,run)-fmin0);
end
%displaying relative results
if ite == 1
fprintf('Iteration Best particle Objective function\n');
end
fprintf('%8g %8g %8.4f\n', ite, index, fmin0);
ite = ite+1;
end
% PSO algorithm end
gbest;
fvalue=10*(gbest(1)-1)^2+20*(gbest(2)-2)^2+30*(gbest(3)-
3)^2;
fff(run)=fvalue;
rgbest(run,:)=gbest;
fprintf('-----------------------------------------\n');
%end
end
% PSO main program is ended
fprintf('\n');
fprintf('*********************************************
*\n');
fprintf('Final result\n');
[bestfun, bestrun]=min(fff)
best_variables=rgbest(bestrun,:)
fprintf('*********************************************
*\n');
toc
% PSO Convergence characteritics
plot(ffmin(1:ffite(bestrun),bestrun),'-k');
xlabel('Fitness function value');
title('PSO Convergence characteritics')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%
%% Q8 Tune the gains of the PID controller using the PSO by
Dr Endalew Ayenew
%% Objective function for the PID controller parameters tuning
by PSO
function J= EvalObj_PSOPID(x)
% Retrieve PID controller parameters
Kp= x(1); Ki= x(2); Kd= x(3);
% Plant
s= tf('s');
Gc= Kp+Ki/s+Kd*s; % PID controller TF
K= 1; T= 5; L= 2;
Gp= K*exp(-L*s)/(T*s+1); % Plant TF
% closed loop TF
G= feedback(Gp*Gc,1); % Overall TF from yr to y
% Unit step response for 40seconds
h= 0.01; t= 0:h:40; r= 1;
y= step(G,t);
% Calculate IAE using the trapezoidal integration rule
e= abs(r-y); % IAE
% e= (r-y).^2; % ISE
% e= t'.*abs(r-y); % ITAE (t must be transposed)
J= 0.5*h*(e(1)+2*sum(e(2:end-1))+e(end));
end
% main function
rng default % Initialize the random number generator
options = optimoptions('particleswarm');
options.PlotFcn= 'pswplotbestf';
options.MaxIterations= 100; % default 200 x nvar
nvar= 3;
LB= [0 0 0];
UB= [10 10 10];
[x,fval]=
particleswarm(@EvalObj_PSOPID,nvar,LB,UB,options)
% Simulate the PID control system tuned by PSO
clf
Kp= 2.2560; Ki= 0.3856; Kd= 1.8355;
s= tf('s'); K= 1; T= 5; L= 2; h= 0.01;
Gc= Kp+Ki/s+Kd*s;
Gp= K*exp(-L*s)/(T*s+1);
G= Gc*Gp/(1+Gc*Gp);
t= 0:h:40;
figure
y= step(G,t);
plot(t,y,'linewidth',2), axis([0 40 0 1.5])
xlabel('t'), ylabel('y(t)')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%
%% Q9 ACO tuned PID controller parameters by Dr Endalew
Ayenew
clc; clear all; close all
% plant
n=[1];
d=[1 2 3];
G=tf(n,d);
Gf = feedback(G,1);
step(Gf) % System step response with out controller
hold on
% ACO Parameters
N=10; ph_decay=0.6; sca_fact=2; itr=10; var=3; h_size=5;
%search space
a=0; b=1000; GBv=10^12; c_cf=0;
for cc=1:var
% initial value
vv=a;
to_ph_mt(cc)=0;
% assign permissible disccrete values & amounts 0f pheromone
for m=1:b
x(cc,m)=vv;
vv=vv+h_size;
ph_mt(cc,m)=1;
to_ph_mt(cc)= to_ph_mt(cc)+ph_mt(cc,m);
end
end
for tt=1:itr
for cc=1:var
% probability of selected path
for m =1:b
prop(cc,m)=ph_mt(cc,m)/to_ph_mt(cc);
end
%Roulette Wheel selection
for m=1:b
if m==1
rws(cc,m)=0;
rwe(cc,m)=rws(cc,m)+prop(cc,m);
else
rws(cc,m)= rwe(cc,m-1);
rwe(cc,m)=rws(cc,m)+prop(cc,m);
end
end
for k=1:N
% select random position for each ant
rr(cc,k)=rand(1);
% put each ant on the Roulette Wheel
for m=1:b
if (rr(cc,k)> rws(cc,m))&&(rr(cc,k)< rwe(cc,m))
ff(cc,k)=x(cc,m);
end
end
end
end
% objective function values corresponding to the paths
chosen
for k=1:N
%Model parameters
kp=ff(1,k);
ki=ff(2,k);
kd=ff(3,k);
% simulate model
Gc= pid(kp,ki,kd);
Gcf = feedback(Gc*G,1);
y= step(Gcf);
% ITAE
ff_cc=0;
sim2=size(y);
for m_err = 1:sim2;
ff_cc=ff_cc+((abs(y(m_err)-1))*m_err);
end
ITAE(k)=ff_cc;
end
% Rank objective function
[BFv, BFl]= min(ITAE);
WFv =max(ITAE);
if BFv < GBv
GBv = BFv;
for cc = 1:var
bkc(cc)=ff(cc,BFl);
end
end
c_cf = c_cf+1;
best_cf_ac(c_cf)=GBv;
% updat the pheromone array
up_ph =(sca_fact*BFv)/WFv;
for cc = 1:var
to_ph_mt(cc)=0;
for m=1:b
if x(cc,m)==bkc(cc)
ph_mt(cc,m)=ph_mt(cc,m)+up_ph;
else
ph_mt(cc,m)=ph_mt(cc,m)*ph_decay;
end
to_ph_mt=to_ph_mt+ph_mt(m);
end
end
end
min_ITAE =GBv
kp=bkc(1)
ki=bkc(2)
kd=bkc(3)
Gc= pid(kp,ki,kd);
Gcf = feedback(Gc*G,1);
step(Gcf); % System step response with optimized PID
controller
t_cf = 1:c_cf
figure
plot(t_cf, best_cf_ac,'r--','LineWidth',2)
xlabel('iteration')
ylabel('cost function (ITAR)')
legend('ITAE for ACO-PID')
