A New Global Optimization Method based on Clustering

by

 

The clustering technique and parabolic approximation are used in this method.

GOBC_PA.m
function [Best_point,Best_result,perf,broken_epoch,performance]=GOBC_PA(func,rang,target)
% GOBC-PA:  A new unconstrained global optimization method based on clustering
%           and parabolic approximation. GOBC-PA gets the
%           minimum values of objective function using with parameter populations at each step.
%           To determine the best samples of the population, the samples are clustered. It is assumed
%           that the cluster centers can be located nearby the minima, maxima or between
%           of them. In this case, the global optimum closes to the cluster centers. If some samples
%           of the next generation are created randomly around the cluster centers, it is possible
%           to reach the global optimum.
%
% Arguments: (input)
% rang         - range of search space
% target       - target function value. That value used for measuring error.
% func         - objective function as text
% Arguments: (output)
% Best_point   - Global minimum points
% Best_result  - Global minimum value on global minimum points
% perf         - errors
% broken_epoch - epoch number of code break
% performance  - best points and their objective value on each step
%
% Example usage:
% (Sixhump_Camelback function)
%
% str='4*(x(1).^2)-2.1*(x(1).^4)+(1/3)*(x(1).^6)+(x(1).*x(2))-4*(x(2).^2)+4*(x(2).^4)';
% str=['@(x) (',str,')'];
% test=str2func(str);
% range=[-5 5];
% target=-1.0316285;
%
% [Best_point,Best_result,err,iteration,performance]=GOBC_PA(test,range,target);
%
% Author: Bayram Cetisli and Ihsan Pene
% e-mail address: bayramcetisli@sdu.edu.tr ihsanpence@mehmetakif.edu.tr
% Release date: 07/07/13
warning off;
if(nargin==2)
    target=0;
    disp('Error value can not be measure because of missing target value');
end

param_n=Calc_dimension(func); % Estimate dimension number of objective function
epoch=1000;                   % Epoch size
p_size=60;                    % Population size
limit=1e-8;                   % Parameter for breaking code
kume=ceil(p_size/3);          % Number of cluster size
size_kum=ceil(kume/3);        % Number of best clusters size
size_pop=2;                   % Number of best population size for creating new populations around them
ma=rang(2);                   % max range
mi=rang(1);                   % min range

popu=mi+(ma-mi)*rand(p_size,param_n);   % First population created randomly

performance=zeros(epoch,param_n+1);
[evali]=Calc_func_values(popu,func);    % Estimate populations objective values
[zz,vv]=sort(evali);
broken_epoch=1000;
for i=1:epoch
    [t,c]=fcm([popu evali],kume,[2;2;1e-5;0]);      % two step Fuzzy C-Means clustering
    t(:,end)=Calc_func_values(t(:,1:end-1),func);   % Estimate populations objective values
    [zz,kk]=sort(t(:,end));
    [zz,mm]=max(c);
    n=size_kum;
    for ii=1:n
        [zz,rr]=find(mm==kk(ii));
        if isempty(rr)==1
            [s,zz]=find(kk==kk(ii));
            kk(s)=[];
            kk(end+1)=kk(ii);
        end
    end
    tr_temp=[];
    for ik=1:p_size
        [zz,bk]=max(c(:,ik));
        tr_temp=[tr_temp bk]; % Finding the cluster names for each populations
    end
    pop_new=[];
    
    for j=1:size_kum
        % Parabolic approximation starts
        x_yer= (tr_temp==kk(j))';
        if sum(x_yer)>2
            A=[popu(x_yer,:).^2 popu(x_yer,:) ones(sum(x_yer==1),1)];
            param=A\evali(x_yer);
            tepe=[-param(param_n+1:end-1)./(2*param(1:param_n))]';
            ind=param(param_n+1:end-1)==0;
            tepe(ind)=0;
            ind=sign(param(1:param_n))<0;
            tepe(ind)=t(kk(j),ind);
            if sum(isnan(tepe))==0 && sum(isinf(tepe))==0
                t(kk(j),1:end-1)=tepe;
            end
        end
        % Parabolic approximation ends
        kume_eleman_std=std(popu(x_yer,:)); % Estimate the standart derivation of populations which are members of clusters
        if(size(kume_eleman_std,2)==1)      % If the cluster has one member fix the standart derivation
            kume_eleman_std=ones(1,param_n);
        end
        pop_new=[pop_new;((1/i)*[kume_eleman_std;kume_eleman_std;kume_eleman_std].*(randc(3,param_n)*2))+[t(kk(j),1:end-1);t(kk(j),1:end-1);t(kk(j),1:end-1)]];
    end
    pop_best=[];
    for j=1:size_pop
        pp_best_temp=[((max(popu)-min(popu))/(2*kume));((max(popu)-min(popu))/(2*kume));((max(popu)-min(popu))/(2*kume))];
        pop_best=[pop_best;((pp_best_temp).*(randc(3,param_n)*2))+[popu(vv(j),:);popu(vv(j),:);popu(vv(j),:)]];
    end
    
    popu=[popu(vv(1),:);...
        popu(vv(2),:);...
        t(kk(1:size_kum),1:param_n);...
        pop_new;...
        pop_best;...
        ];
    popu=[popu;rang(1)+(rang(2)-rang(1))*rand(p_size-size(popu,1),param_n)];
    popu(end,:)=round(popu(vv(1),:));

    ind=popu>ma;    % constrained the population in search space
    popu(ind)=ma;
    ind=popu<mi;
    popu(ind)=mi;
    
    ind=sum(isnan(popu),2)>0;
    ind=logical(ind*ones(1,param_n));
    if sum(sum(ind))>0
        popu(ind)=rang(1)+(rang(2)-rang(1))*rand(1);
    end
    [evali]=Calc_func_values(popu,func);
    [zz,vv]=sort(evali);
    
    performance(i,:)=[popu(vv(1),:) evali(vv(1))];
    perf(i,1)=abs(target-evali(vv(1),1));
    if(perf(i,1)<limit) % Breaking code
        broken_epoch=i;
        break;
    end
end
disp('Best points');
disp(popu(vv(1),:));
disp('Best result for best points');
disp(evali(vv(1),1));
disp('target');
disp(target);
Best_point=popu(vv(1),:);
Best_result=evali(vv(1),1);
disp('error');
disp(perf(i));
broken_epoch

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [output]=Calc_func_values(popu,func)
for i=1:size(popu,1)
    output(i,1)=func(popu(i,:));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function sonuc=Calc_dimension(test)
a=func2str(test);
v=strfind(a,'x(');
rr=[];
for i=1:size(v,2)
    if(str2num(a(v(i)+2))>=0)
        if(str2num(a(v(i)+3))>=0)
            rr=[rr a(v(i)+2) a(v(i)+3) ' '];
        else
            rr=[rr a(v(i)+2) ' '];
        end
    end
end
sonuc=max(str2num(rr));

Contact us