Code covered by the BSD License

# Harmony Search Algorithm

27 Sep 2010 (Updated )

Meta-heuristic Optimization Method

[BestGen,BestFitness,gx]=HarmonySearch
```function [BestGen,BestFitness,gx]=HarmonySearch

% This code has been written with Matlab 7.0
% You can modify the simple constraint handlening method using more efficient
% methods. A good review of these methods can be found in:
% Carlos A. Coello Coello,"Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: a survey of the state of the art"

clc
clear;
global NVAR NG NH MaxItr HMS HMCR PARmin PARmax bwmin bwmax;
global HM NCHV fitness PVB BW gx;
global BestIndex WorstIndex BestFit WorstFit BestGen currentIteration;

NVAR=4;         %number of variables
NG=6;           %number of ineguality constraints
NH=0;           %number of eguality constraints
MaxItr=5000;    % maximum number of iterations
HMS=6;          % harmony memory size
HMCR=0.9;       % harmony consideration rate  0< HMCR <1
PARmin=0.4;      % minumum pitch adjusting rate
PARmax=0.9;      % maximum pitch adjusting rate
bwmin=0.0001;    % minumum bandwidth
bwmax=1.0;      % maxiumum bandwidth
PVB=[1.0 4;0.6 2;40 80;20 60];   % range of variables

% /**** Initiate Matrix ****/
HM=zeros(HMS,NVAR);
NCHV=zeros(1,NVAR);
BestGen=zeros(1,NVAR);
fitness=zeros(1,HMS);
BW=zeros(1,NVAR);
gx=zeros(1,NG);
% warning off MATLAB:m_warning_end_without_block

MainHarmony;

% /**********************************************/

function sum =Fitness(sol)

sum = 0.6224*sol(1)*sol(3)*sol(4)+1.7781*sol(2)*sol(3)^2+3.1661*sol(1)^2*sol(4)+19.84*sol(1)^2*sol(3)+ eg(sol);  %F(x) = f(x) + penalty

end

% /*********************************************/
function sum=eg(sol)

% constraints g(x) > 0
gx(1)=sol(1)-0.0193*sol(3);     %  x1 - 0.0193 x3 > 0
gx(2)=sol(2)-0.00954*sol(3);
gx(3)=3.14*sol(3)^2*sol(4)+(4/3)*3.14*sol(3)^3 - 1296000;
gx(4)=-sol(4)+240;
gx(5)=sol(1) - 1.1;
gx(6)=sol(2) - 0.6;

% we use static penalty function to handle constraints
sum = 0;
for i=1:NG
if(gx(i)<0)
sum = sum - 1000 * gx(i);
end
end
end

% /*********************************************/

function initialize
% randomly initialize the HM
for i=1:HMS
for j=1:NVAR
HM(i,j)=randval(PVB(j,1),PVB(j,2));
end
fitness(i) = Fitness(HM(i,:));
end
end

%/*******************************************/

function MainHarmony
% global NVAR NG NH MaxItr HMS HMCR PARmin PARmax bwmin bwmax;
% global HM NCHV fitness PVB BW gx currentIteration;

initialize;
currentIteration  = 0;

while(StopCondition(currentIteration))

PAR=(PARmax-PARmin)/(MaxItr)*currentIteration+PARmin;
coef=log(bwmin/bwmax)/MaxItr;
for pp =1:NVAR
BW(pp)=bwmax*exp(coef*currentIteration);
end
% improvise a new harmony vector
for i =1:NVAR
ran = rand(1);
if( ran < HMCR ) % memory consideration
index = randint(1,HMS);
NCHV(i) = HM(index,i);
pvbRan = rand(1);
if( pvbRan < PAR) % pitch adjusting
pvbRan1 = rand(1);
result = NCHV(i);
if( pvbRan1 < 0.5)
result =result+  rand(1) * BW(i);
if( result < PVB(i,2))
NCHV(i) = result;
end
else
result =result- rand(1) * BW(i);
if( result > PVB(i,1))
NCHV(i) = result;
end
end
end
else
NCHV(i) = randval( PVB(i,1), PVB(i,2) ); % random selection
end
end
newFitness = Fitness(NCHV);
UpdateHM( newFitness );

currentIteration=currentIteration+1;
end
BestFitness = min(fitness);
end
% /*****************************************/

function UpdateHM( NewFit )
% global NVAR MaxItr HMS ;
% global HM NCHV BestGen fitness ;
% global BestIndex WorstIndex BestFit WorstFit currentIteration;

if(currentIteration==0)
BestFit=fitness(1);
for i = 1:HMS
if( fitness(i) < BestFit )
BestFit = fitness(i);
BestIndex =i;
end
end

WorstFit=fitness(1);
for i = 1:HMS
if( fitness(i) > WorstFit )
WorstFit = fitness(i);
WorstIndex =i;
end
end
end
if (NewFit< WorstFit)

if( NewFit < BestFit )
HM(WorstIndex,:)=NCHV;
BestGen=NCHV;
fitness(WorstIndex)=NewFit;
BestIndex=WorstIndex;
else
HM(WorstIndex,:)=NCHV;
fitness(WorstIndex)=NewFit;
end

WorstFit=fitness(1);
WorstIndex =1;
for i = 1:HMS
if( fitness(i) > WorstFit )
WorstFit = fitness(i);
WorstIndex =i;
end
end

end
end % main if
end %function

% /*****************************************/
function val1=randval(Maxv,Minv)
val1=rand(1)*(Maxv-Minv)+Minv;
end

function val2=randint(Maxv,Minv)
val2=round(rand(1)*(Maxv-Minv)+Minv);
end
% /*******************************************/

function val=StopCondition(Itr)
global MaxItr;
val=1;
if(Itr>MaxItr)
val=0;
end
end

% /*******************************************/

```