amebsa.m

本程序实际上是把模拟退火过程引入单纯形法来求多元函数的极值

function [p,iter]=amebsa(p,y,ndim,pb,ftol,temptr,tt,idum,funk)
%用模拟退火与Nelder和Mead的下降单纯形法相结合的方法求多元函数funk(x)的极小值，
%其中x(1...ndim)为一个ndim维向量。作为输入参数矩阵p（1...ndim+1,1...ndim)共有
%ndim+1行，每行均为一个ndim维向量，分别代表初始单纯形的各个顶点。amebsa的输入项
%还包括矢量y(1...ndim+1)、浮点数ftol和temptr，其中y的各个元素应当被初始化为函数
%funk在p的ndim＋1个顶点的值；ftol为函数值所要达到的相对收敛容许限，一旦满足要求，
%程序将尽早返回；iter和temptr的意义分别是函数值的计算次数和温度。程序在某退火
%温度temptr处进行iter次函数值计算后返回。
%接下来的工作是根据退火进程调低温度temptr，重置iter并再次调用该程序（每次调用时
%其他参数保持不变）。如果iter以正值返回，则说明程序正常收敛。如果第一次调用时yb
%被初始化为一个很大的数，则yb和pb(1...dim)将依次返回已遇到的最佳函数值和最佳点。
global tt;
global idum;
global funk;
funk=@RGDMI2;
psum=[];
tt=-temptr;
mpts=ndim+1;
for n=1:ndim   %置psum
    sum=0.0;
    for m=1:mpts
        sum=sum+p(m,n);
    end
    psum(n)=sum;
end
LOGR=tt*log(ran1(idum))
while(ndim~=0)
    ilo=1;  %确定最高点（即差点）、次高点和最低点（即最佳点）
    lhi=2;
    ylo=y(1)+LOGR;%我们所“看到”的顶点总是处于随机的热起伏状态
    ynhi=ylo;
    yhi=y(2)+LOGR;
    if ylo>yhi
        ihi=1;
        ilo=2;
        ynhi=yhi;
        yhi=ylo;
        ylo=ynhi;
    end
    for i=3:mpts  %对单纯形中的点进行循环
         yt=y(i)+LOGR; %更多的热起伏运动
        if yt<=ylo
            ilo=i;
            ylo=yt;
        end
        if yt>yhi
            ynhi=yhi;
            ihi=i;
            yhi=yt;
        end
        if yt>ynhi
            ynhi=yt;
        end
    end
    rtol=2.0*abs(yhi-ylo)/(abs(yhi)+abs(ylo));
    %计算从最高点到最低点的范围，若合乎要求，则返回
    if rtol<ftol|iter<0 %若返回，将最佳点和最佳值放入槽1中
        swap=y(1);
        y(1)=y(ilo);
        y(ilo)=swap;
        for n=1:ndim
            swap=p(1,n);
            p(1,n)=p(ilo,n);
            p(ilo,n)=swap;
        end
        break;
    end
    iter=iter-2;
    %开始新一轮的迭代，首先从高点通过单纯形的表面，以因子-1做外插
    %即从高点对单纯形进行反射
    [ytry,yhi,pb]=amotsa(p,yb,ihi,yhi,ndim,psum,-1.0);  
    if ytry<=ylo
        [ytry,yhi,pb]=amotsa(p,yb,ihi,yhi,ndim,psum,2.0); 
        %得到了比最佳点还要好的结果，因此以因子2做外推
    elseif ytry>=ynhi
        %反射后的点不如次高点，因此需要找一个中间的次低点，做一次一维收缩
        ysave=yhi;
        [ytry,yhi,pb]=amotsa(p,yb,ihi,yhi,ndim,psum,0.5);
        if ytry>=ysave %似乎无法摆脱高点，最好围绕最低点也即最佳点进行收缩
            for i=1:mpts
                if i~=ilo
                    for j=1:ndim
                        psum(j)=0.5*(p(i,j)+p(ilo,j));
                        p(i,j)=psum(j);
                    end
                    y(i)=funk(psum);
                end
            end
            iter=iter-nim;
            for n=1:ndim    %重新计算psum
                 sum=0.0;
                 for m=1:mpts
                      sum=sum+p(m,n);
                 end
                 psum(n)=sum;
            end 
        end
    else
        iter=iter+1;  %纠正计数器
    end
end