精华内容
下载资源
问答
  • 基于大量的试验数据,分析结合面面压、粗糙度、材料、介质等因素对固定结合面法向动态特性参数的影响规律,并验证以往相关理论的正确性,提出改进固定结合面法向动态特性参数的措施,为机床的模态分析及优化设计提供了...
  • 在直接数字域设计中,我们常常需要用到PID算法,而PID算法投入单片机使用时,往往需要硬件的支持,在调试时非常麻烦。本文通过Matlab仿真的手段实现PID,方便了开发者对系统的设计和实时调试。

    0.符号说明

    1. y(k)——系统响应输出的离散值
    2. u(k)——数字PID控制输出的离散值
    3. r(k)——期望输出的离散值(事先已知),在本例中为常数(即阶跃输入)
    4. e(k)——e(k)=r(k)-y(k),为期望值-实际值,是单位负反馈的误差比较信号
      图片来源于百度百科
      注:图片来源于百度百科

    1.如何根据连续系统建立差分方程

    1.1.获取连续系统的传递函数

    线性定常系统的控制中,PID是个非常常见的控制方式,如果可以通过Matlab仿真出PID的控制效果图,那么对系统设计时的实时调试将会容易得多。在这里我们将会以一个利用系统辨识参数的PID设计为为例展示Matlab仿真PID的过程。
    首先需要对一个未知的系统的参数进行辨识,以延迟环节可以忽略不计的电机调速系统为例。将时间戳导入xdata向量,对应的时刻转速导入ydata向量,进行系统辨识

    链接:Matlab的系统辨识

    我们就以上文链接中辨识的系统传递函数为例:
    G ( s ) = 0.998 0.021 s + 1 G(s)=\frac{0.998}{0.021s+1} G(s)=0.021s+10.998因此通过tf函数建立系统结构体如下:

    sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数,即式4.1
    

    1.2.获取离散系统的传递函数

    由于是数字PID仿真,我们需要选取一个采样时间,本案例选用的是0.005s(注意,采样周期应该小于系统纯滞后时间的0.1倍)。在对其进行数字PID控制前,我们需要将这个系统离散化:

    ts=0.005;  %采样时间=0.005s
    dsys=c2d(sys,ts,'z');      %离散化
    

    dsys即我们根据采样周期离散化的Z变换系统。首先我们需要提取这个Z变化d那系统的参数方便后面的计算:

    [num,den]=tfdata(dsys,'v');%'v'代表强制以向量的格式(默认为元胞数组)输出num和den
    

    1.3.转换为差分方程

    求解出的Z变换表达式为 d s y s = n u m ( 1 ) ⋅ z + n u m ( 2 ) d e n ( 1 ) ⋅ z + d e n ( 2 ) = 0.2114 z − 0.7881 dsys=\frac{num(1)\cdot z +num(2)}{den(1)\cdot z+den(2)}=\frac{0.2114}{z-0.7881} dsys=den(1)z+den(2)num(1)z+num(2)=z0.78810.2114
    在PID仿真的过程中我们需要求解出时域表达式 ,因此需要借助差分方程解决,对于以下的Z变换:

    \begin{equation}
    Y(z)=dsys\cdot U(z)=\frac{num(2)}{den(1)\cdot z+den(2)}\cdot U(z)
    \label{eq:Sample1}
    \end{equation}

    \begin{equation}
    zY(z)+den(2)Y(z)=num(1)zU(z)+num(2)U(z)
    \label{eq:Sample2}
    \end{equation}
    对上式进行反Z变换,可以得到以下的差分方程:

    \begin{equation}
    y(k+1)+den(2)y(k)=num(1)u(k+1)+num(2)u(k)
    \label{eq:Sample3}
    \end{equation}

    \begin{equation}
    y(k+1)=-den(2)y(k)+num(1)u(k+1)+num(2)u(k)
    \label{eq:Sample4}
    \end{equation}
    位置型PID仿真时实际上可以不需要保存前一个数据(u(k)和y(k)),增量型PID必须要保存前一个数据。这里我们使用了位置型PID,但仍然利用 u 1 u_1 u1 y 1 y_1 y1保存了上一个数据,仅仅是为了演示这一过程。\begin{equation}
    y(k+1)=-den(2)y(k)+num(1)u(k+1)+num(2)u(k)
    \end{equation}
    可以转换为下面的式子:
    \begin{equation}
    y(k)=-den(2)y_1+num(1)u(k)+num(2)u_1
    \label{eq:Sample5}
    \end{equation}
    我们的差分方程就这样建立完毕。注意,此差分方程仅仅是描述系统模型的运算规律的,和我们的控制无关。因此是y(k)和u(k)的映射关系。我们下面的控制则是利用负反馈信号e(k)导出u(k)的输出,求解的是控制器u(k)的序列值。

    2.基本PID控制原理

    以位置型PID控制为例。将连续的PID控制转换为数字式时,微分环节被用差分代替,积分环节被累加和代替,比例环节则保持不变。差分的实现非常简单,只需要用 e ( k + 1 ) − e ( k ) e(k+1)-e(k) e(k+1)e(k) e ( k ) − e 1 e(k)-e_1 e(k)e1等效即可。积分的实现在每一次运算的后面都累加原来的误差,即Ee=Ee+e_1;即可。PID的控制器输出 u ( k ) = K p ⋅ e ( k ) + K d ⋅ ( e ( k ) − e 1 ) + K i ⋅ E e u(k)=Kp\cdot e(k)+Kd\cdot (e(k)-e_1)+Ki\cdot Ee u(k)=Kpe(k)+Kd(e(k)e1)+KiEe
    PID控制器构造完毕,我们需要通过r(k)和y(k)得到e(k),再通过e(k)得出u(k),进而再求解出y(k),再结合r(k)求解出e(k),…以此循环,求解出离散的响应点。
    详细的代码如下:

    ts=0.005;  %采样时间=0.005s
    sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数,即式4.1
    dsys=c2d(sys,ts,'z');      %离散化
    [num,den]=tfdata(dsys,'v');   %
    e_1=0;      %前一时刻的偏差      
    Ee=0;       %累积偏差
    u_1=0.0;    %前一时刻的控制量
    y_1=0;       %前一时刻的输出
    %PID参数
    kp=0.22;    
    ki=0.13;
    kd=0;
    u=zeros(1,1000);%预先分配内存
    time=zeros(1,1000);%时刻点(设定1000个)
    for k=1:1:1000
        time(k)=k*ts;   %时间参数
        r(k)=1500;      %期望值
        y(k)=-1*den(2)*y_1+num(2)*u_1+num(1)*u(k);%系统响应输出序列
        e(k)=r(k)-y(k);   %误差信号
        u(k)=kp*e(k)+ki*Ee+kd*(e(k)-e_1); %系统PID控制器输出序列
        Ee=Ee+e(k);    %误差的累加和
        u_1=u(k);    	%前一个的控制器输出值
        y_1=y(k);    	%前一个的系统响应输出值
        e_1=e(k);		%前一个误差信号的值
    end
    %(仅绘制过渡过程的曲线,x坐标限制为[0,1])
    p1=plot(time,r,'-.');xlim([0,1]);hold on;%指令信号的曲线(即期望输入)
    p2=plot(time,y,'--');xlim([0,1]);%不含积分分离的PID曲线
    hold on;
    

    输出的PID控制曲线如下:
    PID控制

    3.比较PID输出,分析参数产生的影响

    一个基本的PID就完成了。下面如果我们想要知道修改PID的三个参数kp,ki,kd会带来什么效果,只需要在程序中修改即可。为了方便起见,我们建立一个PID的数组,kp,ki,kd每次都取数组的一个值,然后设定一个大循环开始循环仿真。再利用subplot输出子图的方式将所有的PID效果都输出到一个图进行对比。该代码根据上述代码修改已经很容易,PID比较图的代码如下:

    close all
    PID=[0.22,0.13,0;
        0.4,0.13,0;
        0.4,0.25,0;
        0.8,0.23,0.4;
        0.8,0.2,1;
        0.7,0.2,0.9];%初始化PID参数
    for pid=1:1:6
    ts=0.005;  %采样时间=0.005s
    sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数,即式4.1
    dsys=c2d(sys,ts,'z');      %离散化
    [num,den]=tfdata(dsys,'v');   %
    e_1=0;      %前一时刻的偏差      
    Ee=0;       %累积偏差
    u_1=0.0;    %前一时刻的控制量
    y_1=0;       %前一时刻的输出
    %PID参数
    kp=PID(pid,1);    
    ki=PID(pid,2);
    kd=PID(pid,3);
    u=zeros(1,1000);
    time=zeros(1,1000);
    for k=1:1:1000
        time(k)=k*ts;   %时间参数
        r(k)=1500;      %给定量
        y(k)=-1*den(2)*y_1+num(2)*u_1+num(1)*u(k);
        e(k)=r(k)-y(k);   %偏差
        u(k)=kp*e(k)+ki*Ee+kd*(e(k)-e_1);   
        Ee=Ee+e(k);    
        u_1=u(k);    
        y_1=y(k);    
        e_1=e(k);
    end
    subplot(2,3,pid);
    p1=plot(time,r,'-.');xlim([0,1]);hold on;
    p2=plot(time,y,'--');xlim([0,1]);
    title(['Kp=',num2str(kp),' Ki=',num2str(ki),' Kd= ',num2str(kd)]);
    hold on;
    end
    

    输出的子图矩阵如下:
    PID子图矩阵
    可以发现,修改Kp会造成上升时间的缩短,但是有可能也会带来较大的超调。积分的增加是一个严重的滞后环节,会减小相位裕度,也会带来超调(超调量并不是绝对的,相对于较小的Kp可能会产生较大的超调,而Kp较大时超调会减小(例如第一行的1图和2图的对比))。然而积分的引入也是必要的,否则将会很长时间无法削弱误差e(k)(例如第二行第二个图)。微分的引入相当于一个超前校正,会减少超调,但是过渡的微分很可能会造成尾部振荡,系统逐渐变得不稳定。因此微分和积分之间需要一个平衡,当满足这个平衡的时候,系统几乎没有振荡,同时响应速度也较快。(第一行的图3是积分过多,产生超调,第二行的图1和图3就比较理想)
    综合上述,PID的调节经验可以归结为以下几点:

    • Kp较小时,系统对微分和积分环节的引入较为敏感,积分会引起超调,微分可能会引起振荡,而振荡剧烈的时候超铁也会增加。
    • Kp增大时,积分环节由于滞后产生的超调逐渐减小,此时如果想要继续减少超调可以适当引入微分环节。继续增大Kp系统可能会不太稳定,因此在增加Kp的同时引入Kd减小超调,可以保证在Kp不是很大的情况下也能取得较好的稳态特性和动态性能。
    • Kp较小时,积分环节不宜过大,Kp较大时积分环节也不宜过小(否则调节时间会非常地长),在下面这个例子中我们还会介绍到,当使用分段PID,在恰当的条件下分离积分,可以取得更好的控制效果。原因在于在稳态误差即将满足要求时,消除了系统的滞后。因此系统超调会明显减少。本例中采样的抗积分饱和的方法是遇限削弱积分法。

    4.改进PID算法(遇限削弱积分法)

    遇限削弱积分法的原理是
    u ( k ) > u m a x u(k)>u_{max} u(k)>umax时,若e(k)>0即输出值还未到达指定值,则认为积分会带来滞后,不再积分。
    u ( k ) < 0 u(k)<0 u(k)<0时,若e(k)<0即输出值超过了指定值,则认为积分会带来滞后,不再积分。
    在本案例中认为 u m a x = r ( k ) u_{max}=r(k) umax=r(k)
    改进PID算法如下(需要些两个循环,当然也可以用一个循环,将其中的PID设为一个子过程调用):

    close all
    ts=0.005;  %采样时间=0.005s
    sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数,即式4.1
    dsys=c2d(sys,ts,'z');      %离散化
    [num,den]=tfdata(dsys,'v');   %
    e_1=0;      %前一时刻的偏差      
    Ee=0;       %累积偏差
    u_1=0.0;    %前一时刻的控制量
    y_1=0;       %前一时刻的输出
    %PID参数
    kp=0.22;    
    ki=0.13;
    kd=0;
    u=zeros(1,1000);
    time=zeros(1,1000);
    for k=1:1:1000
        time(k)=k*ts;   %时间参数
        r(k)=1500;      %给定量
        y(k)=-1*den(2)*y_1+num(2)*u_1+num(1)*u(k);
        e(k)=r(k)-y(k);   %偏差
        u(k)=kp*e(k)+ki*Ee+kd*(e(k)-e_1);   
        Ee=Ee+e(k);    
        u_1=u(k);    
        y_1=y(k);    
        e_1=e(k);
    end
    p1=plot(time,r,'-.');xlim([0,1]);hold on;
    p2=plot(time,y,'--');xlim([0,1]);
    hold on;
    a=1;%控制积分分离的二值数
    e_1=0;Ee=0;u_1=0.0;y_1=0;%重新初始化       
    for k=1:1:1000
        time(k)=k*ts;   %时间参数
        r(k)=1500;      %给定量
        y(k)=-1*den(2)*y_1+num(2)*u_1;
        e(k)=r(k)-y(k);   %偏差
        u(k)=kp*e(k)+ki*Ee+kd*(e(k)-e_1);   
         if ((u(k)>r(k)) && (e(k)>0))||((u(k)<0) && (e(k)<0))
             a=0;
         else 
             a=1;
         end     
        Ee=Ee+a*e(k);    
        u_1=u(k);    
        y_1=y(k);    
        e_1=e(k);
    end
    p3=plot(time,y,'-');xlim([0,1]);
    title('含积分分离与不含积分分离的对比');
    legend([p1,p2,p3],'指令信号','不含积分分离','含积分分离');
    

    输出的曲线对比图如下:
    积分分离之后的改进PID
    可以发现,系统的超调量明显减少了,调节时间也减少了一点。原因在于我们采用了分段PID的手段,既消除了稳态误差还削弱了积分环节带来的滞后影响。

    5.simulink仿真

    需要的模块名称(不区分大小写)如下:

    • gain(参数分别为0.22和0.13/0.005)
    • sum(参数分别为"|±"和"|++")
    • integrator
    • scope
      注意:本文使用的是离散PID仿真,而simulink使用的是连续系统仿真,转换PID参数时P参数不变,I参数应该除以仿真间隔Ts=0.005,D参数应该乘Ts。

    以表中第一组PI参数为例:
    在这里插入图片描述
    得到的示波器曲线如下:
    在这里插入图片描述

    希望本文对您有帮助,谢谢阅读。

    展开全文
  • 本文提供粒子群算法简介和一个算法举例,提供粒子群算法仿真PID的M文件代码...另外,本文还提供了一种动态simulink仿真方法,可以让M文件和simulink文件之间互相交换数据,实现仿真与程序的反馈,增加了仿真的灵活度。


    阅读前必看:

    1. 本代码基于MATLAB2017a版本,如果版本不同可能会报错
    2. 请从set_para.m文件开始运行,其他M文件(+下载的资源包里面的slx文件)放在一个文件夹
    3. 每次迭代都会打印出来,如果运行时间过长怀疑程序死机可以观察迭代次数是否有变化
    4. 最后会输出3幅图(figure)
    5. slx文件下载链接:https://download.csdn.net/download/weixin_44044411/11609837

    混合粒子群算法链接:https://blog.csdn.net/weixin_44044411/article/details/103638611?spm=1001.2014.3001.5502
    MATALB的自带的粒子群优化器来了:particleswarm

    0.背景

    \qquad 调整PID参数时,我们常常需要按照各种方法查各种表,还要根据建立的数学模型和之前的结果做人工调整。这个过程常常是漫长的,而且最终得到的很可能也不是最优的参数。若系统的参数有改动,又要人工再调一次。而3个参数6个方向的调整,如果输入的变量很多,那么调节的可能性组合就更多了,这个工作常常不是非专业人士所能胜任的。
    \qquad 现在我给大家介绍一种不需要知道内部原理也可以调整PID参数的方法,而且调整的思路是所需即所得——你需要的参数是什么样子的,这个算法就会往这个方向调整PID参数。如果这个系统是多输入多输出的,我需要让A输出的超调小一些,而B输出的调节时间可以尽量大一些,只需要调整我们的评价函数就可以让算法自动调参以满足客户需求了。那这样的算法使用起来是不是很容易呢?

    1.粒子群算法

    1.1.算法简介

    \qquad 粒子群算法(PSO)于1995年提出,和遗传算法一样,也是一种群体迭代算法,不同的是,粒子群算法需要整定的参数更少,不存在交叉和变异过程,所以收敛速度更快。
    \qquad 由于是群体迭代算法,因此该算法一般做成并行的。如果做成串行的,也是可以运行的,就是会随着种群规模的增大越来越慢。对于Python运用熟练的朋友可以利用Python做成并行的,对于不熟悉并行算法的朋友,直接使用MATLAB就好了,因为MATLAB的矩阵运算已经集成了并行算法。
    \qquad 为了让大家理解,我们做一个很形象的比喻。将我们的粒子群比作一群可以互相通讯的鸟,每只鸟不知道食物的具体位置,但是可以根据嗅觉知道离食物的距离和自己在地图中的位置。种群中的任意2只鸟均可以知道对方离食物的距离和对方的具体位置。假设每只鸟的飞行都是有惯性的,而每只鸟都会根据自己距离食物的距离和位置及其他鸟距离食物的距离和位置选择下一次的飞行方向。最终,所有鸟都会在信息交换的基础上,汇聚在食物地点。
    \qquad 这只是对粒子群算法的一个形象的比喻,不必过多计较该模型的现实性,闲言少叙,让我们进入 算 法 步 骤 \color{red}算法步骤

    1.2.算法步骤

    【step1】确定参数维度 N \color{blue}{N} N、参数范围(即每只鸟的初始位置),确定惯性系数 c 1 , c 2 , w c_1,c_2,w c1,c2,w,确定种群规模m,迭代次数n。
    每个粒子的初始位置是随机的,设输入向量 x = ( x 1 , x 2 , . . . , x N ) T x=(x_1,x_2,...,x_N)^T x=(x1,x2,...,xN)T必须满足参数范围限制为:
    { x m i n ( 1 ) < x 1 < x m a x ( 1 ) x m i n ( 2 ) < x 2 < x m a x ( 2 ) . . . x m i n ( N ) < x 1 < x m a x ( N ) \begin{cases}x_{min}^{(1)}<x_1<x_{max}^{(1)} \\x_{min}^{(2)}<x_2<x_{max}^{(2)} \\... \\x_{min}^{(N)}<x_1<x_{max}^{(N)} \end{cases} xmin(1)<x1<xmax(1)xmin(2)<x2<xmax(2)...xmin(N)<x1<xmax(N)
    X m i n = ( x m i n ( 1 ) , x m i n ( 2 ) , . . . x m i n ( N ) ) , X m a x = ( x m a x ( 1 ) , x m a x ( 2 ) , . . . x m a x ( N ) ) \color{blue}X_{min}=(x_{min}^{(1)},x_{min}^{(2)},...x_{min}^{(N)}),X_{max}=(x_{max}^{(1)},x_{max}^{(2)},...x_{max}^{(N)}) Xmin=(xmin(1),xmin(2),...xmin(N)),Xmax=(xmax(1),xmax(2),...xmax(N))
    则输入向量 x x x应满足 X m i n < x < X m a x X_{min}<x<X_{max} Xmin<x<Xmax
    每个粒子的初速度定为0,即 v 0 = 0 v_0=0 v0=0,第 j j j个粒子( 1 ≤ j ≤ m 1≤j≤m 1jm)的下一次迭代的速度 v ( j ) v^{(j)} v(j)由三部分组成:
    v ( j ) = w ⋅ v 0 + c 1 ⋅ r a n d ⋅ ( P ( j ) − X ( j ) ) + c 2 ⋅ r a n d ⋅ ( P G − X ( j ) ) v^{(j)}=w\cdot v_0+c_1\cdot rand\cdot (P^{(j)}-X^{(j)})+c_2\cdot rand\cdot(P_G-X^{(j)}) v(j)=wv0+c1rand(P(j)X(j))+c2rand(PGX(j))
    注:rand是(0,1)的随机数, v 0 v_0 v0代表上一次粒子的速度。
    第一部分为自身惯性因子,因为下一次的迭代次数保留了上一次的速度信息;
    第二个部分为自身最优因子, P ( j ) \color{blue}{P^{(j)}} P(j)为第 j j j个因子在之前所有迭代次数中自适应度最高的位置,可以理解为历史上自身离食物最近的位置。
    第三部分为社会因子, P G \color{blue}{P_G} PG为种群在之前所有迭代次数中自适应度最高的位置,可以理解为历史上所有粒子离食物最近的位置中的最近的一个。
    一般情况下,取 w = 1 , c 1 = c 2 = 2 w=1,c_1=c_2=2 w=1,c1=c2=2,当种群规模较大时,可以让 w w w的值随迭代次数减小以增加收敛速度。
    【step2】按照step1的速度公式计算下一次的速度,并计算下一次的粒子位置。对于第 j j j个粒子,第 k + 1 k+1 k+1次迭代(第 k + 1 k+1 k+1代)的位置 X k + 1 ( j ) \color{blue}{X_{k+1}^{(j)}} Xk+1(j)与第 k k k次迭代的位置 X K ( j ) \color{blue}{X_K^{(j)}} XK(j)与速度 v k ( k + 1 ) \color{blue}{v_k^{(k+1)}} vk(k+1)关系为:
    X k + 1 ( j ) = X k ( j ) + v k ( j + 1 ) ⋅ d t X^{(j)}_{k+1}=X^{(j)}_{k}+v_k^{(j+1)}\cdot dt Xk+1(j)=Xk(j)+vk(j+1)dt
    其中 d t \color{blue}{dt} dt是仿真间隔,一般情况下可以取1,如果希望仿真得慢一点,搜索平滑一点,可以适当减小 d t \color{blue}{dt} dt
    【step3】计算每个粒子的自适应度 F k + 1 ( j ) \color{blue}{F^{(j)}_{k+1}} Fk+1(j),为了计算出step1公式中的 P ( j ) 和 P G \color{blue}{P^{(j)}和P_G} P(j)PG。为了方便起见,我们记前k次计算得到了的 P G P_G PG P G ( k ) P_G^{(k)} PG(k),则第k+1次迭代中我们将适应度最高的粒子位置记为 P G ( k + 1 ) P_G^{(k+1)} PG(k+1),则最终的 P G P_G PG为:
    P G = { P G ( k ) , F ( P G ( k ) ) > F ( P G ( k + 1 ) ) P G ( k + 1 ) , F ( P G ( k ) ) < F ( P G ( k + 1 ) ) P_G=\begin{cases}P_G^{(k)},\qquad F(P_G^{(k)})>F(P_G^{(k+1)}) \\[2ex]P_G^{(k+1)},\quad F(P_G^{(k)})<F(P_G^{(k+1)}) \end{cases} PG=PG(k)F(PG(k))>F(PG(k+1))PG(k+1)F(PG(k))<F(PG(k+1))
    同样,我们记前k次计算得到的第 j j j个粒子的位置为 P k ( j ) P^{(j)}_{k} Pk(j),第k+1次计算得到的第 j j j个粒子的位置为 P k + 1 ( j ) P^{(j)}_{k+1} Pk+1(j),则最终的第 j j j个粒子的历史最优解 P ( j ) P^{(j)} P(j)为:
    P ( j ) = { P k ( j ) , F ( P k ( j ) ) > F ( P k + 1 ( j ) ) P k + 1 ( j ) , F ( P k ( j ) ) < F ( P k + 1 ( j ) ) P^{(j)}=\begin{cases}P_{k}^{(j)},\quad F(P_{k}^{(j)})>F(P_{k+1}^{(j)}) \\[2ex]P_{k+1}^{(j)},\quad F(P_{k}^{(j)})<F(P_{k+1}^{(j)}) \end{cases} P(j)=Pk(j)F(Pk(j))>F(Pk+1(j))Pk+1(j)F(Pk(j))<F(Pk+1(j))
    【step4】更新每个粒子的信息。
    由于我们在step2的位置迭代中已经更新过粒子的位置信息,在step1中的速度迭代中已经更新过粒子的速度信息,而在step3中又更新了每个粒子的历史最优位置 P ( j ) P^{(j)} P(j)及种群最优位置 P G P_G PG,因此实际上如果仅需要知道最优解的情况下我们不需要做这一步。
    但如果需要作出参数随迭代次数的改变的图的话,每次迭代产生最优解 P G ( k ) P_G^{(k)} PG(k)及最优适应度 F ( P G ( k ) ) F(P_G^{(k)}) F(PG(k))需要用数组保存下来。

    1.3.算法举例

    \qquad 根据以上的算法步骤,本人编写了一段粒子群算法的MATLAB代码,这个代码很简单,是无约束问题下的全局粒子群算法。全局粒子群算法收敛速度快,占用内存小,但是容易陷入局部最优解,每次迭代的最优解可能存在细微差异,当然如果这个无约束问题本身是单解的话,这个差异会随着种群规模的扩大和迭代次数增加而逐渐消除。

    function [Pg,fmin]=PSO(f,dimension,n,m,xmax,xmin,vmax,vmin)
    %全局粒子群算法,f为目标函数,dimension为维度,n为代数,m为种群规模
    %位置限幅为[xmin,xmax],速度限幅为[vmin,vmax]
        Savef=zeros(n+1,1);
        SaveData=zeros(m,dimension,10);%记录11代种群的位置
        w=1;c1=2;c2=2;%速度惯性系数
        dt=0.3;%位移仿真间隔
        v=zeros(m,dimension);%初始速度为0
        X=(xmax-xmin)*rand(m,dimension)+xmin;%初始位置满足(-xmax,xmax)内均匀分布
        P=X;%P为每个粒子每代的最优位置
        last_f=f(X);
        [fmin,min_i]=min(last_f);%Pg为所有代中的最优位置 
        Pg=X(min_i,:);
        Savef(1)=fmin;
        N=0;
        for i=1:n
            v=w*v+c1*rand*(P-X)+c2*rand*(ones(m,1)*Pg-X);
            v=(v>vmax).*vmax+(v>=vmin & v<=vmax).*v+(v<vmin).*vmin;
            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);%更新函数最小值
             Savef(i)=fmin;
             if mod(i,floor(n/10))==0%10代记录一次种群参数
                 N=N+1;
                SaveData(:,:,N)=X;
             end
    %         w=w-i/n*0.7*w;
        end
        for j=1:10
            figure(j)
            plot(SaveData(:,1,j),SaveData(:,2,j),'o');
            xlim([-xmax,xmax])
            ylim([-xmax,xmax])
            axis tight
        end
        figure
        plot(Savef','b-')
        disp(Pg)
        disp(fmin)
    end
    

    这是一个粒子群算法的函数,可以求函数在一定范围内的极小值,下面我们来定义一个测试函数测试一下:
    f = ( 1 − x 1 ) 2 + ( x 1 2 − 2 x 2 ) 2 f=(1-x_1)^2+(x_1^2-2x_2)^2 f=(1x1)2+(x122x2)2
    我们很容易知道这个函数的极小值为 x ∗ = ( 1 , 0.5 ) T x^*=(1,0.5)^T x=(1,0.5)T
    我们将X和Y方向上的限制都设定为[-40,40],种群规模200,仿真代数100次

    f=@(x)(1-x(:,1)).^2+(x(:,1).^2-2*x(:,2)).^2;
    [Pg,fmin]=PSO(f,2,100,200,40,-40,40,-40)
    

    运行结果为:
    Alt

    序号粒子群位置随迭代次数的变化
    1Alt
    2在这里插入图片描述
    3在这里插入图片描述
    4在这里插入图片描述
    5在这里插入图片描述
    6在这里插入图片描述

    \qquad 读者们可以发现粒子群算法在迭代的过程中部分粒子甚至有反向运动的趋势,这是由于粒子的惯性速度过大导致的,在迭代点附近,惯性因子过大,会导致很长时间无法收敛到一个点,同时也很难找到相比上一次迭代更优的解。
    \qquad 我们只需要随着迭代次数增加减少我们的惯性因子 w w w即可,在这里作者给出的表达式为
    w k + 1 = w k − i / n ∗ 0.7 w k w_{k+1}=w_k-i/n*0.7w_k wk+1=wki/n0.7wk
    i i i为当前迭代代数, n n n为总代数。
    调用同样的测试函数,种群规模、迭代次数及其他参数不变,我们得到结果如下:
    在这里插入图片描述
    和上一次对比我们可以发现fmin小了很多,即在相同的计算规模下我们的最优解更加接近于极小值,可谓一个突破性的改变,我们再看看粒子随迭代次数的分布:

    序号粒子群位置随迭代次数的变化
    1在这里插入图片描述
    2在这里插入图片描述
    3在这里插入图片描述

    \qquad 可以发现粒子群的分布相比之前稳定得多,几乎全部收敛到了中心点,这便是减小惯性因子的作用。

    2.PID自整定

    2.1.基于M文件编写的PID参数自整定

    \qquad 利用M文件可以制作一个简单的PID引擎,用于测试我们的PSO算法是否有效,书写简单PID引擎的代码可以参考如下链接:
    【MATLAB】仿真PID控制
    \qquad 我们仅需要在此基础上做一些改变并做成并行算法即可,当然,我们还需要设置一个自适应度函数,在自适应度函数的作用下,PSO算法会选择到我们满意的参数。在调整PID参数时,我们通常会关注响应曲线的超调、上升时间、调节时间、峰值时间等等。利用这些参数组合为评价函数,让评价函数越小,仿真得到的参数越好,以此作为解的好坏的对比依据(相当于鸟离食物的距离)为了简便起见,我们只选取了调节时间 t s t_s ts和超调量 σ \sigma σ作为我们的比较依据,我们将评价函数设定为:
    f a c c e s s = l n ( t s 5 ⋅ 1 0 − 2 + 1 ) + l n ( σ 1 0 − 4 + 1 ) f_{access}=ln(\frac{t_s}{5\cdot 10^{-2}}+1)+ln(\frac{\sigma}{10^{-4}}+1) faccess=ln(5102ts+1)+ln(104σ+1)
    如果对该函数求导就可以发现:
    ∂ f ∂ t s ∣ t s = 5 ⋅ 1 0 − 2 = 1 2 ⋅ 5 ⋅ 1 0 − 2 \frac{\partial f}{\partial t_s}_{|t_s=5\cdot 10^{-2}}=\frac{1}{2}\cdot 5\cdot 10^{-2} tsfts=5102=215102
    ∂ f ∂ σ ∣ σ = 1 0 − 4 = 1 2 ⋅ 1 0 − 4 \frac{\partial f}{\partial \sigma}_{|\sigma=10^{-4}}=\frac{1}{2}\cdot 10^{-4} σfσ=104=21104
    \qquad 因此 f a c c e s s f_{access} faccess t s = 5 ⋅ 1 0 − 2 t_s=5\cdot 10^{-2} ts=5102的调整率为 2.5 × 1 0 − 2 2.5×10^{-2} 2.5×102,在 σ = 1 0 − 4 \sigma=10^{-4} σ=104时的调整率为 0.5 × 1 0 − 4 0.5×10^{-4} 0.5×104如果将期望参数定为 t s ∗ = 5 ⋅ 1 0 − 2 , σ ∗ = 1 0 − 4 t_s^*=5\cdot 10^{-2},\sigma ^*=10^{-4} ts=5102,σ=104则评价函数在小于任意一个参数的情况下就会下降,但是下降的速度会越来越大,而大于任意一个参数就会上升,但上升速率会越来越慢。由于两个参数均为正数,因此下降的速率不会大于 ∂ f ∂ t s ∣ t s = 0 = 1 \frac{\partial f}{\partial t_s}_{|t_s=0}=1 tsfts=0=1 ∂ f ∂ σ ∣ σ = 0 = 1 \frac{\partial f}{\partial \sigma}_{|\sigma=0}=1 σfσ=0=1,而上升的速率可以无限小。若将评价函数小视为自适应度高,则调整过后的系统参数总是在期望值附近调整,并尽量满足$所有参数都小于期望值的情况。因为在期望值附近上升速率是越来越慢而下降速率是越来越快的,而参数最小只能取到0。因此,有多个参数大于期望值而有少个参数小于期望值的和会比有少个参数大于期望值而多个参数小于期望值的和要大得多。
    \qquad 写好评价函数之后,我们可以设自适应度函数 F = − f a c c e s s F=-f_{access} F=faccess,因为我们的评价函数是越小越好的,而自适应度是越高越好。
    \qquad 根据1.2的算法步骤我们可以编写出粒子群算法,而PID并行仿真引擎的编写需要参考上文提供的连接。(simulink的PID仿真请参考下一小节2.2.)粒子群算法每迭代出一代的参数(PSO.m),就把这一代里面最好的参数 P G P_G PG返还给PID引擎进行仿真(PID_sim.m),再编写一个参数自识别的程序(parameter_cal.m),将需要识别的指标(调节时间、超调)计算出来,传递给 f a c c e s s f_{access} faccess进行计算(PID_access.m),将计算得到的参数

    【PID_sim.m】(PID并行仿真引擎)

    function [f_infty,tp,ts,sigma]=PID_sim(kp,ki,kd,debug)
        %kp,ki,kd为PID参数,T0为采样时间,total_t为仿真时间
        Tt=5e-4;%仿真采样周期
        T0=1e-2;%控制采样周期
        Tf=1e-3;%微分环节的滤波器系数
        alp=Tf/(Tt+Tf);
        total_t=1;
        N=floor(total_t/Tt);%样本总数
        M=numel(kp);%仿真序列数
        kp=reshape(kp,M,1);
        ki=reshape(ki,M,1);
        kd=reshape(kd,M,1);
        sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数,即式4.1
        dsys=c2d(sys,Tt,'z');      %离散化
        [num,den]=tfdata(dsys,'v');   %
        e_1=zeros(M,1);      %前一时刻的偏差      
        Ee=zeros(M,1);       %累积偏差
        u_1=zeros(M,1);    %前一时刻的控制量
        y_1=zeros(M,1);       %前一时刻的输出
        ud_1=zeros(M,1);     %前一时刻的微分输出
        %预分配内存
        r=1500;
        y=zeros(M,N);
        u=zeros(M,N);
        for k=1:1:N
            y(:,k)=-1*den(2).*y_1+num(2)*u_1+num(1).*u(:,k);%系统响应输出序列
            e0=r-y(:,k);   %误差信号
            ud=alp.*ud_1+(1-alp)/T0.*kd.*(e0-e_1);
            u(:,k)=kp.*e0+T0*ki.*Ee+ud; %系统PID控制器输出序列
            Ee=Ee+e0;    %误差的累加和
            u_1=u(:,k);    	%前一个的控制器输出值
            y_1=y(:,k);    	%前一个的系统响应输出值
            e_1=e0;		%前一个误差信号的值
            ud_1=ud;
        end
        if debug %非调试模式下不显示也不打印图像
            plot(linspace(0,total_t,N),y)
        end
        [f_infty,tp,ts,sigma]=parameter_cal(y,Tt,2e-2,debug);%取得阶跃响应指标
    end
    

    【parameter_cal.m】(计算仿真曲线参数)

    function [f_infty,tp,ts,sigma]=parameter_cal(y,T0,delt_err,debug)
    %y是系统输出序列
    %T0是采样时间,可以结合时间序列点序号算出实际时间
    %delt_err是调节时间的误差精度
        M=size(y,1);N=size(y,2);%M为运算维度,N为时间序列长度
        f_infty=y(:,N);%稳态值序列
        err=y-f_infty*ones(1,N);%通过稳态值计算误差序列
        ferr=fliplr(abs(err));%倒序并取绝对值
        [~,ts_i]=max(ferr>delt_err*f_infty,[],2);
        ts_i=N*ones(M,1)-ts_i;
        ts=ts_i*T0;%调节时间
        [fp,tp]=max(y,[],2);%峰值和函数峰值
        tp=tp*T0;
        tp(abs(fp-f_infty)<1e-5)=NaN; %过阻尼无超调,没有峰值时间
        sigma=(fp-f_infty)./f_infty;
        if debug && M==1 %非调试模式下不显示
            disp(['系统稳态值为',num2str(f_infty)])
            disp(['系统超调量为',num2str(sigma*100),'%'])
            if isnan(tp)
                disp('系统不存在峰值时间')
            else
                disp(['系统峰值时间为',num2str(tp),'s'])
            end
            disp(['系统的调节时间为',num2str(ts),'s'])
        end
    end
    

    【PID_access.m】(根据参数得出评价度的函数,程序中为评价度越低越好)

    function y=PID_access(para)
    %PID调节性能的指标参数
    kp=para(:,1);ki=para(:,2);kd=para(:,3);
    [~,~,ts,sigma]=PID_sim(kp,ki,kd,false);
    y=log(ts/5e-2+1)+log(sigma/1e-2+1);
    end
    

    【PSO.m】(PSO算法优化函数)

    function [Pg,fmin]=PSO(dimension,n,m)
    %全局粒子群算法,f为目标函数,dimension为维度,n为代数,m为种群规模
        w=1;c1=2;c2=2;%速度惯性系数
        sigma_data=zeros(1,n);
        ts_data=zeros(1,n);
        dt=0.3;%位移仿真间隔
        vmax=1;%速度限幅
        xmax=[15,50,2];%位置限幅
        xmin=[0.2,0,0];
        v=zeros(m,dimension);%初始速度为0
        X=(xmax-xmin).*rand(m,dimension)+xmin;%初始位置满足(xmin,xmax)内均匀分布
        P=X;%P为每个粒子每代的最优位置
        last_f=PID_access(X);
        [fmin,min_i]=min(last_f);%Pg为所有代中的最优位置 
        Pg=X(min_i,:);
        N=0;
        for i=1:n
            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=PID_access(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,sigma]=PID_sim(Pg(1),Pg(2),Pg(3),true)
            sigma_data(1,i)=sigma;
            ts_data(1,i)=ts;
            hold on
        end
        legend('1','2','3','4','5','6','7','8','9','10','11','12','13','14','15')
        title('响应曲线随迭代次数的变化')
        figure
        plot(ts_data)
        title('调节时间随迭代次数的变化')
        figure
        plot(sigma_data)
        title('超调量随迭代次数的变化')
    end
    

    【main.m】(主函数,从这里开始)

    dimension=3;n=20;m=200;%PID为3位参数,n是迭代次数,m为种群规模
    [Pg,fmin]=PSO(dimension,n,m)
    

    仿真结果如下:
    在这里插入图片描述
    注:响应曲线可能要在左上角局部放大才能获得该图像。
    在这里插入图片描述
    在这里插入图片描述
    \qquad 这是一个简单的一阶惯性系统的PID调参仿真,用于模拟直流电机的转速模型,可以看出直流电机的转速响应是越来越快的(调节时间越来越小),超调量在第一次就达到了0,因此系统不存在超调。

    *2.2.复杂系统的PID自整定(基于simulink仿真)

    \qquad 我们以双闭环直流调速系统为例来说明PSO对PID自整定参数的作用。直流调速系统由外环(转速环)和内环(电流环)构成,转速环是一个PI调节器,调节转速为期望转速,用于抵抗负载扰动、电机参数变化带来的扰动;电流环也是一个PI调节器,与转速环不同的是,电流调节器的目的是在于让电流快速跟随给定值,以获得更好的动态特性。我们以电机的启动过程为例仿真这个系统,直流电机启动时,电流首先达到容许最大值,电机以最大加速度进行加速,到达期望转速时让电流迅速减小并达到负载转矩所需电流,电机维持期望转速不变。

    2.2.1.PSO优化PID的过程详解

    \qquad 下面我们设定我们需要达成的调节指标——由于启动时的电机电流会迅速上升到过载电流(允许的最大值),是一个恒定值,故会有一个电流平台,记上升到电流平台产生的超调为电流超调 σ i \sigma_i σi,从超调回到电流平台所需的时间为电流的调节时间 t s i t_{si} tsi,最终转速稳定前的超调记为转速超调 σ n \sigma_n σn,转速到达稳定值区域的时间为调节时间 t s n t_{sn} tsn
    【搭建simulink模型】
    双闭环调速系统的结构框图如下:
    在这里插入图片描述
    搭建simulink仿真模型如下:
    在这里插入图片描述
    注:上图可能不够清楚,根据以下参数搭建simulink模型,并且命名为double_circle.slx(和PSO的所有程序同名)
    电机的参数如下:

    • 额定电压/电流/功率:800V/2164A 1500KW
    • 额定转速:44r/min
    • 过载系数: λ i = 1.5 \lambda_i=1.5 λi=1.5
    • 励磁:110V/191A
    • 电枢电阻 R a R_a Ra:0.1 Ω \Omega Ω
    • 主回路电阻R:0.1014 Ω \Omega Ω
    • 电枢(主回路)电感L:3.25mH
    • 变流器放大倍数 K s K_s Ks=80
    • 变流器通态压降:3V
    • 系统飞轮惯量: 8.6 × 1 0 5 N ⋅ m 2 8.6×10^5N\cdot m^2 8.6×105Nm2

    \qquad 根据工程设计方法当然可以得到一套电流和转速的PI参数,计算的相关参数为: [ K n p , K n I , K i p , K i I ] E N = [ 7.3660 , 1.5658 , 0.8257 , 30.2773 ] [K_{np},K_{nI},K_{ip},K_{iI}]_{EN}=[7.3660 ,1.5658 ,0.8257, 30.2773] [Knp,KnI,Kip,KiI]EN=[7.3660,1.5658,0.8257,30.2773]
    通过仿真,我们得到的转速和电流曲线如下:

    额定负载启动空载启动
    在这里插入图片描述在这里插入图片描述

    调用我们的参数计算程序计算一下参数可以得到: σ n = 4.07 × 1 0 − 4 % , t s n = 0.7875 , σ i = 20.38 % , t s i = 0.0172 \sigma_n=4.07×10^{-4}\%,t_{sn}=0.7875,\sigma_i=20.38\%,t_{si}=0.0172 σn=4.07×104%,tsn=0.7875,σi=20.38%,tsi=0.0172
    \qquad 这个参数当然是比较理想的参数,如果我们不知道工程设计方法或者说工程设计方法不适用于我们的模型时,那么就可以使用PSO整定PID参数的方法。
    【修改simulink模型以适应PSO的M文件调用】
    \qquad 这里需要调整的参数为转速调节器ASR的 K n p , K n I K_{np},K_{nI} KnpKnI和电流调节器ACR的 K i p , K i I K_{ip},K_{iI} Kip,KiI,如果我们希望能实时调整simulink的参数进行动态仿真,我们需要直接将图中的几个比例环节依次设为符号( K n p , K n I , K i p , K i I K_{np},K_{nI},K_{ip},K_{iI} Knp,KnI,Kip,KiI),并在工作区实时修改这4个参数再进行仿真。如果我们要将实时调整PI参数的代码写在M文件中仿真,必须要将变量声明为全局(global),此时变量会在工作区里面,simulink就可以调用了。另存为PSO_double_circle.slx
    【编写辅助的M文件】
    cal_n.m:单组PID参数仿真结果的参数计算

    function [n_sigma,ts_n,I_sigma,ts_transient_i,n_err]=cal_n(parameters)
        global Kpn KIn Kpi KIi
        %将simulink所需参数导入工作区,以便仿真使用
        obj_slx='PSO_double_circle.slx';
        Kpn=parameters(1);
        KIn=parameters(2);
        Kpi=parameters(3);
        KIi=parameters(4);
        sim(obj_slx);
        I_data=I.signals.values;n_data=n.signals.values;time=I.time;
        n_infty=mean(n_data(numel(n_data)-10:numel(n_data)));
        n_err=abs(n_infty-10/0.227);%期望与输出的转速误差
        [i_max,tr_index]=max(I_data);
        diffI=[0;diff(I_data)];%电流微分
        pos=find(abs(diffI)<1);%找到电流波动较小的部分
        pos=pos(pos>tr_index);%排除启动时电流不变的情况
        min_pos=min(pos);
        I_transient=mean(I_data(min_pos+10:min_pos+20));%恒流升速阶段的暂稳态电流值
        I_sigma=(i_max-I_transient)/I_transient;
        ts_transient_i=time(min_pos-tr_index);%电流暂稳态的调节时间
        n_max=max(n_data);
        n_sigma=(n_max-n_infty)/n_infty;
        ts_n_index=numel(n_data)-find(abs(flip(n_data)-n_infty)/n_infty>2e-2,1);
        ts_n=time(ts_n_index);
    end
    

    PID_n_access.m:PID参数的仿真,以代为单位的

    function [y,ts_n,n_sigma,ts_transient_i,I_sigma]=PID_n_access(parameters_list)
        M=size(parameters_list,1);%计算需要仿真的组数(simulink无法并行仿真)
        ts_n0=0.8;n_sigma0=1e-2;I_sigma0=0.02;ts_i0=1e-4;%期望值
        %将simulink所需参数导入工作区,以便仿真使用
        y=zeros(M,1);
        n_sigma=zeros(M,1);ts_n=zeros(M,1);
        for sim_i=1:M
            parameters=parameters_list(sim_i,:);
            [n_sigma,ts_n,I_sigma,ts_transient_i,n_err]=cal_n(parameters);
            y(sim_i)=log(ts_n./ts_n0+1)+log(n_sigma./n_sigma0+1)+log(I_sigma/I_sigma0+1)+log(ts_transient_i/ts_i0+1)+exp(n_err*10);
        end
    end
    

    para_cal.m:参数计算,以代为单位的

    obj_slx='double_circle.slx';
    sim(obj_slx);
    I_data=I.signals.values;n_data=n.signals.values;time=I.time;
    I_infty=mean(I_data(numel(I_data)-10:numel(I_data)));
    n_infty=mean(n_data(numel(n_data)-10:numel(n_data)));
    [i_max,tr_index]=max(I_data);
    diffI=[0;diff(I_data)];%电流微分
    pos=find(abs(diffI)<1);%找到电流波动较小的部分
    pos=pos(pos>tr_index);%排除启动时电流不变的情况
    min_pos=min(pos);
    I_transient=mean(I_data(min_pos+10:min_pos+20));%恒流升速阶段的暂稳态电流值
    I_sigma=(i_max-I_transient)/I_transient;
    if I_sigma<1e-4
        tr_i=NaN;    
    else
        tr_i=time(tr_index);%电流的上升时间
    end
    ts_transient_i=time(min_pos-tr_index);%电流暂稳态的调节时间
    ts_i_index=numel(I_data)-find(abs(flip(I_data)-I_infty)/I_infty>2e-2,1);
    ts_i=time(ts_i_index);%转速稳态时电流的调节时间
    [n_max,tr_index]=max(n_data);
    n_sigma=(n_max-n_infty)/n_infty;
    %参数展示
    %转速参数展示
    if n_sigma<1e-4
        tr_n=NaN;
    else
        tr_n=time(tr_index);
    end
    ts_n_index=numel(n_data)-find(abs(flip(n_data)-n_infty)/n_infty>2e-2,1);
    ts_n=time(ts_n_index);
    if isnan(tr_n)
        disp('转速无超调')
    else
        disp(['转速上升时间为',num2str(tr_n),'s'])
        disp(['转速超调量为',num2str(100*n_sigma),'%'])
    end
    disp(['转速调节时间为',num2str(ts_n),'s'])
    %电流参数展示
    if isnan(tr_i)
        disp('电流无超调')
    else
        disp(['电流上升时间为',num2str(tr_i),'s'])
        disp(['电流超调量为',num2str(100*I_sigma),'%'])
        disp(['电流暂稳态调节时间为',num2str(ts_transient_i),'s'])
    end
    disp(['电流的最终调节时间为',num2str(ts_i),'s'])
    

    draw.m:仅根据一组PID参数画出转速波形图

    function draw(parameters)
        global Kpn Kin KpI KiI
        Kpn=parameters(1);
        Kin=parameters(2);
        KpI=parameters(3);
        KiI=parameters(4);
        sim('PSO_double_circle')
        n_data=n.signals.values;time=n.time;%获取转速序列及对应时刻
        plot(time,n_data)
    

    draw_n_I.m:根据一组PID参数绘制电流和转速曲线

    function draw_n_I(parameters)
        global Kpn KIn Kpi KIi
        Kpn=parameters(1);
        KIn=parameters(2);
        Kpi=parameters(3);
        KIi=parameters(4);
        sim('PSO_double_circle')
        n_data=n.signals.values;time=n.time;%获取转速序列及对应时刻
        I_data=I.signals.values;
        figure
        subplot(211)
        plot(time,I_data,'r-','LineWidth',1.5)
        title('电流波形')
        subplot(212)
        plot(time,n_data,'b-','LineWidth',1.5)
        title('转速波形')
    end
    

    PSO.m:PSO优化程序

    function [Pg,fmin]=PSO(n,m,xmax,xmin,cal_f,f)
    %全局粒子群算法,f为目标函数,dimension为维度,n为代数,m为种群规模
        dimension=numel(xmax);
        w=1;c1=2;c2=2;%速度惯性系数
        sigma_data=zeros(1,n);
        ts_data=zeros(1,n);
        dt=0.3;%位移仿真间隔
        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;
        legend_str=cell(0);
        legend_i=1;
        figure(1)
        legend_str{legend_i}=num2str(1);
        draw(Pg)
        for i=1:n
            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);%更新函数最小值
            [sigma,ts,~,~]=cal_f(Pg);
            sigma_data(1,i)=sigma;
            ts_data(1,i)=ts;
            if last_Pg~=Pg
                legend_i=legend_i+1;
                figure(1)
                legend_str{legend_i}=num2str(i);
                draw(Pg)
                hold on
            end
            last_Pg=Pg;
            disp(['迭代次数:',num2str(i)])
        end
        figure(1)
        legend(legend_str)
        title('响应曲线随迭代次数的变化')
        figure(2)
        plot(ts_data)
        title('调节时间随迭代次数的变化')
        figure(3)
        plot(sigma_data)
        title('超调量随迭代次数的变化')
    end
    

    set_para.m(设置参数程序,从这里开始)

    clear
    best_parameters=[8.4725,97,72,1.794,56.06];
    global Kpn KIn Kpi KIi
    Kpn=0;KIn=0;Kpi=0;KIi=0;
    n=10;m=100;xmax=[20,120,10,60];xmin=[0,0,0,0];cal_f=@cal_n;f=@PID_n_access;
    [Pg,fmin]=PSO(n,m,xmax,xmin,cal_f,f);
    
    

    我们将期望的参数设定为
    σ n ∗ = 0 , 1 % , t s n ∗ = 1 , σ i ∗ = 20 % , t s i ∗ = 0.005 \sigma_n^*=0,1\%,t_{sn}^*=1,\sigma_i^*=20\%,t_{si}^*=0.005 σn=0,1%,tsn=1,σi=20%,tsi=0.005利用PSO算法构造的评价函数为
    f a c c e s s = l n ( t s n 1 + 1 ) + l n ( σ n 0.001 + 1 ) + l n ( σ i 0.2 + 1 ) + l n ( t s i 0.005 + 1 ) ; f_{access}=ln(\frac{ts_n}{1}+1)+ln(\frac{\sigma_n}{0.001}+1)+ln(\frac{\sigma_i}{0.2}+1)+ln(\frac{ts_i}{0.005}+1); faccess=ln(1tsn+1)+ln(0.001σn+1)+ln(0.2σi+1)+ln(0.005tsi+1);
    得到的最佳参数为:
    [ K n p , K n I , K i p , K i I ] P S O = [ 18.2613 , 73.0965 , 1.6843 , 16.5703 ] P S O [K_{np},K_{nI},K_{ip},K_{iI}]_{PSO}=[18.2613,73.0965 ,1.6843,16.5703]_{PSO} [Knp,KnI,Kip,KiI]PSO=[18.2613,73.0965,1.6843,16.5703]PSO
    下的仿真曲线如下:
    超调量和调节时间随迭代次数变化
    在这里插入图片描述在这里插入图片描述
    PSO与工程设计方法电机启动曲线对比(额定负载启动)

    工程设计方法PSO方法
    在这里插入图片描述在这里插入图片描述

    通过计算程序得到的参数为: σ n = 6.7556 × 1 0 − 4 % , t s n = 0.7862 , σ i = 18.83 % , t s i = 0.01341 \sigma_n=6.7556×10^{-4}\%,t_{sn}=0.7862,\sigma_i=18.83\%,t_{si}=0.01341 σn=6.7556×104%,tsn=0.7862,σi=18.83%,tsi=0.01341
    工程设计方法的参数为:
    σ n = 4.07 × 1 0 − 4 % , t s n = 0.7875 , σ i = 20.38 % , t s i = 0.0172 \sigma_n=4.07×10^{-4}\%,t_{sn}=0.7875,\sigma_i=20.38\%,t_{si}=0.0172 σn=4.07×104%,tsn=0.7875,σi=20.38%,tsi=0.0172
    \qquad 可以发现除了最后一项指标之外,其余指标均达到了期望值要求,并且和工程设计方法得出的指标相比甚至还略优于它。

    2.2.2.在PSO优化过程中修改参数价值权重

    \qquad 转速的超调和调节时间,电流的超调和调节时间,我们最后评估的指标除了转速无静差的基本要求外,这4项指标就是我们的评价标准,如果我们更希望牺牲某一个或多个参数去优化另一个或另一些参数时,该怎么调节我们的算法呢?
    \qquad 答案也非常简单,调节的方法有2个,一个是期望值,一个是评价函数的各项系数。这里我们用调节期望值的方法做一个示范:
    \qquad 之前的期望值为:
    σ n ∗ = 0.1 % , t s n ∗ = 1 , σ i ∗ = 20 % , t s i ∗ = 0.005 \sigma_n^*=0.1\%,t_{sn}^*=1,\sigma_i^*=20\%,t_{si}^*=0.005 σn=0.1%,tsn=1,σi=20%,tsi=0.005
    修改期望值为:
    σ n ∗ = 10 % , t s n ∗ = 0.8 , σ i ∗ = 2 % , t s i ∗ = 0.0001 \sigma_n^*=10\%,t_{sn}^*=0.8,\sigma_i^*=2\%,t_{si}^*=0.0001 σn=10%,tsn=0.8,σi=2%,tsi=0.0001
    重新仿真,对比二者的结果如下:

    PSO1PSO2
    在这里插入图片描述在这里插入图片描述

    \qquad 可以看出转速超调量略微增大了一些,但是转速和电流的调节时间明显缩短,电流的超调量基本不变。这一点也表明多项优化指标是相互耦合的,可以牺牲一者满足另一者,只要参数在这方面是可调的。
    \qquad 需要注意的是,simulink中无法使用串行的算法,因此代码中使用时是串行的,可能仿真比较慢,请大家耐心等待(每迭代一次都会打印一次)
    \qquad 希望本文对您有帮助,谢谢阅读!

    展开全文
  • LDO参数解读、特性、参考设计

    万次阅读 多人点赞 2019-12-05 18:34:19
    这篇博客讲述LDO相关知识,包括基本原理,参数解释,特性及参考设计等。
    展开全文
  • 3 PI控制器参数整定 3.1从PMSM电机的数学模型出发。 dq 轴 电压方程: dq 轴 轴磁链方程: dq 轴 转矩方程: dq 轴 运动方程: 分析上述方程,如果我们能够控制 id=0 那么电压...

    注:
    1:此为永磁同步控制系列文章之一,应大家的要求,关于永磁同步矢量控制的系列文章已经在主页置顶,大家可以直接去主页里面查阅,希望能给大家带来帮助,谢谢。
    2:矢量控制的六篇文章后。弱磁、MTPA、位置控制系列讲解已经补充,也放在主页了,请大家查阅。
    3: 恰饭一下,也做了一套较为详细教程放在置顶了,内含基本双闭环、MTPA、弱磁、三闭环、模糊PI等基本控制优化策略,也将滑模,MRAS等无速度控制课题整理完成,请大家查看_
    **
    4、文章对应资料附件放在了文章末尾

    1 电流内环调节器设计

    矢量控制系统的电流环是对 iq进行控制,控制的是定子电流,进而控制电机转矩。
    电流内环的作用是在电机启动过程中能够以最大电流启动,同时在外部扰动是能够快速恢复,加快动态跟踪响应速度,提高系统的稳定性。
    这里写图片描述
    上图为电流内环的流程图,电流内环的输入为电流信号的误差值,输出为参考电压,控制电动机转矩。第一个环节是PI调节器,第二个环节是延迟环节,第三个环节是PWM环节。
    其中电机传递函数可通过近似处理为:
    这里写图片描述
    在开关频率为10KHZ时,由于开关频率较高,就可以把延迟环节和PWM环节合并处理,记 td = Ts ,并将 Kpwm看成 1 来处理,可得以下流程图:
    这里写图片描述
    对以上流程图分析,将电流环按照典型的 I 型系统来整定。
    则开环传递函数:
    在这里插入图片描述

    若使得 taoi = Lq / R 可以得到 整定后开环传函:
    在这里插入图片描述

    与典型一型环节对比,(实际典型一型环节是一个二阶系统)
    ![这里写图片描述](http://latex.codecogs.com/png.latex?%5Clarge%20G%28s%29%20=%20%5Cfrac%7BK%7D%7Bs%28Ts+1%29%7得到
    可对K和T进行求解,
    在这里插入图片描述
    一阶系统按 KT = 0.5 计算得出

    这里写图片描述

    这里写图片描述

    2 转速外环调节器的设计

    转速外环设计合理的话,可以减少扰动对系统的影响、减小转速波动,使得系统工作在稳定状态。
    这里写图片描述
    在研究转速外环的时候,将电流环视为一节环节:
    这里写图片描述
    由二阶系统自身性能,在阻尼比为0.707时性能最佳,即可推:
    这里写图片描述
    同电流环,将延时环节与简化的电流环合并处理得
    这里写图片描述
    流程图进一步简化为:
    这里写图片描述
    将转速环按二阶典型环节整定,
    设转速环 PI 调节器为:
    这里写图片描述
    可得一下开环传函:
    这里写图片描述
    整理后得:
    这里写图片描述
    按照典型的二型系统的参数关系,应有
    这里写图片描述
    这里写图片描述
    由典型二阶系统整定理论得,h = 5 时 系统性能最佳。
    经过整理得到:
    这里写图片描述
    即可得PI调节器参数为
    这里写图片描述

    这里写图片描述
    由于传函很多细节部分还有得可讲,也值得从自动控制原理的角度去探究PI参数的整定,特整理了一系列专门针对PI参数整定的文章,这一系列文章可以帮助大家从头到尾理解PI参数到底从何而来,也传函框图中每个环节的由来,应该是值得大家阅读的文章。
    需要文章资料与仿真模型的同学请博客下评论留一下邮箱,看到就会发过去。

    参数整定以及自动控制原理系列文章:

    如何用matlab画bode图——自动控制原理基础补充(一)

    一阶惯性环节的性能分析——自动控制原理基础补充(二)

    二阶系统的性能分析(开环相幅和阶跃响应)——自动控制原理基础补充(三

    转速环PI参数整定详解(一)——电机传递函数的来源

    转速环PI参数整定详解(二)——转速环各个环节传递函数的来源

    转速环PI参数整定详解(三)——转速环开环传函特性及其整定策略

    整理不易,希望大家帮忙点个赞呀,谢谢啦~_

    系列文章链接:

    永磁同步电机矢量控制到无速度传感器控制学习教程(PMSM)
    永磁同步电机矢量控制(一)——数学模型
    永磁同步电机矢量控制(二)——控制原理与坐标变换推导
    永磁同步电机矢量控制(四)——simulink仿真搭建
    永磁同步电机矢量控制(五)——波形记录及其分析
    永磁同步电机矢量控制(六)——MTPA最大转矩电流比控制
    永磁同步电机矢量控制(七)——基于id=0的矢量控制的动态解耦策略
    永磁同步电机矢量控制(八)——弱磁控制(超前角弱磁)
    永磁同步电机矢量控制(九)——三闭环位置控制系统
    永磁同步电机矢量控制(十)——PMSM最优效率(最小损耗)控制策略

    展开全文
  • 天线基础与HFSS天线设计流程

    万次阅读 多人点赞 2019-04-28 15:10:10
    1.2 天线的性能参数 1. 方向图 2. 辐射强度 3.方向性系数 4. 效率 5. 增益 6. 输入阻抗 7. 天线的极化 HFSS天线设计流程 2.1 HFSS天线设计流程概述 1.设置求解类型 2.创建天线的结构模型 3.设置边界...
  • 本节主要解决5个问题: 1.动态系统建模要素 2.开放式动态系统建模 3.动态系统分类 4.建立方程模型 5.Simulink建模提示
  • 结合Matlab编程进行动态仿真分析,根据活塞下腔容积、工作载荷、螺旋槽断面尺寸和活塞半径4个主要参数变化获得的动态响应曲线,讨论了数字缸的结构参数和工作参数动态特性的影响,进行了修正系统设计;结果表明,该结构...
  • MyBatis面试题(2020最新版)

    万次阅读 多人点赞 2019-09-24 16:40:33
    整理好的MyBatis面试题库,史上最全的MyBatis面试题,...MyBatis 避免了几乎所有的 JDBC 代码和手动设置参数以及获取结果集。MyBatis 可以使用简单的 XML 或注解来配置和映射原生类型、接口和 Java 的 POJO(Plai...
  • MATLAB数字信号处理(2)LFM脉冲雷达回波处理仿真

    万次阅读 多人点赞 2019-03-13 15:41:24
    将上学期的“气象雷达原理与系统”课程报告放到blog上。 摘要 线性调频(LFM)信号是应用...本设计实现了对线性调频(LFM)脉冲压缩雷达的工作原理仿真,在MATLAB 平台中模拟一个叠加的线性调频回波信号,对该...
  • 下面三幅图分别是滤波前的时频图,滤波器的滤波特性曲线图和滤波后的时频图,通过图可以看出成功留下了400Hz的高频成分而把不要的低频成分100Hz去除了。 3.带通滤波器 function y=bandp(x,f1,f3,fsl,fsh,rp,rs,Fs) %...
  • 带通滤波器的设计

    万次阅读 多人点赞 2019-08-09 12:17:18
    一、滤波器: ...对于并联谐振电路可知,L越小C越大,衰减特性越陡。 谐振电路的谐振频率:     串联谐振与并联谐振的组合电路: 五、介质滤波器: 六、微带线制作滤波器:
  • 入门学习Linux常用必会60个命令实例详解doc/txt

    千次下载 热门讨论 2011-06-09 00:08:45
    所以这个选项当然没有时间参数,但是可以输入一个用来解释的讯息,而这信息将会送到每位使用者。 -F:在重启计算机时强迫fsck。 -time:设定关机前的时间。 -m: 将系统改为单用户模式。 -i:关机时显示系统...
  • 2、动态特性 设置参数 α=100,a=50,β=0.5,μ=1,ωsw\alpha=100,a=50,\beta=0.5, \mu=1,\omega_{sw}α=100,a=50,β=0.5,μ=1,ωsw​取不同值时的输出曲线: 设置参数 α=100,a=50,ωsw=2π,μ=1,β\alpha=100,a=50,...
  • 滤波器主要参数特性

    万次阅读 2017-12-29 11:37:09
    当滤波器幅频特性满足设计要求时,为保证输出信号失真度不超过允许范围,对其相频特性∮(w)也应提出一定要求。在滤波器设计中,常用群时延函数d∮(w)/dw*价信号经滤波后相位失真程度。群时延函数d∮(w)/dw越接近常数...
  • 针对如何确定隔振器低频动态力学特性的问题,设计了一种利用材料试验机进行加载测试以获,取其力学特性参数的试验方案。依据谐波振动理论,推导了待测参数表达式,进而提出了一种隔振器加载时力学特性参数测试新方法。以...
  • 电压电流双环控制PI参数计算01

    千次阅读 多人点赞 2020-03-07 18:38:28
    外环截止频率 2 开环截止频率和闭环截止频率 开环截止频率也称为剪切频率,是开环幅频特性中,幅频特性曲线穿越0dB线的频率,记为ωc;闭环截止频率也称为带宽频率,是指当闭环幅频特性下降到频率为零时的分贝值以下...
  • 基于对化工生产过程参数的精确控制的目的,采用建立机理动态数学模型的方法,该控制方法不仅能给出过程稳态特性设计,同时能给出过程动态特性设计,必将使设计工作进入到一个新的更高的水平。通过水位、温度控制...
  • Java面向对象

    千次阅读 多人点赞 2019-05-08 23:33:50
    Java面向对象1.Java面向对象是什么2.面向过程和面向对象的优缺点1.... OOP(面向对象编程)、OOD(面向对象设计)、OOA(面向对象的分析) Object Oriented Programming     &nb...
  • 半导体二极管正向特性参数测试电路如图13.2-1所示。表13.2-1是正向测试的数据,从仿真数据可以看出:二极管电阻值不是固定值,当二极管两端正向电压小,处于“死区”,正向电阻很大、正向电流很小,当二极管两端正向...
  • 天线设计相关性能参数

    万次阅读 多人点赞 2018-10-06 09:50:46
    天线的指标参数分为两部分:1.电路参数,天线高效率辐射的保证,天线应用的必要条件 2.辐射参数,天线应用的本质,天线应用的充分条件 下面将每部分的参数做介绍(终端小天线大多数时候较强调电路参数)。电路参数:...
  • MOS管的主要参数与重要特性

    千次阅读 2019-05-23 14:23:04
    双极性晶体管:NPN和PNP管; 单极性晶体管:场效应管(MOSFET和JFET); MOS管相对三极管具有速度快、输入阻抗高、...静态特性参数 BVDSS:漏源击穿电压,为正温度系数。 BVGS:栅源击穿电压 VGS(th):阈值电...
  • PID参数整定一些总结

    千次阅读 2020-06-12 17:34:01
    PID参数整定法总结 PID控制规律 2、PID传递函数 3、各环节的作用 比例环节作用: 系统一旦出现偏差,比例环节立即产生调节作用以减小系统偏差,比例作用大 ,可以加快调节,减少误差,但是过大的比例,使系统的...
  • 自适应控制基本思想

    万次阅读 多人点赞 2017-09-22 10:18:20
    自适应控制所讨论的对象,一般是指对象的结构已知,仅仅是参数未知,而且采用的控制方法仍是基于数学模型的方法 但实践中我们还会遇到结构和参数都未知的对象,比如一些运行机理特别复杂,目前尚未被人们充分理解的...
  • JAVA上百实例源码以及开源项目

    千次下载 热门讨论 2016-01-03 17:37:40
    例如,容易实现协议的设计。 Java EJB中有、无状态SessionBean的两个例子 两个例子,无状态SessionBean可会话Bean必须实现SessionBean,获取系统属性,初始化JNDI,取得Home对象的引用,创建EJB对象,计算利息等;...
  • 关注、星标公众号,不错过精彩内容直接来源:21ic电子网之前给大家分享过PID基础理论的文章:重温经典PID算法PID原理和参数调试今天进一步分享一些PID相关细节内容。在过程控制中,按...
  • 研究、开发了真空断路器机械特性测试系统,克服了以往的一些应用难题,能实现包括机械、电气、触头温度等多种断路器运行特性的综合监测;该系统采用上下位机结构,下位机是一个基于TI公司的VC33数字信号处理器的数据...
  • 基于MATLAB的语音信号处理

    万次阅读 多人点赞 2018-07-15 01:21:20
     语音信号分析是语音信号处理的前提和基础,只有分析出可表示语音信号本质特征的参数,才有可能利用这些参数进行高效的语音通信、语音合成和语音识别等处理。而且,语音合成的音质好坏,语音识别率的高低,也都取决...
  • Python动态特性

    千次阅读 2015-11-22 02:18:19
    长期习惯了C/C++系的静态语言后,切换到Python中往往仍习惯使用静态办法解决问题,而不能充分利用Python强大的动态特性。大多数时候,这使得代码变得不必要地长且难以理解。 希望在进阶Python的过程中,逐步掌握...
  • 利用matlab来设计FIR滤波器参数

    千次阅读 2019-01-23 10:02:05
    MATLAB自带fdatool,可以用来生成参数 输入信号:采样率为48KHz的单声道立体声信号(Fs) 通带频率:15KHz(Fpass) 阻带频率:16KHz(Fstop) 通带平坦度:0.1dB(Apass) 阻带衰减:80dB(Astop) 滤波器...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 505,367
精华内容 202,146
关键字:

动态特性参数设计