粒子群算法优化PID参数实例

Posted 2022-03-14 studyer_domi

%% 清空环境
clear all
clc
tic; %记录算法运行时间
%% 参数设置
%标准PSO算法
 w = 0.6;      % 惯性因子 
 c1 = 2;       % 加速常数
 c2 = 2;       % 加速常数

%%经验值
% w = 0.7298;      % 惯性因子 
% c1 = 1.49618;       % 加速常数
% c2 = 1.49618;       % 加速常数

%%随机权重策略
% mean_max = 0.8;  %随机权重均值最大值
% mean_min = 0.5;  %随机权重均值最小值
% sigma = 0.2;          %方差

%以下参数值都可以改

Dim = 3;            % 维数
SwarmSize = 30;    % 粒子群规模
ObjFun = @PSO_PID;  % 待优化函数句柄，来自于PSO_PID.m

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

Vmax = 1;
Vmin = -1;
% Ub = [300 300 300];  %变量上限
 Ub = [40 1 40];
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,:) = feval(ObjFun,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) )
    for j=1:SwarmSize
%         miu = mean_min+(mean_max-mean_min)*rand;
%         w = miu +sigma*randn();
        
        % 速度更新
        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,:) = feval(ObjFun,Swarm(j,:));
        % 个体最优更新     
        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

elapsed_time=toc % time taken for this algorithm

%% 绘图输出
figure(1)      % 绘制性能指标ITAE的变化曲线
plot(y_fitness,'LineWidth',2)
title('ITAE fitness curve','fontsize',12);
xlabel('No of Iteration','fontsize',12);ylabel('value of fitness','fontsize',12);
set(gca,'Fontsize',12);
%下面代码读取当前电脑时间命名并自动figure(1)保存为.fig格式，在当前文件位置
h1=figure(1);
fn1=['ITAE by PSO-',datestr(clock,'yy-mm-dd-HH-MM-SS'),'.fig'];
saveas(h1,fn1);


%下面图2是Kp Ki Kd的值变化过程，可以不画不保存
figure(2)      % 绘制PID控制器参数变化曲线
plot(K_p);
hold on
plot(K_i,'k','LineWidth',3);
plot(K_d,'--r');
title('Kp、Ki、Kd','fontsize',12);
xlabel('No of Iteration','fontsize',12);
ylabel('value of KpKiKd','fontsize',12);
set(gca,'Fontsize',12);
legend('Kp','Ki','Kd');
%下面代码读取当前电脑时间命名并自动figure(2)保存为.fig格式，在当前文件位置
h2=gcf;
fn2=['Kpid by PSO-',datestr(clock,'yy-mm-dd-HH-MM-SS'),'.fig'];
saveas(h2,fn2);