单纯形法求解函数极值问题 matlab代码

最近整理以前的代码，将以前老师上课的作业代码重新整理，分享出来，作业的内容是编写单纯形法，对测试函数进行寻优(极大值或者极小值)。

首先介绍一下单纯形法：将上课的ppt转化为图片。ppt蓝色背景，眼睛快看瞎了

按照ppt的描述编写算法如下：

clear all;
clc;
% mode可以选择测试函数
% mode = 'exp_test';  
% mode = 'Schaffer';
mode = 'Rosen_Brock';
if strcmp(mode,'Schaffer')
    x=-4:0.1:4;
    y=-4:0.1:4;
    [X,Y]=meshgrid(x,y);
    Z = 0.5-((sin(sqrt(X.^2+Y.^2)).^2)-0.5)./(1+0.001.*(X.^2+Y.^2)).^2;
    figure(1);
    surf(X,Y,Z);
    title('Schaffer Function');
    xlabel('X-轴');
    ylabel('Y-轴');
    zlabel('Z-轴');
end
 
if strcmp(mode,'exp_test')
    x=-4:0.1:4;
    y=-4:0.1:4;
    [X,Y]=meshgrid(x,y);
%     Z = 1400-(20*(X.^2-Y.^2).^2-(1-Y).^2-3*(1+Y).^2+0.3);
%     Z = 5*sin(X.*Y)+X.^2+Y.^2;
    Z = X.*exp(-X.^2-Y.^2)+1;
     figure(1);
    surf(X,Y,Z);
    title('exp_test Function');
    xlabel('X-轴');
    ylabel('Y-轴');
    zlabel('Z-轴');
end
 
if strcmp(mode,'Rosen_Brock')
    x=-2.048:0.1:2.048;
    y=-2.048:0.1:2.048;
    [X,Y]=meshgrid(x,y);
%     Z = 1400-(20*(X.^2-Y.^2).^2-(1-Y).^2-3*(1+Y).^2+0.3);
%     Z = 5*sin(X.*Y)+X.^2+Y.^2;
%     Z = X.*exp(-X.^2-Y.^2)+1;
    Z = 100.*(Y-X.^2).^2+(1-X).^2;
    figure(1);
    surf(X,Y,Z);
    title('Rosen Brock valley Function');
    xlabel('X-轴');
    ylabel('Y-轴');
    zlabel('Z-轴');
end
 
if strcmp(mode,'Rosen_Brock')
    variable_min_x = -2.048;   % 变量的上下限
    variable_max_x = 2.048;
    variable_min_y = -2.048;   % 变量的上下限
    variable_max_y = 2.048;
else
    variable_min_x = -4;   % 变量的上下限
    variable_max_x = 4;
    variable_min_y = -4;   % 变量的上下限
    variable_max_y = 4;
end
 
% variable_min_x = -4;   % 变量的上下限
% variable_max_x = 4;
% variable_min_y = -4;   % 变量的上下限
% variable_max_y = 4;
 
yita = 0.5;
alpha = 1.3;
precision = 10^(-3);
% Schaffer function
% x=variable_min_x:0.1:variable_max_x;
% y=variable_min_y:0.1:variable_max_y;
% [X,Y]=meshgrid(x,y);
% Z = 0.5-((sin(sqrt(X.^2+Y.^2)).^2)-0.5)./(1+0.001.*(X.^2+Y.^2)).^2;
% surf(X,Y,Z);
 
% 复合型法寻找schafer functiong 最大值
 
%% 产生初始的顶点 sige
% M=4;
% N=3;
x=zeros(4,2);
for i=1:4
    x_1 = variable_min_x+rand*(variable_max_x-variable_min_x);
    x_2 = variable_min_y+rand*(variable_max_y-variable_min_y);
    x(i,:)=[x_1,x_2];
end
max_iteraton = 250; % 最大的迭代次数
%%
% 判断索取点是否满足要求：
 index_error = 1;
 error = [];   
 optimization_path = [];
 
