用粒子群解决有约束的最优解问题

  前一段师姐帮学姐用 粒子群算法计算带线性等式约束的最优问题。matlab自带的粒子群算法并没有这个功能，所以对原有的粒子群增加了两个步骤来进行线性等式约束。

  没理解基础粒子群算法的童鞋请先去别的文章处弄懂，我是在这里看的 推荐一下  https://blog.csdn.net/nightmare_dimple/article/details/74331679

  增加的两个步骤。

  1. 引入惩罚因子，在离线性等式较远的地方让他适应度（函数值）无穷大（我找的是最小值）。

  2.速度分情况给值，如果不在那条线附近（线性等式约束），那么他的速度只是群体学习速度。

我觉得这种思想（惩罚因子和速度选择）也是可以解决非线性约束的，有需要的童鞋可以验证一下。

clear
% 初始化种群
f_1 = @(x)0.01011.*x.^2-1.49559.*x+230.38805;
f_2 = @(x)0.0103.*x.^2-1.52432.*x+233.45944;
f_min = @(x,y)2*f_1(x/16800)*x+2*f_2(y/9600)*y;

tic;
N = 5000;                         % 初始种群个数 粒子数为1000个
d = 2;                          % 空间维数 
threshold =45800;             %   x1+x2+x3+x4≥45800 的阈值（线性不等式约束）
ger = 180;                      % 最大迭代次数    （我这个点比较多，迭代次数不用那么多） 
limit = [14000, 15220;7700,8700];      % X 的上下限
           
vlimit = [- 7.5, 7.5;-7.5, 7.5];          % 种群移动速度的上下限
c_1 = 0.8;                        % 惯性权重
c_2 = 0.5;                       % 自我学习因子
c_3 = 0.5;                       % 群体学习因子 
 for i = 1:d
    x(:,i) = limit(i, 1) + (limit(i, 2) - limit(i, 1)) * rand(N, 1);%初始种群的位置
end    
v = rand(N, d);                  % 初始种群的速度
xm = x;                          % 每个个体的历史最佳位置
ym = zeros(1, d);                % 种群的历史最佳位置
fxm = zeros(N, 1)+inf;               % 每个个体的历史最佳适应度
fym = inf;                      % 种群历史最佳适应度

iter = 1;                   % 初始迭代次数为一
times = 1;                   % 进度条计数用的
bar = waitbar(0,'读取数据中...');           % 查看进度用的进度面板
ezplot('(96/168)*x',[10000,16000]);
  hold on
while iter <= ger
     str=['计算中...',num2str(100*iter/ger),'%'];    % 百分比形式显示处理进程
     waitbar(iter/ger,bar,str)                       % 更新进度条bar，配合bar使用
     for i = 1:N
        fx(i) = f_min(x(i,1),x(i,2)) ; % 个体当前适应度  
            if ((2*x(i,1)+2*x(i,2))<threshold)         % 若不满足x1+x2+x3+x4≥45800要求则惩罚因子无穷大
                  fx(i) = inf;
            end
            if((abs(9600*x(i,1)-16800*x(i,2)))>=500)            % 若不满足|9600*x1-16800*x2|>=500 则惩罚因子无穷大
                fx(i) = inf;
            end
    end
      
     for i = 1:N      
        if (fxm(i) > fx(i))         
            fxm(i) = fx(i);     % 更新个体历史最佳适应度(即：看谁的f_min更小，并把每个粒子历史上最小的f_min存起来）
            xm(i,:) = x(i,:);   % 更新个体历史最佳位置
        end 
     end
     if fym > min(fxm)
        [fym, nmin] = min(fxm);   % 更新群体历史最佳适应度 找最小值    
        ym = xm(nmin, :);      % 更新群体历史最佳位置
     end
     
  for i = 1:N
    if(((2*x(i,1)+2*x(i,2))<threshold))
        v(i,:) =  c_3 * rand *(ym - x(i,:));  %如果不满足x1+x2+x3+x4≥45800，则速度方向取决于群体学习
    elseif((9600*x(i,1)~=16800*x(i,2)))        %如果不满足|9600*x1-16800*x2|>=500，则速度方向取决于群体学习
         v(i,:) =  c_3 * rand *(ym - x(i,:));
    else
         v(i,:) = v(i,:) * c_1 + c_2 * rand *(xm(i,:) - x(i,:)) + c_3 * rand *(ym - x(i,:));% 速度更新
    end
  end
    % 边界速度处理
    for i=1:d 
        for j=1:N
        if  v(j,i)>vlimit(i,2)
            v(j,i)=vlimit(i,2);
        end
        if  v(j,i) < vlimit(i,1)
            v(j,i)=vlimit(i,1);
        end
        end
    end      
    if iter % 100 == 0 
        
     ezplot('(96/168)*x',[10000,16000]);
     hold on
    end
    x = x + v;% 位置更新
    for i=1:d            % 边界位置处理
        for j=1:N
        if  x(j,i)>limit(i,2)
            x(j,i)=limit(i,2);
        end
        if  x(j,i) < limit(i,1)
            x(j,i)=limit(i,1);
        end
        end
    end     
  
    if times >= 10              %这是画动态图用的
        cla;
        scatter(x(:,1),x(:,2));title('状态位置变化');
        xlabel(" X1 ");
        ylabel(" X3 ");
        pause(0.5);           % 暂停0.5秒
        times=0;
        
    end
   iter = iter+1;
    times=times+1;
end


toc;