Matlab算法粒子群算法求解二维线性优化问题（附MATLAB代码）

Posted 2023-04-02 Albert_Lsk

MATLAB求解二维线性优化问题

• 前言
• 正文
• 函数实现
• 可视化结果

前言

二维线性优化问题指的是在二维空间中，对于一个由线性函数构成的目标函数，通过限制自变量的范围或满足特定的约束条件，寻找一个最优解（最小值或最大值）。这个问题可以用以下形式的数学模型来描述：

$$\begin{array}{c}\begin{array}{r}min/maxz=c1\ast x1+c2\ast x2\\ subjectto:a1\ast x1+a2\ast x2<=b1\\ a3\ast x1+a4\ast x2<=b2\\ lb<=x1<=ub\\ lb<=x2<=ub\end{array}\end{array}$$

其中，$x_1$和$x_2$是自变量，$c_1$和$c_2$是常数，$a_1$、$a_2$、$a_3$、$a_4$、$b_1$和$b_2$也是常数，表示约束条件和目标函数的系数。最小化$z$表示寻找目标函数的最小值，最大化$z$表示寻找目标函数的最大值。在约束条件下，$x_1$和$x_2$的取值必须满足上下界（$lb$和$ub$）的限制和不等式约束。在这个问题中，优化的目标是找到一组满足约束条件的自变量值，使得目标函数的值最小或最大。

粒子群优化（PSO）算法是一种基于群体智能的优化算法，它可以用于优化二维线性问题。以下是使用PSO算法解决二维线性优化问题的基本步骤：

1. 定义适应度函数：将目标函数定义为适应度函数，即 $f(x_1,x_2)=c_1{x}_1+c_2{x}_2$。
2. 初始化种群：设置粒子数量、位置范围、速度范围等参数，随机生成粒子的位置和速度。
3. 计算适应度：根据粒子的位置计算适应度值。
4. 更新全局最优解和局部最优解：根据适应度值更新全局最优解和局部最优解。
5. 更新粒子速度和位置：根据当前位置和速度、全局最优解和局部最优解更新粒子的速度和位置。
6. 检查终止条件：检查粒子的位置是否达到预设的最大迭代次数或适应度值是否达到预设的最小值，若满足则停止迭代，输出最优解。

在优化过程中，我们可以使用Matlab进行可视化操作，将优化过程和最终结果以图形化的方式展示出来，更直观地观察算法的运行情况和优化效果。

正文

下面我们使用粒子群优化算法求解下列二维线性函数的优化问题：
$$minf(x)=x{1}^2+x{2}^2-2x1-4x2+5$$

函数实现

% 设置PSO参数
num_particles =  100; % 粒子数量
num_iterations = 200; % 迭代次数
c1 = 1.5; % 学习因子1
c2 = 1.5; % 学习因子2
w = 0.7; % 惯性权重
lb = -10; % 自变量下限
ub =  10; % 自变量上限
n = 2; % 自变量个数

% 初始化粒子群
particles = rand(num_particles, n) .* (ub - lb) + lb;
velocities = zeros(num_particles, n);
pbest = particles;
pbest_fitness = inf(num_particles, 1);
gbest = zeros(1, n);
gbest_fitness = inf;

% 迭代PSO算法
for i = 1:num_iterations
    % 更新粒子速度和位置
    for j = 1:num_particles
        % 计算学习因子
        r1 = rand;
        r2 = rand;
        phi_p = c1 * r1;
        phi_g = c2 * r2;
        
        % 更新速度和位置
        velocities(j,:) = w * velocities(j,:) ...
            + phi_p * (pbest(j,:) - particles(j,:)) ...
            + phi_g * (gbest - particles(j,:));
        particles(j,:) = particles(j,:) + velocities(j,:);
        
        % 确保粒子位置在边界内
        particles(j,:) = max(particles(j,:), lb);
        particles(j,:) = min(particles(j,:), ub);
    end
    
    % 更新个体最优解和全局最优解
    for j = 1:num_particles
        fitness = particles(j,1)^2 + particles(j,2)^2 - 2*particles(j,1) - 4*particles(j,2) + 5;
        if fitness < pbest_fitness(j)
            pbest(j,:) = particles(j,:);
            pbest_fitness(j) = fitness;
        end
        
        if fitness < gbest_fitness
            gbest = particles(j,:);
            gbest_fitness = fitness;
        end
    end
    
    % 可视化当前迭代的粒子群和最优解
    scatter(particles(:,1), particles(:,2), 'k.');
    hold on;
    scatter(gbest(1), gbest(2), 'r*');
    hold off;
    title(sprintf('Iteration %d, Best Fitness: %g', i, gbest_fitness));
    drawnow;
