基于粒子群算法的电机控制器参数整定研究

文章目录

1. 1. 智能优化算法概述
2. 2. 粒子群算法
   1. 2.1. 粒子群算法介绍
   2. 2.2. 粒子群算法原理
   3. 2.3. 改进粒子群算法
3. 3. PMSM仿真模型搭建
4. 4. 参数整定程序
   1. 4.1. 入口程序(main.m)
   2. 4.2. 粒子群算法函数(PSO.m)
   3. 4.3. 电机控制效果计算函数(cal_n.m)
   4. 4.4. 评价函数计算函数(PID_n_access.m)
   5. 4.5. 作图函数(draw.m)
5. 5. 仿真结果分析
6. 6. 参考资料

本科毕业设计课题，本文将对之进行重新修订整理。
本文主要介绍粒子群算法原理及参数整定原理，搭建双闭环PI控制的PMSM仿真模型，并建立粒子群算法电机控制器参数整定程序。将使用粒子群算法整定与传统工程整定方法得到的PI控制器参数分别进行仿真分析，对比二者的动态与稳态性能，来证明采用粒子群算法可以有效改善系统性能。

智能优化算法概述

智能优化算法是为解决信号解析、图像辨析、模式识别、任务调度、生产规划分配、控制系统设计等各行各业的复杂优化类问题，受到人类智能、生物种群的社会行为或者是自然界的固有现象而启发，通过现代化的计算机技术手段所设计出的算法。
常见的智能优化算法可以分为以下6类：传统智能优化算法（模糊PID控制算法、专家经验PID控制算法）、进化类智能算法（遗传算法、差分进化算法、免疫算法）、群体智能类算法（粒子群算法、蚁群算法、人群搜索算法、人工蜂群算法等）、模拟退火算法、禁忌搜索算法与神经网络算法。
在电机控制领域，引入了智能优化算法的电机 PID 控制器在参数整定上相较于传统参数整定方式具有更好的整定效果，可以有效提高电机的控制精度，更好地满足电机控制指标。

粒子群算法

粒子群算法介绍

粒子群优化算法(Particle Swarm Optimization, PSO)是由肯尼迪和埃伯哈特提出的一种新型智能优化算法，属于演化类智能优化算法，其算法灵感来源于自然界中鸟群觅食时的群体行为特点。粒子群算法将鸟群中的一只只觅食的鸟抽象为一个个粒子，这些粒子只具备坐标与速率这两个属性，鸟群中的鸟不知道食物的具体位置，但是知道其余各个鸟的具体位置，以及它们各自离食物的距离。也就是粒子群中各个粒子存在的信息共享机制，虽然单独的鸟未必每次运动都会离目标越近，但鸟群整体会越来越靠近目标食物位置。
粒子群算法易于理解，快速高效，在实时性强的复杂函数优化问题的最优化求解上效果显著，目前广泛应用于模式识别、神经深度学习、控制系统参数整定、目标函数最优化等各个场景。其缺点是运算过程中容易陷入局部最优解，导致所求解的不全面，使问题处理效果不尽人意。

粒子群算法原理