for gen=1:max_iteraton
     % 判断点是否在可行域内
     alpha = 1.3;
     for i=1:4
        if ((x(i,1)<variable_min_x)||(x(i,1)>variable_max_x))||((x(i,2)<variable_min_y)||(x(i,2)>variable_max_y))
            error(index_error) = i;  % 不满足的点
            index_error = index_error+1;
        end
     end
     %  error中存放不满足条件的点的索引值
     if isempty(error) == 0  %  error 非空
         % 对不满足的点进行处理
        norm = [];  %  满足约束条件的点
        abnorm = [];  % 不满足约束条件的点
        %%%%%%%%%%%-----
        len_error = length(error);
        for i=1:len_error
            abnorm(i,:) = x(error(i),:);   % 坏点
        end
        kk = 1;
        for m = 1:4
            %kk = kk+1;
            for n = 1:len_error
              if x(m,:) == abnorm(n,:)
                  break;
              end
              if n >= len_error
                  norm(kk,:) = x(m,:);
                  kk = kk+1;
              end
            end
        end
        % 将不满足条件的点进行调整
        % 正常点的形星
        norm_center = 0;
        for m = 1:size(norm,1)
            norm_center = norm_center + norm(m,:);
        end
        norm_center = norm_center/(size(norm,1));
        for j = 1:size(abnorm,1)
            abnorm_new(j,:) = norm_center + yita*(abnorm(j,:)-norm_center);
        end
        %  组成新的点
        x = [norm;abnorm_new];
         
     end 
    %  计算函数值
    for k = 1:4
        f(k) = test_function(x(k,:),mode);
    end
    [f_min,index_1] = min(f);   %  寻找最好的点
    [f_max,index_2] = max(f);
    %  最差的点X_h     次差的点X_g   最好的点X_l
    X_h = x(index_1,:);   
    X_l = x(index_2,:);
    cnt = 1;
    left_points = [];
    left_two_points = [];
    for i=1:4
        if (i == index_1)||(i == index_2)
            continue;
        else
            left_points(cnt) = f(i);    % 除过最好的点和最差的点
            left_two_points(cnt,:) = x(i,:);
            cnt = cnt + 1;
        end
    end
    [left_points_min,index_3] = min(left_points);
    X_g = left_two_points(index_3,:);
    % 求形星：
    cnt = 1;
    for i=1:4
        if i==index_1
            continue;
        else
            points_except_worst(cnt,:) = x(i,:);
            cnt = cnt + 1;
        end
    end
    % 计算形星
    sum = 0;
    for i = 1:length(points_except_worst)
        sum = sum + points_except_worst(i,:);
    end
    X_c = sum/length(points_except_worst);
    %  判断形星是否在可行域内
    if ((X_c(1)>variable_min_x)&&(X_c(1)<variable_max_x))&&((X_c(2)>variable_min_y)&&(X_c(2)<variable_max_y))  % 在可行域内
        % 求反射点
        X_r = X_c + alpha*(X_c - X_h);
        while 1
            if ((X_r(1)>variable_min_x)&&(X_r(1)<variable_max_x))&&((X_r(2)>variable_min_y)&&(X_r(2)<variable_max_y))  % 反射点在可行域内
                break;
            else
                if alpha < precision
                    break;
                end
                alpha = alpha*0.5;
                X_r = X_c + alpha*(X_c - X_h);
            end
        end
    else   % 形星不在可行域内
        % 重新确定上下界
        if X_l(1) <= X_c(1)
            variable_min_x = X_l(1);
            variable_max_x = X_c(1);
        else
            variable_min_x = X_c(1);
            variable_max_x = X_1(1);
        end
        
        if X_l(2) <= X_c(2)
            variable_min_y = X_l(2);
            variable_max_y = X_c(2);
        else
            variable_min_y = X_c(2);
            variable_max_y = X_l(2);
        end
        %   产生新的复合型
        for i=1:4
            x_1 = variable_min_x+rand*(variable_max_x-variable_min_x);
            x_2 = variable_min_y+rand*(variable_max_y-variable_min_y);
            x(i,:)=[x_1,x_2];
        end
        for k = 1:4
            f(k) = Schaffer(x(k,:));
        end
        [f_min,index_1] = min(f);    %  最差的点
        [f_max,index_2] = max(f);    %  寻找最好的点
        %  最差的点X_h     次差的点X_g   最好的点X_l
        X_h = x(index_1,:);   
        X_l = x(index_2,:);
        
        cnt = 1;
        left_two_points = [];
        for i=1:4
           if (i == index_1)||(i == index_2)
               continue;
           else
               left_points(cnt) = f(i);    % 除过最好的点和最差的点
               left_two_points(cnt,:) = x(i,:);
               cnt = cnt + 1;
           end
        end
        [left_points_min,index_3] = min(left_points);
        X_g = left_two_points(index_3,:);
        % 求形星：
        cnt = 1;
        for i=1:4
           if i==index_1
               continue;
           else
               points_except_worst(cnt,:) = x(i,:);
               cnt = cnt + 1;
           end
        end
        % 由新的点产生的形星
        sum = 0;
        for i = 1:length(points_except_worst)
           sum = sum + points_except_worst(i,:);
        end
        X_c = sum/length(points_except_worst);
        % 求反射点
        X_r = X_c + alpha*(X_c - X_h);
        while 1
            if ((X_r(1)>variable_min_x)&&(X_r(1)<variable_max_x))&&((X_r(2)>variable_min_y)&&(X_r(2)<variable_max_y))  % 反射点在可行域内
                break;
            else
                if alpha < precision
                    break;
                end
                alpha = alpha*0.5;
                X_r = X_c + alpha*(X_c - X_h);
            end
        end
    % 反射点计算结束
    end
    %   x X-r X_h X_l X_g 
    f_X_r = test_function(X_r,mode);
    flag_1 = 0;
    if f_X_r > f_min
        x(index_1,:) = X_r;   % 替换最差的点
    else
        % 重新计算反射点
        while test_function(X_r,mode) < f_min
            alpha = alpha*0.5;
            X_r = X_c + alpha*(X_c - X_h);
            if alpha < precision
                flag_1 = 1;
                break;
            end
        end
        %  重新替换  
        if flag_1 == 0   % 新的反射点满足约束条件
            x(index_1,:) = X_r;
        else     %  alpha 减小到足够的范围任然没有找到合适的点
            % 使复合型向最好的点收缩   再重新计算一次新的复合型
            %Nop();
            %disp('nono');
            x(index_1,:) = X_g;
        end
    end
    %  x 更新完毕    如果想画成凸多边形，要处理x
    figure(2);
    for i=1:4-1
        plot([x(i,1),x(i+1,1)], [x(i,2),x(i+1,2)]);
        hold on;
    end
    plot([x(4,1),x(1,1)],[x(4,2),x(1,2)]);
    hold off;
    pause(0.01);
    disp([gen, alpha]);
    function_value(gen) = test_function(X_l,mode);
    max_value = test_function(X_l,mode);
    optimization_path(gen,:) = X_l;