end

在上述代码中，我们设置了粒子数量为$100$，迭代次数为$200$。学习因子$c1$和$c2$分别设置为$1.5$，惯性权重$w$设置为$0.7$。

可视化结果

通过上图可知求解出的$x_1,x_2$分别为$1$和$2$。

优化求解基于非线性动态自适应惯性权重粒子群算法(IPSO)Matlab源码

 一、粒子群算法

粒子群算法是在1995年由Eberhart博士和Kennedy博士一起提出的，它源于对鸟群捕食行为的研究。它的基本核心是利用群体中的个体对信息的共享从而使整个群体的运动在问题求解空间中产生从无序到有序的演化过程，从而获得问题的最优解。设想这么一个场景：一群鸟进行觅食，而远处有一片玉米地，所有的鸟都不知道玉米地到底在哪里，但是它们知道自己当前的位置距离玉米地有多远。那么找到玉米地的最佳策略，也是最简单有效的策略就是搜寻目前距离玉米地最近的鸟群的周围区域。

在PSO中，每个优化问题的解都是搜索空间中的一只鸟，称之为"粒子"，而问题的最优解就对应于鸟群中寻找的"玉米地"。所有的粒子都具有一个位置向量（粒子在解空间的位置）和速度向量（决定下次飞行的方向和速度），并可以根据目标函数来计算当前的所在位置的适应值（fitness value）,可以将其理解为距离"玉米地"的距离。在每次的迭代中，种群中的例子除了根据自身的经验（历史位置）进行学习以外，还可以根据种群中最优粒子的"经验"来学习，从而确定下一次迭代时需要如何调整和改变飞行的方向和速度。就这样逐步迭代，最终整个种群的例子就会逐步趋于最优解。

上面的解释可能还比较抽象，下面通过一个简单的例子来进行说明

在一个湖中有两个人他们之间可以通信，并且可以探测到自己所在位置的最低点。初始位置如上图所示，由于右边比较深，因此左边的人会往右边移动一下小船。

现在左边比较深，因此右边的人会往左边移动一下小船

一直重复该过程，最后两个小船会相遇

得到一个局部的最优解

将每个个体表示为粒子。每个个体在某一时刻的位置表示为，x(t),方向表示为v(t)

p（t）为在t时刻x个体的自己的最优解，g(t)为在t时刻所有个体的最优解，v(t)为个体在t时刻的方向，x(t)为个体在t时刻的位置

下一个位置为上图所示由x,p,g共同决定了

种群中的粒子通过不断地向自身和种群的历史信息进行学习，从而可以找到问题的最优解。

但是，在后续的研究中表表明，上述原始的公式中存在一个问题：公式中V的更新太具有随机性，从而使整个PSO算法的全局优化能力很强，但是局部搜索能力较差。而实际上，我们需要在算法迭代初期PSO有着较强的全局优化能力，而在算法的后期，整个种群应该具有更强的局部搜索能力。所以根据上述的弊端，shi和Eberhart通过引入惯性权重修改了公式，从而提出了PSO的惯性权重模型：

每一个向量的分量表示如下

