%使用模擬退火法求函數f(x,y)=5sin(x*y)+x^2+y^2的最小值
format long
XMAX=4; %搜索的最大區間
YMAX=4; %搜索的最大區間
MarkovLength=10000; %馬可夫鏈長度
DecayScale=0.95; %衰減參數0.95
StepFactor=0.02; %步長因子
Temperature=100; %初始溫度
Tolerance=1e-8; %容差
AcceptPoints=0.0; %Metropolis過程中總接受點
PreX=-XMAX*rand; %初始的搜索值
PreY=-YMAX*rand; %初始的搜索值
PreBestX=PreX; %上一個最優解
PreBestY=PreY; %上一個最優解
BestX=PreX; %最終解
BestY=PreY; %最終解
while(1)
   Temperature=Temperature*DecayScale; %每迭代一次退火一次（降溫），直到滿足迭代條件為止
   AcceptPoints=0.0;
   %在當前溫度下迭代（即MARKOV鏈長度）次
   for i=0:1:MarkovLength 
   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
   %判斷新產生點與原來最優點哪個更優
       if(minfunction(BestX,BestY)>minfunction(NextX,NextY))
          PreBestX=BestX; %保留上一個最優解
          PreBestY=BestY;
          BestX=NextX; %新的最優解
          BestY=NextY;
       end
       %接受新產生的點為下一迭代的開始點
       if(minfunction(PreX,PreY)-minfunction(NextX,NextY)>0)
          PreX=NextX;
          PreY=NextY;
          AcceptPoints=AcceptPoints+1;
       else
          change=-1*(minfunction(NextX,NextY)-minfunction(PreX,PreY))/Temperature;
          %不接受，保存原解
          if(exp(change)>rand)
             PreX=NextX;
             PreY=NextY;
             AcceptPoints=AcceptPoints+1;
          end
       end
   end
   %結束條件為根據上一個最優解與最新的一個最優解的之差小于某個容差
   if(~(abs(minfunction(BestX,BestY)-minfunction(PreBestX,PreBestY))>Tolerance))
   break
       
   end
end
   a=BestX
   b=BestY
   c=minfunction(BestX,BestY)
   %經測試，已成功運行a=1.0767,b=-1.0767,c=-2.2640