end 
figure(3);
plot(function_value);
title(['函数值变化曲线',' 最大值：',num2str(max_value)]);
xlabel('迭代次数');
ylabel('函数值');
 
figure(4);
plot(optimization_path(:,1),optimization_path(:,2));
title('寻优轨迹');
xlabel('最优点-X');
ylabel('最优点-Y');
 
 
% function f = Schaffer(buf)
%    f = 0.5-((sin(sqrt(buf(1).^2+buf(2).^2)).^2)-0.5)./(1+0.001.*(buf(1).^2+buf(2).^2)).^2;
% end
 
% function f = test_function(buf,mode)
%    if strcmp(mode,'Schaffer')
%        f = 0.5-((sin(sqrt(buf(1).^2+buf(2).^2)).^2)-0.5)./(1+0.001.*(buf(1).^2+buf(2).^2)).^2;
%    end
%    
% %    if mode=='Rastrigin'
% %        f = (buf(1).^2-10*cos(2*pi*buf(1))+10) + (buf(2).^2-10*cos(2.*pi.*buf(2))+10);
% %    end
%    
%    if strcmp(mode,'exp_test')
%        f = buf(1)*exp(-buf(1)^2-buf(2)^2)+1;
%    end
%    
%    if strcmp(mode,'Rosen_Brock')
%         f = 100*(buf(2)-buf(1).^2).^2+(1-buf(1)).^2;     
%    end
%    
% end

test_function.m文件

function f = test_function(buf,mode)
   if strcmp(mode,'Schaffer')
       f = 0.5-((sin(sqrt(buf(1).^2+buf(2).^2)).^2)-0.5)./(1+0.001.*(buf(1).^2+buf(2).^2)).^2;
   end
   
%    if mode=='Rastrigin'
%        f = (buf(1).^2-10*cos(2*pi*buf(1))+10) + (buf(2).^2-10*cos(2.*pi.*buf(2))+10);
%    end
   
   if strcmp(mode,'exp_test')
       f = buf(1)*exp(-buf(1)^2-buf(2)^2)+1;
   end
   
   if strcmp(mode,'Rosen_Brock')
        f = 100*(buf(2)-buf(1).^2).^2+(1-buf(1)).^2;     
   end
   
end

代码运行结果：
1.Schaffer函数（函数是不是很有特点呢）

寻优的值：

 

2. Rosen_Brock函数（是不是更有特点呢）

寻优的值：

寻优点在x-y平面上的寻优轨迹：

从寻优的结果发现，算法陷入了局部最优而没有跳出来，这是一大问题。

2018-12-2  周日。重度雾霾，闲来无事

------------------------------------------------------------end----------------------------------------------------------------------

期待后续有版本更新