其中w称为是PSO的惯性权重，它的取值介于【0,1】区间，一般应用中均采用自适应的取值方法，即一开始令w=0.9，使得PSO全局优化能力较强，随着迭代的深入，参数w进行递减，从而使的PSO具有较强的局部优化能力，当迭代结束时，w=0.1。参数c1和c2称为学习因子，一般设置为1,4961；而r1和r2为介于[0,1]之间的随机概率值。

整个粒子群优化算法的算法框架如下：

step1种群初始化，可以进行随机初始化或者根据被优化的问题设计特定的初始化方法，然后计算个体的适应值，从而选择出个体的局部最优位置向量和种群的全局最优位置向量。

step2 迭代设置：设置迭代次数，并令当前迭代次数为1

step3 速度更新：更新每个个体的速度向量

step4 位置更新：更新每个个体的位置向量

step5 局部位置和全局位置向量更新：更新每个个体的局部最优解和种群的全局最优解

step6 终止条件判断：判断迭代次数时都达到最大迭代次数，如果满足，输出全局最优解，否则继续进行迭代，跳转至step 3。

　　对于粒子群优化算法的运用，主要是对速度和位置向量迭代算子的设计。迭代算子是否有效将决定整个PSO算法性能的优劣，所以如何设计PSO的迭代算子是PSO算法应用的研究重点和难点。

二、竞争学习的粒子群优化算法(CLPSO)

      传统的时变惯性权重粒子群优化算法对于求解一般的全局最优问题具有良好的效果,而对于复杂高维的优化问题易陷入局部收敛、存在早熟等缺点。针对以上存在的缺点,提出了种群进化离散度的概念,并考虑Sigmoid函数在线性与非线性之间较好的平衡性能,给出一种非线性动态自适应惯性权重的粒子群优化算法。该算法充分考虑进化过程中种群粒子之间进化差异,自适应地赋予不同的惯性权重因子,满足粒子群优化算法在不同进化时期对全局探索和局部开发能力的需求,仿真实例测试结果验证了该算法的有效性。 

三、部分代码

function [gbest,gbestval,fitcount]= CLPSO_new_func(fhd,Max_Gen,Max_FES,Particle_Number,Dimension,VRmin,VRmax,varargin)
%[gbest,gbestval,fitcount]= CLPSO_new_func('f8',3500,200000,30,30,-5.12,5.12)
rand('state',sum(100*clock));
me=Max_Gen;
ps=Particle_Number;
D=Dimension;
cc=[1 1];   %acceleration constants
t=0:1/(ps-1):1;t=5.*t;
Pc=0.0+(0.5-0.0).*(exp(t)-exp(t(1)))./(exp(t(ps))-exp(t(1)));
% Pc=0.5.*ones(1,ps);
m=0.*ones(ps,1);
iwt=0.9-(1:me)*(0.7/me);
% iwt=0.729-(1:me)*(0.0/me);
cc=[1.49445 1.49445];
if length(VRmin)==1
    VRmin=repmat(VRmin,1,D);
    VRmax=repmat(VRmax,1,D);
end
mv=0.2*(VRmax-VRmin);
VRmin=repmat(VRmin,ps,1);
VRmax=repmat(VRmax,ps,1);
Vmin=repmat(-mv,ps,1);
Vmax=-Vmin;
pos=VRmin+(VRmax-VRmin).*rand(ps,D);

for i=1:ps;
    e(i,1)=feval(fhd,pos(i,:),varargin{:});
end

fitcount=ps;
vel=Vmin+2.*Vmax.*rand(ps,D);%initialize the velocity of the particles
pbest=pos;
pbestval=e; %initialize the pbest and the pbest's fitness value
[gbestval,gbestid]=min(pbestval);
gbest=pbest(gbestid,:);%initialize the gbest and the gbest's fitness value
gbestrep=repmat(gbest,ps,1);

stay_num=zeros(ps,1);

