simplex.m

遗传算法与单纯形算法复合用于光纤与波导芯片对准仿真技术分析

%Simplex method
%对于光纤对准问题，是寻找函数最大点为优化目标。
%通常最优化目标是寻找函数最小值，所以程序要注意针对这一点的改动。
%Ref. Book(Su Tashan 2001) P97
%fun:name of the function f(x)
%x0:nx1 vector of initial point
%option=[beita,gama,eps]'
%beita:0<beita<1
%gama: gama>1
%eps:precision control
%x:returned maximum x
%五维单纯形法，即自变量只包含五个变量
%simplex([3 2 3 4 20]',[0.5 2 0.999]',[6 -4 1])
function [x,counts,y,xi]=simplex(x0,option,scale)
if nargin<=1 || isempty(option)
    beita=0.5;                  %缩减因子,0.5
    gama=2;                     %扩张因子,2
    eps=0.99;                  %判别条件
    lateral=1;                    %横向调整边长
    gap=-4;                       %纵向调整边长
    ang=1;                        %角度调整边长
    %lateral=6;
else
    beita=option(1);        %缩减因子,0.5
    gama=option(2);         %扩张因子,2
    eps=option(3);          %判别条件
    lateral=scale(1);              %单纯形边长,%横向调整边长
    gap=scale(2);                       %纵向调整边长
    ang=scale(3);                        %角度调整边长 
end

counts=0;
fhandle = str2func('funn');
n=length(x0);
s=x0;
xi(1,:)=x0';                %记录自变量的变化
y(1,1)=feval(fhandle,x0);   %初始点函数值
                            %构造初始单纯形，初始点再加n个点
                            %对于光纤对准问题，线性与角度值，横向与纵向取值范围差别甚大，所以如何构造单纯形尚需仔细考虑！
for j=2:3                   %横向
    s(:,j)=x0;              %第j列元素都等于初始值x0
    s(j-1,j)=s(j-1,j)+lateral;    %第j列第j-1行元素加lateral length(单纯形边长)
end
for j=4:5                   %角度
    s(:,j)=x0;              %第j列元素都等于初始值x0
    s(j-1,j)=s(j-1,j)+ang;  %第j列第j-1行元素加lateral length(单纯形边长)
end
 j=6;                       %纵向
 s(:,j)=x0;                 %第j列元素都等于初始值x0
 s(j-1,j)=s(j-1,j)+gap;     %第j列第j-1行元素加lateral length(单纯形边长)
