% 使用模擬退火法求函數 f(x,y)=5sin(xy)+ x^2 + y^2 的最小值, 

% 根據題意,我們設計冷卻表進度表為:
% 初溫度為100
% 衰減參數為0.95
% 馬可夫鏈長度為10000
% 使用Metropolis接受準則進行模擬,Metropolis的步長為0.02
% 結束條件為根據上一個最優解與最新的一個最優解的之差小于某個容差.
function r=SimulateAnnealing()
  format long
  tic
  %搜索的最大區間
  XMAX=4;
  YMAX=4;
  
  %冷卻表參數
  MarkovLength=10000;   %馬可夫鏈長度
  DecayScale=0.95;      %衰減參數
  StepFactor=0.02;      %步長因子
  Temperature=100;      %初始溫度
  Tolerance=1e-8;       %容差
  
  %PreX,NextX  分別為prior and next value of x
  %PreY,NextY  分別為prior and next value of y
  %PreBextX,PreBestY 是上一個最優解
  %BestX,BestY   是最終解
  AcceptPoints=0;     %Metropolis過程中總接受點
  
  %隨機選點
  PreX=-XMAX*rand;
  PreY=-YMAX*rand;
  PreBestX=PreX;
  PreBestY=PreY;
  BestX=PreX;
  BestY=PreY;
  
  %每迭代一次退火一次(降溫),直到滿足迭代條件為止
 %Decay=Temperature*DecayScale;
  while(1)
     Temperature=Temperature*DecayScale;
     AcceptPoints=0;
     
     %當前溫度下迭代loop(即MARKO鏈長度)次
     for i=1:MarkovLength
         
        % 1) 在此點附近隨機選下一點
        while(1)
            NextX=PreX+StepFactor*XMAX*(rand-0.5);
            NextY=PreY+StepFactor*YMAX*(rand-0.5);
            if (NextX>=-XMAX & NextX<=XMAX & NextY>=-YMAX & NextY<=YMAX)
                break;
            end
        end    

        % 2) 是否全局最優解
        if ObjectFunction(BestX,BestY)>ObjectFunction(NextX,NextY)
            %保留上一個最優解
            PreBestX=BestX;
            PreBestY=BestY;
            %此為新的最優解
            BestX=NextX;
            BestY=NextY;
        end              
        
        % 3) Metropolis過程              
        if ObjectFunction(PreX,PreY)-ObjectFunction(NextX,NextY)>0
            % 接受,此處lastPoint即下一個迭代的點以新接受的點開始
            PreX=NextX;
            PreY=NextY;
            AcceptPoints=AcceptPoints+1;
        else
            change=-1*(ObjectFunction(NextX,NextY)-ObjectFunction(PreX,PreY))/Temperature;
            if exp(change)>rand
                PreX=NextX;
                PreY=NextY;
                AcceptPoints=AcceptPoints+1;
            end %不接受,保存原解
        end   
     end
     %符合終止條件退出
     if (abs(ObjectFunction(PreBestX,PreBestY)-ObjectFunction(BestX,BestY))<Tolerance)
            BestX
            BestY
            BestF=ObjectFunction(BestX,BestY)
            break;
     end
  end    
  toc
  
  %graph
  
  %[X,Y]=meshgrid(-4:4,-4:4);
  %Z=5*sin(X*Y)+X^2+Y^2;
  
  ezmesh('f')
  grid on
  hold on
  plot3(BestX,BestY,BestF,'*r')
  hold off
  
 
  figure;clf;
  grid on
  ezsurf('f',[-4,4])
  hold on
  plot3(BestX,BestY,BestF,'*r')
  hold off
  
function f=ObjectFunction(x,y)
f=5*sin(x*y)+x^2+y^2;