有约束情况下的粒子群算法函数寻优

对无约束粒子群算法介绍
(统一回复一下，代码链接在此，这个代码和博客有所不同，对有点基础的同学来说问题不大的，祝诸位学业有成~
新链接：https://pan.baidu.com/s/1VUs872FqWeuFbGat7p4-Gg?pwd=4c9m 提取码: 4c9m 复制这段内容后打开百度网盘手机App，操作更方便哦)

约束条件

对目标函数中的变量进行的约束。即一个函数中的各个变量有的会有限制条件，只有在约束条件下进行的求解才能是有效的求解。

以上图为例，f(x) 为目标函数，s.t. 为约束条件

等式约束

例如上图中的h(x) =0 ，我们称之为等式约束。

对于函数 y = x1 +x2*x3

x1 + x2 = x3 就是等式约束

不等式约束

例如上图中g(x) <= 0 ，我们称之为不等式约束。

对于函数 y = x1 +x2*x3

x1 + x2 <= x3 就是等式约束

在算法中引入约束

我们知道，在一般的粒子群算法中，我们并没有引入变量之间约束的概念，所以要把约束条件引入算法，需要将约束条件转化为限制粒子飞行的因素。

在算法执行的过程中，我们通过粒子的适应度（一般为粒子代表的函数值）来代表这个粒子的优势度。如果这个粒子不满足我们引入的约束条件，我们可以给这个粒子的适应度赋予一个很大的值（适应度越小越好），说明这个粒子不能满足约束条件，因此，我们引入了罚函数的概念。

罚函数

百度解释

罚函数是指在求解最优化问题(无线性约束优化及非线性约束优化)时，在原有目标函数中加上一个障碍函数，而得到一个增广目标函数，罚函数的功能是对非可行点或企图穿越边界而逃离可行域的点赋予一个极大的值，即将有约束最优化问题转化为求解无约束最优化问题。

我们将约束条件进行整理，转化为右边为0的形式

以上面为例：x1 + x2 - x3 <= 0

function Z=getconstraints(fnonlin,u)
% 一个很大的数，用于惩罚没有满足约束条件的粒子
PEN=10^15; 
lam=PEN; lameq=PEN;
Z=0;
% 非线性约束
[g,geq]=fnonlin(u);

%通过不等式约束建立罚函数
for k=1:length(g),
    Z=Z+lam*g(k)^2*getH(g(k));
end
% 等式条件约束
for k=1:length(geq),
   Z=Z+lameq*geq(k)^2*geteqH(geq(k));
end

代码实现

等式约束判断

% 满足约束为0 反之为1
function H=geteqH(g)
if g==0,
    H=0;
else
    H=1; 
end

不等式约束判断

% 满足约束为0 反之为1
function H=getH(g)
if g<=0, 
    H=0; 
else
    H=1; 
end

罚函数

function Z=getconstraints(fnonlin,u)
% 一个很大的数，用于惩罚没有满足约束条件的粒子
PEN=10^15;
lam=PEN; lameq=PEN;

Z=0;
% 非线性约束
[g,geq]=fnonlin(u);

%通过不等式约束建立罚函数
for k=1:length(g),
    Z=Z+ lam*g(k)^2*getH(g(k));
end
% 等式条件约束
for k=1:length(geq),
   Z=Z+lameq*geq(k)^2*geteqH(geq(k));
end

随机性衰减因子

用于控制粒子更新的位置的随机性

% 随机性衰减因子
function alpha=newPara(alpha,gamma);
alpha=alpha*gamma;

寻找全局最优个体

  for i=1:n,   
    fval=Fun(fhandle,fnonlin,best(i,:)); 
    % 更新最有个体
    if fval<=fbest, 
        gbest=best(i,:);
        fbest=fval;
    end
  end

约束条件的建立

function [g,geq]=constraint(x)
% g(1) = x1 +x2 -x3 的形式
% g 用于存储不等式约束
% geq 用于存储不等式约束

带有约束条件的求解

function [gbest,fbest]=pso_mincon(fhandle,fnonlin,Lb,Ub,para)
if nargin<=4,
    para=[20 150 0.95];
end
n=para(1);% 粒子种群大小
time=para(2); %时间步长，迭代次数
gamma=para(3); %gama参数
 scale=abs(Ub-Lb); %取值区间
% 验证约束条件是否合乎条件
if abs(length(Lb)-length(Ub))>0,
    disp('Constraints must have equal size');
    return
end

  alpha=0.2; % alpha=[0,1]粒子随机衰减因子
  beta=0.5;  % 收敛速度(0->1)=(slow->fast);

% 初始化粒子群
best=init_pso(n,Lb,Ub);

fbest=1.0e+100;
% 迭代开始
for t=1:time,     
   
%寻找全局最优个体
  for i=1:n,   
    fval=Fun(fhandle,fnonlin,best(i,:)); 
    % 更新最有个体
    if fval<=fbest, 
        gbest=best(i,:);
        fbest=fval;
    end
        
  end

% 随机性衰减因子
 alpha=newPara(alpha,gamma);

% 更新粒子位置 
  best=pso_move(best,gbest,alpha,beta,Lb,Ub);  
 
% 结果显示
	str=strcat('Best estimates: gbest=',num2str(gbest));
	str=strcat(str,'  iteration='); str=strcat(str,num2str(t));
	disp(str);

    fitness1(t)=fbest;
    plot(fitness1,'r','Linewidth',2)
    grid on
    hold on
    title('适应度')
end

运行例子结果如下：