while 1
    counts=counts+1;        %迭代次数计数 
    for j=1:n+1             %求n+1个顶点函数值
        f(j,1)=feval(fhandle,s(:,j));
    end
    fabs=abs(f);            %选择Xl,Xh,Xg,Xh对应最好点,Xl对应最坏点,Xg对应次坏点。 不管有多少个顶点，都只选“最好，最坏，次坏”三个点。        
    [minf,l]=min(fabs);     %l对应最小函数值，在光纤对准问题中应为最坏点。
    [maxf,h]=max(fabs);     %h标号对应函数最大值，在光纤对准问题中应为最好点。
    
    xh=s(:,h);
    xl=s(:,l);
    
    xi(counts+1,:)=xh';                 %每次迭代后的自变量（最好点）
    y(counts+1,1)=feval(fhandle,xh);    %每次迭代后的函数值
            
    fabs(l)=[];             %把向量中l对应的元素删掉，后面的元素前移。即：删掉当前最小，再求剩余元素最小，即得“次坏点”。
    [bigerf,g]=min(fabs);   %g标号对应中间值，在光纤对准问题中应为次坏点。
    xg=s(:,g);
    average=zeros(n,1);     %准备计算形心与误差
    for j=1:n+1
        if j~=l             %计算形心时不包括最坏点。
            average=average+s(:,j);%顶点相加求平均
        end
    end
    average=average/n;      %计算形心与误差
    %error=sqrt(sum((f-feval(fhandle,average)*ones(n+1,1)).^2)/(n+1));%ones语句表示从1到n+1累加。该判别终止条件不适合用于光纤对准数学模型。
    %不如用自变量最小来做终止条件
    %为了目标的一致性，还是以目标函数值即耦合效率达到0.999为判别终止条件,at the same time,the endless iterations should be prevented.
    goal=y(counts+1,1);
    if goal>eps || counts>500
        x=s(:,h);
        break;
    else
        xr=2*average-s(:,l);    %开始迭代，xr为反射点,注意，是两倍均值减最坏点。反射系数取为1。
        fr=feval(fhandle,xr);   %fr为反射点xr对应的函数值
        if fr>maxf              %第一种情况，反射点xr比最好点Xh还好，趁机扩张，得到新点
            xe=average+gama*(xr-average);%计算扩张点
            fe=feval(fhandle,xe);
            if fe>fr          %如果扩张点比反射点还好，则用它取代最坏点构成新的单纯形。重要！
                s(:,l)=xe;
            else                %否则，扩张失败，用反射点取代最坏点构成新的单纯形。
                s(:,l)=xr;
            end
        else
            if fr>bigerf        %第二种情况，反射点不优于“最好点”，但优于“次坏点”。
                s(:,l)=xr;      %用反射点取代最坏点构成新的单纯形。
            else
                if fr>minf      %第三种情况，反射点不优于“次坏点”，但优于“最坏点”。
                    xc=average+beita*(xr-average);%计算收缩点Xc
                    fc=feval(fhandle,xc);
                else            %第四种情况，反射点比“最坏点”Xl还要坏。
                    xc=average+beita*(s(:,l)-average);%计算收缩点Xc，计算公式与第三种情况稍异。
                    fc=feval(fhandle,xc);
                end
                if fc>minf      %如果压缩点优于“最坏点”Xl，则以压缩点取代最坏点构成新的单纯形。
                    s(:,l)=xc;
                else            %否则，将当前单纯形各顶点向“最好点”Xh紧缩,即缩边。
                    xh=s(:,h);
                    for j=1:n+1                        
                        s(:,j)=xh+0.5*(s(:,j)-xh);
                    end
                end
            end
        end
    end
end                 

figure(2);
plot([0:1:counts]',y,'b');xlabel('steps');ylabel('normalized coupling efficiency');grid on;
%plot([0:1:counts]',y,'s');xlabel('steps');ylabel('normalized coupling efficiency');grid on;hold off;
figure(3);
plot([0:1:counts]',xi(:,1),'b');xlabel('steps');ylabel('x(um)');grid on;hold on;
plot([0:1:counts]',xi(:,1),'o');xlabel('steps');ylabel('x(um)');grid on;hold off;
figure(4);
plot([0:1:counts]',xi(:,2),'b');xlabel('steps');ylabel('y(um)');grid on;hold on;
plot([0:1:counts]',xi(:,2),'o');xlabel('steps');ylabel('y(um)');grid on;hold off;
figure(5);
plot([0:1:counts]',xi(:,3),'o');xlabel('steps');ylabel('pitch(deg)');hold on;grid on;
plot([0:1:counts]',xi(:,3),'b');xlabel('steps');ylabel('pitch(deg)');hold off;
figure(6);
plot([0:1:counts]',xi(:,4),'o');xlabel('steps');ylabel('yaw(deg)');hold on;grid on;
plot([0:1:counts]',xi(:,4),'b');xlabel('steps');ylabel('yaw(deg)');hold off;
figure(7);
plot([0:1:counts]',xi(:,5),'o');xlabel('steps');ylabel('z(um)');hold on;grid on;
plot([0:1:counts]',xi(:,5),'b');xlabel('steps');ylabel('z(um)');hold off;
figure(8);
plot(xi(:,1),xi(:,2),'b');xlabel('x(um)');ylabel('y(um)');grid on;%hold on;
%plot(xi,'s');xlabel('x(um)');ylabel('y(um)');grid on;hold off;