ai=zeros(ps,D);
f_pbest=1:ps;f_pbest=repmat(f_pbest',1,D);
for k=1:ps
    ar=randperm(D);
    ai(k,ar(1:m(k)))=1;
    fi1=ceil(ps*rand(1,D));
    fi2=ceil(ps*rand(1,D));
    fi=(pbestval(fi1)<pbestval(fi2))'.*fi1+(pbestval(fi1)>=pbestval(fi2))'.*fi2;
    bi=ceil(rand(1,D)-1+Pc(k));
    if bi==zeros(1,D),rc=randperm(D);bi(rc(1))=1;end
    f_pbest(k,:)=bi.*fi+(1-bi).*f_pbest(k,:);
end

stop_num=0;
i=1;


while i<=me&fitcount<=Max_FES
    i=i+1;
    for k=1:ps
        
        if stay_num(k)>=5
            %     if round(i/10)==i/10%|stay_num(k)>=5
            stay_num(k)=0;
            ai(k,:)=zeros(1,D);
            f_pbest(k,:)=k.*ones(1,D);
            ar=randperm(D);
            ai(k,ar(1:m(k)))=1;
            fi1=ceil(ps*rand(1,D));
            fi2=ceil(ps*rand(1,D));
            fi=(pbestval(fi1)<pbestval(fi2))'.*fi1+(pbestval(fi1)>=pbestval(fi2))'.*fi2;
            bi=ceil(rand(1,D)-1+Pc(k));
            if bi==zeros(1,D),rc=randperm(D);bi(rc(1))=1;end
            f_pbest(k,:)=bi.*fi+(1-bi).*f_pbest(k,:);
        end
        
        for dimcnt=1:D
            pbest_f(k,dimcnt)=pbest(f_pbest(k,dimcnt),dimcnt);
        end
        aa(k,:)=cc(1).*(1-ai(k,:)).*rand(1,D).*(pbest_f(k,:)-pos(k,:))+cc(2).*ai(k,:).*rand(1,D).*(gbestrep(k,:)-pos(k,:));%~~~~~~~~~~~~~~~~~~~~~~
        vel(k,:)=iwt(i).*vel(k,:)+aa(k,:);
        vel(k,:)=(vel(k,:)>mv).*mv+(vel(k,:)<=mv).*vel(k,:);
        vel(k,:)=(vel(k,:)<(-mv)).*(-mv)+(vel(k,:)>=(-mv)).*vel(k,:);
        pos(k,:)=pos(k,:)+vel(k,:);
        
        if (sum(pos(k,:)>VRmax(k,:))+sum(pos(k,:)<VRmin(k,:)))==0;
            e(k,1)=feval(fhd,pos(k,:),varargin{:});
            fitcount=fitcount+1;
            tmp=(pbestval(k)<=e(k));
            if tmp==1
                stay_num(k)=stay_num(k)+1;
            end
            temp=repmat(tmp,1,D);
            pbest(k,:)=temp.*pbest(k,:)+(1-temp).*pos(k,:);
            pbestval(k)=tmp.*pbestval(k)+(1-tmp).*e(k);%update the pbest
            if pbestval(k)<gbestval
                gbest=pbest(k,:);
                gbestval=pbestval(k);
                gbestrep=repmat(gbest,ps,1);%update the gbest
            end
        end
        
    end
    
    % if round(i/100)==i/100
    %     plot(pos(:,D-1),pos(:,D),'b*');hold on;
    %     for k=1:floor(D/2)
    %         plot(gbest(:,2*k-1),gbest(:,2*k),'r*');
    %     end
    %     hold off
    %     title(['PSO: ',num2str(i),' generations, Gbestval=',num2str(gbestval)]);
    %     axis([VRmin(1,D-1),VRmax(1,D-1),VRmin(1,D),VRmax(1,D)])
    %     drawnow
    % end
    
    if fitcount>=Max_FES
        break;
    end
    if (i==me)&(fitcount<Max_FES)
        i=i-1;
    end
end
gbestval

四、仿真结果 

 五、参考文献及代码私信博主

[1]王生亮,刘根友.一种非线性动态自适应惯性权重PSO算法[J].计算机仿真,2021,38(04):249-253+451.