假设粒子群参数维度为\(N\)，粒子群中粒子总数\(m\)，迭代次数为\(n\)，设定迭代速度系数为\(c_1\)、\(c_2\)、\(\omega\)。
为确保每个粒子在限定范围内呈随机分散分布，输入向量\(x=(x_1,x_2,\cdots,x_N )^T\)的限制范围为： \[\left\{ \begin{gathered} x_{\min }^{(1)} < {x_1} < x_{\max }^{(1)} \hfill \\ x_{\min }^{(2)} < {x_2} < x_{\max }^{(2)} \hfill \\ ... \hfill \\ x_{\min }^{(N)} < {x_N} < x_{\max }^{(N)} \hfill \\ \end{gathered} \right.\] 记\(X_{min}=(x_{min}^{(1)},x_{min}^{(2)},\cdots,x_{min}^{(N)} )\)，\(X_{max}=(x_{max}^{(1)},x_{max}^{(2)},\cdots,x_{max}^{(N)})\)
输入向量应满足\(X_{min}<x<X_{max}\)
这里设定每个粒子的初始速度\(v_0 = 0\)，则第\(j\)个粒子\((1 \leqslant j \leqslant m)\)下一步的迭代速度可以用公式表示成：
\[{v^{(j)}} = \omega \cdot {v_0} + {c_1} \cdot rand \cdot ({P^{(j)}} - {X^{(j)}}) + {c_2} \cdot rand \cdot ({P_G} - {X^{(j)}})\]

迭代速度公式由三个部分相加而得，其中\(rand\)表示调用函数生成的0到1之间的随机数。

公式第一部分体现了粒子的自身因素，其包含粒子上一次迭代时的速度信息，\(v_0\)是上一次迭代时的粒子速度；
第二部分体现了粒子的学习因素，\(P^{(j)}\)为第\(j\)个粒子在历次迭代运动过程中最接近于待优化问题结果的那一次迭代时的位置（鸟离食物最近时的位置）；
第三部分体现粒子的交流因素，\(P_G\)表示当前的最接近于待优化问题结果的粒子的所在位置（鸟群离食物最近时的位置）。
粒子迭代速度系数可以根据仿真实际需要加以修改，结合实际情况与工程实践经验，可以取\(\omega=1,c_1=c_2=2\)，对于粒子群中粒子数量较多的情况，为确保运算结果能得到较快收敛，可以适当减小\(\omega\)的值。

根据迭代速度公式计算下一时刻的速度及其粒子位置，其公式可表达为：
\[X_{k + 1}^{(j)} = X_k^{(j)} + v_{k + 1}^{(j)} \cdot dt\] 式中\(X_{k+1}^{(j)}\)表示第\((k+1)\)步迭代的位置，\(X_k^{(j)}\)表示第\(k\)步迭代的位置，\(v_{k+1}^{(j)}\)表示第\((k+1)\)步迭代的速度，\(dt\)是仿真间隔，默认值为1，适当减小该值可以使仿真更慢更平滑。

计算\(P_G\)与\(P^{(j)}\)的公式如下所示，用于迭代过程中的最优粒子与粒子群的更新： \[{P_G} = \left\{ \begin{gathered} P_G^{(k)},F(P_G^{(k)}) > F(P_G^{(k + 1)}) \hfill \\ P_G^{(k + 1)},F(P_G^{(k)}) < F(P_G^{(k + 1)}) \hfill \\ \end{gathered} \right.\] \[P_G^{(j)} = \left\{ \begin{gathered} P_k^{(j)},F(P_k^{(j)}) > F(P_{k + 1}^{(j)}) \hfill \\ P_{k + 1}^{(j)},F(P_k^{(j)}) < F(P_{k + 1}^{(j)}) \hfill \\ \end{gathered} \right.\]

改进粒子群算法

在粒子群算法的参数调整中，惯性因子\(\omega\)以及学习因子\(c_1,c_2\)是关键。
学习因子主要影响粒子的目标方向感，即目标识别能力。
在进化上，\(c_1\)相对越大，粒子行为越独立，“分散”现象越明显，即学习因素比重较大，全局最优解颇具盲目性，算法收敛速度慢;
\(c_2\)相对越大，粒子群体意识越强，“群聚”现象越明显，交流因素比重大，但容易陷入局部最优解。

为解决传统的粒子群算法存在的收敛速度慢、容易陷入局部最优解等问题。需要对算法的惯性权重与学习因子在迭代过程中动态调整。
对于学习因子，在前期应当有较大的搜索空间，随着粒子的进化逐渐缩小搜索区域，向全局最优解加速进化，这就 要求\(c_1\)应当先大后小，\(c_2\)先小后大。
学习因子在迭代过程中的公式如下所示，这里引入了反正切函数进行优化。设定\(c_1\)与\(c_2\)的值域为\([0.1,2]\)。\(n\)为迭代总次数，\(k\)为当前迭代次数。
\(c_1=-\frac{4.265}{2\pi}arctan[\frac{8}{n}(k-\frac{n}{2})]+1\)
\(c_2=\frac{4.265}{2\pi}arctan[\frac{8}{n}(k-\frac{n}{2})]+1\)
注：arctan的结果单位是角度，式中\(2\pi\)为360°

对于惯性因子，应当先大后小，使算法在迭代初期增强全局搜索能力，更大范围遍历整体空间，避免陷入局部最优解；而在后期增强局部搜索能力，在更小范围内搜索到精确解。
惯性因子\(\omega\)在迭代过程中的公式如下所示，其在迭代过程中由权重最大值\(\omega_{max}\)线性衰减至权重最小值\(\omega_{min}\)。\(n\)为迭代总次数，\(k\)为当前迭代次数。
\(\omega=\omega_{max}-\frac{k}{n}(\omega_{max}-\omega_{min})\)
这里取\(\omega_{max}=0.9\)、\(\omega_{min}=0.4\)

PMSM仿真模型搭建

PMSM的数学模型参见鲁棒控制理论与应用02-H∞控制器的设计方法，本文沿用了上述数学模型，这里不再赘述。

电机具体参数如下：
定子电阻\({R_s} = 0.958\Omega\)，极对数\(p_n=4\)，粘滞摩擦系数\(B = 0.008N \cdot M \cdot s\)，d轴电感\({L_d} = 5.25mH\)，q轴电感\({L_q} = 12mH\)，永磁体的励磁磁链\({\Psi _f} = 0.1827wb\)，转动惯量\(J = 0.003kg \cdot {m^2}\)。

首先搭建双闭环PI控制的PMSM控制系统，使用工程整定方法确定转速环与电流环的PI控制器参数：
取转速环PI控制器参数\({k_p} = 0.2\)，\({k_i} = 30\)，输出\(i_q^*\)的限幅范围为-25.7A≤\(i_q^*\)≤25.7A；
取电流环\(i_q\)的PI控制器参数\(K_{pn} = 300\)，\(K_{In} = 23950\)，输出\(u_q^*\)的限幅范围为-161.6V≤\(u_q^*\)≤161.6V；
取电流环\(i_d\)的PI控制器参数\(K_{pi} = 131.25\)，\(K_{Ii} = 23950\)，输出\(u_d^*\)的限幅范围为-161.6V≤\(u_d^*\)≤161.6V。

在双闭环PI控制的PMSM控制基础上，引入改进的粒子群算法，对转速环的PI参数进行优化，限幅值及电流环PI参数保持不变。其控制系统原理图如下图所示： 图1 PMSM控制系统原理图

该系统采用空间电压矢量调制技术(Space Vector Pulse Width Modulation，SVPWM)。空间电压矢量调制又被称作称“磁链跟踪控制”。它以三相正弦波电压供电时交流电机的理想磁通轨迹为基准，通过三相功率逆变器的六个功率开关元件组成的特定开关模式来产生的脉宽调制波，用逆变器不同的开关模式产生的实际磁通去逼近基准磁通圆，从而达到较高的控制性能。

该系统采用\(i_d=0\)的矢量控制策略。该方法实现了对PMSM的解耦控制，具有控制简单、调速范围宽、转矩特性好、有效降低铜耗的优点。

图中，\(\omega^*\)为目标转速，\(\omega\)为实际转速；\(i_a,i_b,i_c\)为三相静止坐标系下的相电流；\(i_\alpha,i_\beta\)为经过CLARKE变换后的两相静止绕组电流；\(i_d,i_q\)为经过PARK变换后的两相旋转绕组电流；\(u_\alpha,u_\beta\)为\(u_d,u_q\)经过PARK逆变换得到的参数。

根据PMSM控制系统原理图，可以搭建出如下图所示的控制系统仿真模型。 图2 PMSM控制系统仿真模型

仿真总时长为0.2s，采样周期为\(10^{-5}s\)。设定电机转速给定值\(\omega_r^* = 1000r/\min\)，使其空载起动，0.1秒后对电机施加12\(N \cdot m\)的负载转矩。

参数整定程序

粒子群算法参数整定程序共由5个MATLAB脚本文件和1个Simulink仿真文件构成，运行入口程序后，将依次执行粒子群算法参数初始化，粒子群算法参数整定，Simulink仿真文件调用及仿真结果数据处理，仿真结果与参数整定结果图像输出等程序功能。

入口程序(main.m)

入口程序是第一个被执行的程序，用于初始化粒子群算法相关参数，包括迭代次数、粒子群规模、粒子限幅值等参数。并导入需要使用的其他MATLAB文件。

clear
global Kpn KIn bar
Kpn=0.2;KIn=30;
bar=0;
n=20;m=50;xmax=[2,500];xmin=[0.05,5];
cal_f=@cal_n;
f=@PID_n_access;
[Pg,fmin]=PSO(n,m,xmax,xmin,cal_f,f);
draw();
disp('OK');

粒子群算法函数(PSO.m)

粒子群算法函数的输入参数为迭代次数、粒子群的粒子总数、粒子限幅值、电机控制效果计算函数和评价函数计算函数，输出参数为评价函数最小值及其对应的电机PI控制参数。

function [Pg,fmin]=PSO(n,m,xmax,xmin,cal_f,f)
%全局粒子群算法,f为目标函数,dimension为维度,n为代数,m为种群规模
    dimension=numel(xmax);  %numel函数，返回数组元素个数   dimension=2
    wmax=0.9;wmin=0.4;%速度惯性系数 先大后小
    dt=1;%位移仿真间隔
    vmax=0.1*(xmax-xmin);%速度限幅
    v=zeros(m,dimension);%初始速度为0
    X=(xmax-xmin).*rand(m,dimension)+xmin;%初始位置满足(xmin,xmax)内均匀分布
    P=X;%P为每个粒子每代的最优位置
    last_f=f(X);
    [fmin,min_i]=min(last_f);%Pg为所有代中的最优位置 
    Pg=X(min_i,:);
    last_Pg=Pg;
    Pg_draw=zeros(n,dimension);%待绘图的PI参数集合
    legend_i=1;
    Pg_draw(legend_i,:)=Pg;
    [ts_n,sigma_n,ts_nb,sigma_nb,n_err]=cal_f(Pg);%计算初始粒子群的性能参数
    disp(['初始粒子群最优粒子 f=',num2str(fmin)]);
    disp(['sigma_n=',num2str(sigma_n),' ts_n=',num2str(ts_n),' sigma_nb=',num2str(sigma_nb),' ts_nb=',num2str(ts_nb),' n_err=',num2str(n_err)]);
    disp(['Kpn=',num2str(Pg(1)),' KIn=',num2str(Pg(2))]);
    for i=1:n
        w=wmax-i*(wmax-wmin)/n;%改进w 先大后小
        c1=-4.265/2/pi*atan(8*(i-n/2)/n)+1;%改进c1 先大后小
        c2=4.265/2/pi*atan(8*(i-n/2)/n)+1;%改进c2 先小后大
        v=w*v+c1*rand*(P-X)+c2*rand*(ones(m,1)*Pg-X);
        v=(v>vmax).*vmax+(v>=-vmax & v<=vmax).*v+(v<-vmax).*(-vmax);
        X=X+v*dt;
        X=(X>xmax).*xmax+(X>=xmin & X<=xmax).*X+(X<xmin).*(xmin);
        new_f=f(X);%新的目标函数值
        update_j=find(new_f<last_f);
        P(update_j,:)=X(update_j,:);%修正每个粒子的历史最优值
        [new_fmin,min_i]=min(new_f);
        new_Pg=X(min_i,:);
        Pg=(new_fmin<fmin)*new_Pg+(new_fmin>=fmin)*Pg;
        last_f=new_f;%保存当前的函数值
        fmin=min(new_fmin,fmin);%更新函数最小值
        [ts_n,sigma_n,ts_nb,sigma_nb,n_err]=cal_f(Pg);
        if last_Pg~=Pg
            legend_i=legend_i+1;
            Pg_draw(legend_i,:)=Pg;
        end
        last_Pg=Pg;
        disp(['第',num2str(i),'次迭代完成 f=',num2str(fmin)]);
        disp(['sigma_n=',num2str(sigma_n),' ts_n=',num2str(ts_n),' sigma_nb=',num2str(sigma_nb),' ts_nb=',num2str(ts_nb),' n_err=',num2str(n_err)]);
        disp(['Kpn=',num2str(Pg(1)),' KIn=',num2str(Pg(2))]);
    end
end

电机控制效果计算函数(cal_n.m)

电机控制效果计算函数的输入参数为一组（单个粒子）的电机转速环PI控制参数(\(K_{pn},K_{In}\))，首先将上述参数导入到双闭环PMSM控制系统仿真模型文件，然后运行仿真程序，提取所得转速数据，分析所得数据，求得转速环调节时间、转速环超调量、突加负载最大转速降落、突加负载转速环恢复时间与系统稳定后的转速波动，作为输出参数用于评价函数计算。

这里取系统从响应至保持在终值± 2%所经过的时间作为调节时间（或恢复时间），系统稳定后的转速波动为0.18s~0.2s时刻转速最大值与最小值的差值。

function [ts_n,sigma_n,ts_nb,sigma_nb,n_err]=cal_n(parameters)
    global Kpn KIn
    %将simulink所需参数导入工作区，以便仿真使用
    obj_slx='PMSM_SVPWM_2015b_Jiazai.mdl';
    Kpn=parameters(1);
    KIn=parameters(2);
    sim(obj_slx);
    n_data=n.signals.values;
    n_max=max(n_data);
    sigma_n=n_max-1000; %转速超调量
    ts_n1=find(n_data(1:10000)>1020);
    ts_n2=find(n_data(1:10000)<980);
    if(numel(ts_n1)==0)
        ts_n=ts_n2(end);
    else
        ts_n=ts_n1(end);
    end
    n_min=min(n_data(10001:20000));
    sigma_nb=1000-n_min; %突加负载最大转速降落
    ts_nb1=find(n_data(10001:20000)<980);
    if(numel(ts_nb1)==0)
        ts_nb=0;
    else
        ts_nb=ts_nb1(end);
    end
    n_infty_max=max(n_data(18000:20000));%转速终值
    n_infty_min=min(n_data(18000:20000));%转速终值
    n_err=n_infty_max-n_infty_min;%转速波动
end

评价函数计算函数(PID_n_access.m)

评价函数是粒子群算法待优化的函数，本程序的输入参数是待处理的所有粒子群参数，程序首先将将粒子群逐一分散为粒子，针对每个粒子所对应的电机PI控制参数逐一单独仿真，带入参数计算程序，计算所得的各个控制效果参数，包括转速环调节时间\(T_{sn}\)、转速环超调量\(\omega_{sn}\)、突加负载最大转速降落\(n_d\)、突加负载转速环恢复时间\(T_{sd}\)、系统稳定后的转速波动\(n_{err}\)。
最后将上述参数带入如下式所示的评价函数，其中\(T_{sn0},\omega_{sn0},n_{d0},T_{sd0}\)为基准值。
\[f=3*\log(\frac{T_{sn}}{T_{sn0}}+1)+3*\log(\frac{\omega_{sn}}{\omega_{sn0}}+1)+\log(\frac{n_d}{n_{d0}}+1)+\log(\frac{T_{sd}}{T_{sd0}}+1)+10*n_{err}\] 上述五个参数都是越小越好，这使得所得函数值越小，代表电机控制效果越好。各项式子前的系数代表着各个参数对于评价函数值的影响权重，可以根据实际控制性能的不同需求修改这些系数。

function [y,ts_n,sigma_n,ts_nb,sigma_nb,n_err]=PID_n_access(parameters_list)
    global bar
    M=size(parameters_list,1);%计算需要仿真的组数（simulink无法并行仿真）
    ts_n0=1200.0;sigma_n0=10.0;ts_nb0=600;sigma_nb0=20.0;%基准值
    %将simulink所需参数导入工作区，以便仿真使用
    %if((mod(bar,4)==1&&bar<=200)||(mod(bar,2)==1&&bar>=200))
    y=zeros(M,1);
    for sim_i=1:M
        bar=bar+1;
        if(mod(bar,2)==1)
            disp('█');
            fprintf('%c',8);
        end
        parameters=parameters_list(sim_i,:);
        [ts_n,sigma_n,ts_nb,sigma_nb,n_err]=cal_n(parameters);
        y(sim_i)=3*log(ts_n/ts_n0+1)+3*log(sigma_n/sigma_n0+1)+log(ts_nb/ts_nb0+1)+log(sigma_nb/sigma_nb0+1)+n_err*10;
    end
end

作图函数(draw.m)

用于绘制参数整定前后的转速曲线图。

function draw()
    global Kpn KIn
    figure(1)
    sim('PMSM_SVPWM_2015b_Jiazai.mdl')
    n_data=n.signals.values;time=n.time;%获取转速序列及对应时刻
    plot(time,n_data)
    hold on
    
    figure(1)
    Kpn=0.2;
    KIn=30;
    sim('PMSM_SVPWM_2015b_Jiazai.mdl')
    n_data=n.signals.values;time=n.time;%获取转速序列及对应时刻
    plot(time,n_data)
    hold on
    
    figure(1)
    legend_str{1}='PSO-PI';
    legend_str{2}='PI';
    legend(legend_str)
    xlabel('时间/s') 
    ylabel('电机转速n/(r/min)') 
    disp('比较曲线绘制完成');

仿真结果分析

图3 粒子群算法参数整定前后转速曲线比较

运行代码后，输出参数整定前后的转速曲线如上图所示。其具体参数见下表：

表1 粒子群算法参数前后的控制性能比较

	PI	PSO-PI
转速环PI参数	\(K_{pn}=0.2,K_{In}=30\)	\(K_{pn}=0.24005,K_{In}=276.6148\)
转速环调节时间/ms	21.57	14.54
转速环超调量/(r/min)	44.167	28.0042
突加负载最大转速降落/(r/min)	40.1785	24.7302
转速环恢复时间/ms	8.48	1.99
稳态转速波动/(r/min)	0.17	0.09

由上表可知，采用粒子群算法整定后的PID参数，无论是在电机启动速度及超调量，还是加入负载转矩后恢复时的属性，都显著优于工程整定方法所得的参数。

参考资料

• CSDN—沉沙丶—永磁同步电机矢量控制系列教程
• 【Simulink】粒子群算法（PSO）整定PID参数（附代码和讲解）
• 刘子健. 基于改进的PSO-BP神经网络算法的永磁同步电机PID参数整定[J]. 兰州文理学院学报（自然科学版），2024，38(1)：77-81.
• 谢金法. 基于改进粒子群优化算法的PMSM矢量控制[J]. 现代制造工程，2019(7)：6-11.
• 曾雄飞. 基于粒子群算法优化BP神经网络的PID控制算法[J]. 电子设计工程，2022，30(11)：69-73.
• 金宝宝. 基于智能算法的无刷直流电机控制系统设计与实现[D]. 兰州交通大学，2023