matlab PSO粒子群算法优化PID参数

Posted 2022-12-16 studyer_domi

1、内容简介

略
633可以交流、咨询、答疑

2、内容说明

粒子群算法是一种智能优化算法。关于智能，个人理解，不过是在枚举法的基础上加上了一定的寻优机制。试想一下枚举法，假设问题的解空间很小，比如一个函数 y = x^2 ，解空间在[-1,1],现在求这个函数的最小值，我们完全可以使用枚举法，比如在这里，在解空间[-1,1]上，取1000等分，也就是步长为0.002，生成1000个x值，然后代入函数中，找到这1000个最小的y就可以了。然而实际情况不是这样的，比如为什么选1000等分，不是1w，10w等分，很显然等分的越大，计算量也就越大，带来的解当然也就越精确，那么实际问题中如何去平衡这两点呢？也就是既要计算量小（速度快），也要准确（精度高），这就是智能算法的来源了，一般的智能算法基本上都是这样的，在很大的搜索空间上，即保证了速度快，也能比较好的找到最优解。

再来看看粒子群算法（也称PSO算法），也是一种进化算法，模拟生物群体的觅食行为，是一种群体智能算法，类似的算法想遗传算法，模拟退火算法等等。PSO是通过当前已知种群寻找到的所有解来决定新的解的寻找方向，也就是新解的生成方式依赖于这些种群历史上寻找的所有解。

形象的理解比如下图： 

开始随机生成一堆种群，那么这些种群之间的每个个体可以相互交流，比如下一时刻，A告诉B说我的解比你好，那么B就往A那个地方飞，也就是B的解朝着A的解方向变化，当然所有粒子间都这样操作，想想一旦粒子群中间有一个粒子找到了一个最优解，是不是所有的粒子会一窝蜂朝着这个方向而去了，而在这个去的过程中，万一某个粒子找到了一个更好的解，那它还会走吗？不会了，它就告诉剩下的所有粒子说我的解更好呀，大家快来呀（很无私的），然后所有粒子又一窝蜂的照着这个粒子方向前进，当然在这个前进的过程中可能又会产生新的解，就这样一步步的迭代，最终慢慢的趋近于一个最优解，这个解是不是全局最优解，不知道，可能是，也可能不是，取决于原始问题的复杂程度，也取决于粒子前进的多少等等。

粒子群算法相对于其他算法来说还是有很多优点的，典型的就是计算速度很快，在每次迭代时，所有粒子同时迭代，是一种并行计算方式，而且粒子的更新方式简单，朝着一个优秀解方向更新。这个优秀解包括两个部分： 
1）一个是朝着自己在迭代的历史上找到的个体最优解gbest前进 
2）一个是朝着群体在得带历史上找到的全体最优解zbest前进 
现在还有一个问题就是每次迭代的时候更新多少呢？也就是自变量的增加步长了，我们用一个速度量V来表示，也就是每个粒子的更新速度了，公式化的表示就是这样的： 

Vid(t+1)=Vidt+c1∗rand∗(Pid−xid(t))+c2∗rand∗(Pgd−xid(t))Vid(t+1)=Vidt+c1∗rand∗(Pid−xid(t))+c2∗rand∗(Pgd−xid(t))

3、仿真分析

略

%% 清空环境
clear
clc

%% 参数设置
wmax = 1;     % 最大惯性因子 
wmin = 0.6;   % 最小惯性因子 

% w = 0.6;      % 惯性因子 
c1 = 2;       % 加速常数
c2 = 2;       % 加速常数

Dim = 3;            % 维数
SwarmSize = 6;    % 粒子群规模
ObjFun = @PSO_PID;  % 待优化函数句柄

MaxIter = 4;      % 最大迭代次数  
MinFit = 0.1;       % 最小适应值 

Vmax = 1;
Vmin = -1;
Ub = [300 300 300];
Lb = [0 0 0];

%% 粒子群初始化
Range = ones(SwarmSize,1)*(Ub-Lb);
Swarm = rand(SwarmSize,Dim).*Range + ones(SwarmSize,1)*Lb;      % 初始化粒子群
VStep = rand(SwarmSize,Dim)*(Vmax-Vmin) + Vmin;                 % 初始化速度
fSwarm = zeros(SwarmSize,1);
for i=1:SwarmSize
    fSwarm(i,:) = PSO_PID(Swarm(i,:));                         % 粒子群的适应值
end

%% 个体极值和群体极值
[bestf,bestindex]=min(fSwarm);
zbest=Swarm(bestindex,:);   % 全局最佳
gbest=Swarm;                % 个体最佳
fgbest=fSwarm;              % 个体最佳适应值
fzbest=bestf;               % 全局最佳适应值

%% 迭代寻优
iter = 0;
y_fitness = zeros(1,MaxIter);   % 预先产生4个空矩阵
K_p = zeros(1,MaxIter);         
K_i = zeros(1,MaxIter);
K_d = zeros(1,MaxIter);
while( (iter < MaxIter) && (fzbest > MinFit) )
    w = wmax-(wmax-wmin)/MaxIter*iter;  %% 惯性权重因子调整
    for j=1:SwarmSize
        % 速度更新
        VStep(j,:) = w*VStep(j,:) + c1*rand*(gbest(j,:) - Swarm(j,:)) + c2*rand*(zbest - Swarm(j,:));
        if VStep(j,:)>Vmax, VStep(j,:)=Vmax; end
        if VStep(j,:)<Vmin, VStep(j,:)=Vmin; end
        % 位置更新
        Swarm(j,:)=Swarm(j,:)+VStep(j,:);
        for k=1:Dim
            if Swarm(j,k)>Ub(k), Swarm(j,k)=Ub(k); end
            if Swarm(j,k)<Lb(k), Swarm(j,k)=Lb(k); end
        end
        % 适应值
        fSwarm(j,:) = PSO_PID(Swarm(i,:));
        % 个体最优更新     
        if fSwarm(j) < fgbest(j)
            gbest(j,:) = Swarm(j,:);
            fgbest(j) = fSwarm(j);
        end
        % 群体最优更新
        if fSwarm(j) < fzbest
            zbest = Swarm(j,:);
            fzbest = fSwarm(j);
        end
    end 
    iter = iter+1;                      % 迭代次数更新
    y_fitness(1,iter) = fzbest;         % 为绘图做准备
    K_p(1,iter) = zbest(1);
    K_i(1,iter) = zbest(2);
    K_d(1,iter) = zbest(3);
end
%% 绘图输出
figure     % 绘制性能指标ITAE的变化曲线
plot(y_fitness,'LineWidth',2)
title('最优个体适应值','fontsize',18);
xlabel('迭代次数','fontsize',18);ylabel('适应值','fontsize',18);
set(gca,'Fontsize',18);

figure      % 绘制PID控制器参数变化曲线
plot(K_p)
hold on
plot(K_i,'k','LineWidth',3)
plot(K_d,'--r')
title('Kp、Ki、Kd 优化曲线','fontsize',18);
xlabel('迭代次数','fontsize',18);ylabel('参数值','fontsize',18);
set(gca,'Fontsize',18);
legend('Kp','Ki','Kd');

4、参考论文

略