funopt.m

用蚂蚁群算法 编写的函数优化 matlab程序

function [ObjFunVal,FunVarVal]=FunOpt()
global ASOption
%%%%%%%%%%%%% the key parameters of Ant System %%%%%%%%%
MaxITime=1e3;
alpha=1;
rho=0.65;
AntNum=20;
VarDim=2;
BlaCat=11;
%%%%%%%%%%%%% the key parameters of Ant System %%%%%%%%%
FunIter=1;
N=10;
Tolerance=0.01;
delt=5;
    VarLowUpp=InitOptFun();  
    h(FunIter,:)=EnCode(VarLowUpp,N,FunIter);
    tempTol0=h(FunIter,:);
    [GBTour,GBObjFun,Option,IBRecord]=AS(AntNum,VarDim,N+1,MaxITime,alpha,rho,VarLowUpp,tempTol0);
    ObjFunVal(FunIter,1)=GBObjFun;
    FunVarVal(FunIter,:)=DeCode(GBTour',h(FunIter,:),VarLowUpp(1,:)); 
while 1      
    FunIter=FunIter+1;
    VarLowUpp=InitVarLowUpp(FunVarVal(FunIter-1,:),h(FunIter-1,:),delt);
    h(FunIter,:)=EnCode(VarLowUpp,N,FunIter);
    tempTol=h(FunIter,:);
    [GBTour,GBObjFun,Option,IBRecord]=AS(AntNum,VarDim,N+1,MaxITime,alpha,rho,VarLowUpp,tempTol);
    ObjFunVal(FunIter,1)=GBObjFun;
    FunVarVal(FunIter,:)=DeCode(GBTour',h(FunIter,:),VarLowUpp(1,:));   
    if FunTerminate(h(FunIter,:),Tolerance)
        break;
    end
end

function VarLowUpp=InitVarLowUpp(OptVarVal,Tol,delt)
global ASOption
VarLowUpp=zeros(3,ASOption.n);
VarLowUpp(2,:)=OptVarVal;
VarLowUpp(1,:)=VarLowUpp(2,:)-Tol*delt;
VarLowUpp(3,:)=VarLowUpp(2,:)+Tol*delt;


function Tol=EnCode(VarLowUpp,N,FunIter)
global ASOption
Tol=(VarLowUpp(3,:)-VarLowUpp(1,:))/N;


function VarVal=DeCode(Tours,Tol,VarLow)
global  ASOption
for i=1:ASOption.n
    VarVal(i)=VarLow(1,i)+(Tours(1,i)-1)*Tol(i);
%     VarVal(i)=BlaCat(Tours(i));
end

function  ObjFunVal=CalObjFun(VarVal)
%  cos_in = cos(2*pi*VarVal);
%  ObjFunVal    = sum((VarVal.^2-10*cos_in + 10), 2);
 ObjFunVal    = (VarVal(1)-1).^2+(VarVal(2)-2.2).^2 + 1;
 
function FunVarLowUpp=InitOptFun()
global ASOption
FunVarLowUpp=zeros(3,2);
FunVarLowUpp(1,:)=[0,1];
FunVarLowUpp(3,:)=[2,3];

 
function flag =FunTerminate(Tol,Tolerance)
flag = false;
if  Tolerance>=max(Tol);
    flag = true;
end

function [GBTour,GBObjFun,Option,IBRecord]=AS(AntNum,VarDim,BlaCat,MaxITime,alpha,rho,VarLowUpp,Toler)
%% (Ant System) date:070501 

global ASOption Problem AntSystem
ASOption = InitParameter(AntNum,VarDim,BlaCat,alpha,rho,MaxITime);
Problem = InitProblem(BlaCat,VarDim,VarLowUpp,Toler);
AntSystem = InitAntSystem();
ITime = 0;
IBRecord = [];
if ASOption.DispInterval ~= 0
    close all
    set(gcf,'Doublebuffer','on');
    hline=plot(1,1,'-o');
end
while 1
    InitStartPoint();
    for step = 2:ASOption.n
        for ant = 1:ASOption.m
            [pheromone,P]= CaculateShiftProb(step,ant);
            qrand=rand;
            if (qrand<=ASOption.qo)
                [maxtao,i]=max(pheromone);
                nextnode=i;
            else              
                nextnode = Roulette(P,1);
            end
            RefreshTabu(step,ant,nextnode);
            currentnode=AntSystem.tours(ant,step-1);
            oldtau=Problem.tau(:,:,step);
            newtau=oldtau(currentnode,nextnode)*(1-ASOption.rho)+ASOption.rho*ASOption.C;   
            Problem.tau(:,:,step)=newtau;
        end
    end
    CloseTours()
%     LocalSearch();ObjFun
    
    CaculateToursObjFun();
    GlobleRefreshPheromone();
    ITime = ITime + 1;
    [GBTour,GBObjFun,IBRecord(:,ITime)] = GetResults(ITime);
%     ShowIterativeCourse(GBTour,ITime,hline);
    if Terminate(ITime)
        break;
    end
end
Option = ASOption;
%% --------------------------------------------------------------
function ASOption = InitParameter(AntNum,VarDim,BlaCat,alpha,rho,MaxITime)
ASOption.l =BlaCat;
ASOption.n =VarDim;
ASOption.m =AntNum;
ASOption.alpha = alpha;
ASOption.rho = rho;
ASOption.C = 100;
ASOption.MaxITime = MaxITime;
ASOption.qo = 0.9;
ASOption.Q = 10;
ASOption.GBObjFun = inf;
ASOption.GBTour = zeros(VarDim+1,1);
ASOption.OptITime = 1;
ASOption.DispInterval = 10;
rand('state',sum(100*clock));
%% --------------------------------------------------------------
function Problem = InitProblem(BlaCat,VarDim,VarLowUpp,Toler)
global ASOption
MatrixTau = (ones(BlaCat,BlaCat,VarDim))*ASOption.C;
VarLower=VarLowUpp(1,:);
Problem = struct('VarDim',VarDim,'tau',MatrixTau,'VarLow',VarLower,'Tol',Toler);
%% --------------------------------------------------------------
function AntSystem = InitAntSystem()
global ASOption
AntTours = zeros(ASOption.m,ASOption.n+1); 
ToursObjFun = zeros(ASOption.m,1);
AntSystem = struct('tours',AntTours,'ObjFuns',ToursObjFun);
%% --------------------------------------------------------------
function InitStartPoint()
global AntSystem ASOption
AntSystem.tours = zeros(ASOption.m,ASOption.n+1); 
rand('state',sum(100*clock));
AntSystem.tours(:,1) = randint(ASOption.m,1,[1,ASOption.l]);
AntSystem.ObjFuns = zeros(ASOption.m,1);
%% --------------------------------------------------------------
function [pheromone,Probs]= CaculateShiftProb(step_j, ant_k)
global AntSystem ASOption Problem
currentnode=AntSystem.tours(ant_k,step_j-1);
tau_i = Problem.tau(currentnode,:,step_j);
pheromone=(tau_i.^ASOption.alpha);
Probs = pheromone;
if sum(Probs) ~= 0
    Probs = Probs/sum(Probs);
end
%% --------------------------------------------------------------
function Select = Roulette(P,num) 
m = length(P);
flag = (1-sum(P)<=1e-5);
Select = zeros(1,num);
rand('state',sum(100*clock));
r = rand(1,num);
for i=1:num
    sumP = 0;
    j = ceil(m*rand); 
    while (sumP<r(i)) && flag
        sumP = sumP + P(mod(j-1,m)+1);
        j = j+1;
    end
    Select(i) = mod(j-2,m)+1;
end
%% --------------------------------------------------------------
function RefreshTabu(step_j,ant_k,nextnode)
global AntSystem
AntSystem.tours(ant_k,step_j) = nextnode;
%% --------------------------------------------------------------
function CloseTours()
global AntSystem ASOption
AntSystem.tours(:,ASOption.n+1) = AntSystem.tours(:,1);
%% --------------------------------------------------------------
function CaculateToursObjFun()
global AntSystem ASOption Problem
ObjFuns = zeros(ASOption.m,1);
for k=1:ASOption.m
    ObjFuns(k)=CalObjFun(DeCode(AntSystem.tours(k,:),Problem.VarLow,Problem.Tol));
end
AntSystem.ObjFuns = ObjFuns;
%% --------------------------------------------------------------
function [GBTour,GBObjFun,Record] = GetResults(ITime)
global AntSystem ASOption
[IBObjFun,AntIndex] = min(AntSystem.ObjFuns);
IMObjFun=mean(AntSystem.ObjFuns);
IBTour = AntSystem.tours(AntIndex,:);
if IBObjFun<=ASOption.GBObjFun 
	ASOption.GBObjFun = IBObjFun;
	ASOption.GBTour = IBTour;
	ASOption.OptITime = ITime;
end
GBTour = ASOption.GBTour';
GBObjFun = ASOption.GBObjFun;
Record = [IMObjFun,IBObjFun]';
%% --------------------------------------------------------------

function GlobleRefreshPheromone()
global AntSystem ASOption Problem
[IBObjFun,AntIndex] = min(AntSystem.ObjFuns);
IBTour = AntSystem.tours(AntIndex,:);
sumdtau=zeros(ASOption.l,ASOption.l,ASOption.n);  
    sumdtau(IBTour(ASOption.n),IBTour(1),1)=sumdtau(IBTour(ASOption.n),IBTour(1),1)+ASOption.Q/IBObjFun;   
for i=2:ASOption.n 
    sumdtau(IBTour(i-1),IBTour(i),i)=sumdtau(IBTour(i-1),IBTour(i),i)+ASOption.Q/IBObjFun;       
end
Problem.tau=Problem.tau*(1-ASOption.rho)+ASOption.rho*sumdtau;
%% --------------------------------------------------------------
function flag = Terminate(ITime)
global ASOption
flag = false;
if  ITime>=ASOption.MaxITime
    flag = true;
end
%% --------------------------------------------------------------

function ShowIterativeCourse(IBTour,ITime,hmovie)
global Problem ASOption
if mod(ITime,ASOption.DispInterval)==0
    title(get(hmovie,'Parent'),['ITime = ',num2str(ITime)]);    
    plot(IBTour);
    pause(0.01);
end