自动控制原理PID参数整定的Matlab实现 |
您所在的位置:网站首页 › pid调节器参数整定方法 › 自动控制原理PID参数整定的Matlab实现 |
以一道题为例介绍调节PID控制器系数的方法,有:试凑法(Trial-and-Error Method)、齐格勒-尼科尔斯校正规则(Ziegler and Nichols First Method、Ziegler and Nichols Second Method),并在Matlab中进行仿真实验。 这是出自《线性控制系统工程/Linear Control Systems Engineering》的21.9题,要求设置PID控制器的参数。 无控制器时的响应先分析没有PID控制器影响下的系统的阶跃响应,此时系统的开环传函为: G H ( s ) = 1 ( s + 3 ) 3 = 1 s 3 + 9 s 2 + 27 s + 27 GH(s)=\frac{1}{(s+3)^3} =\frac{1}{s^3+9s^2+27s+27} GH(s)=(s+3)31=s3+9s2+27s+271 采用以下的Matlab代码,得到系统阶跃响应曲线: num = 1; den = [1 9 27 27]; G = tf(num,den); sys = feedback(G,1); [y,t] = step(sys); plot(t,y); hold on legend('Step Response'); ylabel('Amplitude');xlabel('Time(s)'); grid;简单分析:在阶跃信号作用下,该系统大概经过3.1秒就达到稳态值,没有超调量,稳态值约为0.036,稳态误差为: e s s = 1 1 + K p = 1 1 + lim s → 0 G ( s ) = 27 28 ≈ 0.964 \mathrm{e}_{\mathrm{ss}}=\frac{1}{1+\mathrm{K}_{\mathrm{p}}}=\frac{1}{1+\lim _{\mathrm{s} \rightarrow 0} \mathrm{G}(\mathrm{s})}=\frac{27}{28} \approx 0.964 ess=1+Kp1=1+lims→0G(s)1=2827≈0.964 试凑法(Trial-and-Error Method)如果将PID控制的传递函数用Matlab的行向量表示,由于将 直接代入会产生无穷大使得程序出错,所以先对系统框图进行化简,求出其闭环传递函数: C R = 1 1 + K p ( 1 + 1 T i s + T d s ) ( s + 1 ) 3 = K p ⋅ T d ⋅ T i ⋅ s 2 + K p ⋅ T i ⋅ s + K p T i ⋅ s 4 + 9 T i ⋅ s 3 + ( 27 + K p ⋅ T d ) ⋅ T i ⋅ s 2 + ( 27 + K p ) ⋅ T i ⋅ s + K p \begin{aligned} \frac{C}{R} & =\frac{1}{1+K_{p}\left(1+\frac{1}{T_{i}s}+\mathrm{T}_{\mathrm{d}} \mathrm{s}\right)(\mathrm{s}+1)^{3}} \\ & =\frac{\mathrm{K}_{\mathrm{p}} \cdot \mathrm{T}_{\mathrm{d}} \cdot \mathrm{T}_{\mathrm{i}} \cdot \mathrm{s}^{2}+\mathrm{K}_{\mathrm{p}} \cdot \mathrm{T}_{\mathrm{i}} \cdot \mathrm{s}+\mathrm{K}_{\mathrm{p}}}{\mathrm{T}_{\mathrm{i}} \cdot \mathrm{s}^{4}+9 \mathrm{~T}_{\mathrm{i}} \cdot \mathrm{s}^{3}+\left(27+\mathrm{K}_{\mathrm{p}} \cdot \mathrm{T}_{\mathrm{d}}\right) \cdot \mathrm{T}_{\mathrm{i}} \cdot \mathrm{s}^{2}+\left(27+\mathrm{K}_{\mathrm{p}}\right) \cdot \mathrm{T}_{\mathrm{i}} \cdot \mathrm{s}+\mathrm{K}_{\mathrm{p}}} \end{aligned} RC=1+Kp(1+Tis1+Tds)(s+1)31=Ti⋅s4+9 Ti⋅s3+(27+Kp⋅Td)⋅Ti⋅s2+(27+Kp)⋅Ti⋅s+KpKp⋅Td⋅Ti⋅s2+Kp⋅Ti⋅s+Kp 首先令 T i = 0 T_{i}=0 Ti=0、 T d = 0 T_{d}=0 Td=0,调节 K p K_{p} Kp以满足期望的动态特性,但是此时调节 K p K_{p} Kp会发现:当 K p > 0 K_{p}>0 Kp>0时,阶跃响应均为一条直线,且幅值为1,如图所示。 因为此时闭环传函 G H ( s ) = K p K p = 1 GH(s)=\frac{K_{p}}{K_{p}}=1 GH(s)=KpKp=1,有 c ( t ) = L − 1 [ G H ( s ) ⋅ R ] = L − 1 [ 1 s ] = 1 c(t)=\mathscr{L}^{-1}\left[ GH(s)\cdot {R} \right] =\mathscr{L}^{-1}\left[ \frac{1}{s} \right] =1 c(t)=L−1[GH(s)⋅R]=L−1[s1]=1 所以需要增加一个变量 T i T_{i} Ti来共同调节。设 K p = 100 K_{p}=100 Kp=100、 T i = 1 T_{i}=1 Ti=1、 T d = 0 T_{d}=0 Td=0可得如下所示的阶跃响应。 代码如下 Kp = 100; Ti = 1; Td = 0; sys = tf([Kp*Td*Ti Kp*Ti Kp],[Ti 9*Ti (27+Kp*Td)*Ti (27+Kp)*Ti Kp]); [y,t] = step(sys); %得到阶跃响应达到稳定前的幅值和时间的行向量y,t [Y,T] = max(y); %找出阶跃响应达到稳定前的幅值y中最大值Y和所在时间T %求出阶跃响应达到稳定的时间SteadyStateTime,即求t的长度 SteadyStateTime = length(t); %求出阶跃响应达到稳定的幅值SteadyStateOutput,即t= SteadyStateTime时,y的值 SteadyStateOutput = y(SteadyStateTime); [pks,locs] = findpeaks(y,t);%找出y和t中的峰值pks和达到峰值时间locs figure(1);plot(t,y);hold on if isempty(pks) == 0%如果找到了pks和 locs就进入if条件 PO1 = 100*(pks(1)-SteadyStateOutput)/SteadyStateOutput;%计算第一个超调量 PO2 = 100*(pks(2)-SteadyStateOutput)/SteadyStateOutput;%计算第二个超调量 plot(locs(1),pks(1),'o');hold on;%在图中用“o”标出第一个超调量 plot(locs(2),pks(2),'*');hold on;%在图中用“*”标出第二个超调量 POP = PO2/PO1;%计算两个超调量的比值 POtext = strcat("第二次超调量是第一次超调量的 ",num2str(POP)," 倍"); text(locs(2)*0.8,pks(1)*0.5,POtext); legend(' Step Response',' First Peak',' Second Peak'); else legend(' Step Response'); end ylabel('Amplitude');xlabel('Time(s)'); grid;由于加入微分环节会增大系统的阻尼比,能用来改善系统的动态特性,一定程度上可以用来减少超调量、缩短调节时间。所以在加入微分环节前,虽然目前 K p = 100 K_{p}=100 Kp=100、 T i = 1 T_{i}=1 Ti=1、 T d = 0 T_{d}=0 Td=0时的阶跃曲线,不满足第二次超调量是第一次超调量的1/4的最优性能的条件,但在加入微分环节后可以改善系统的阻尼比,可以使得上述参数下的系统阶跃响应往“好”的响应发展。 令 K p = 100 K_{p}=100 Kp=100、 T i = 1 T_{i}=1 Ti=1、 T d = 0.1 T_{d}=0.1 Td=0.1,此时系统的阶跃响应如图所示。 此时系统的阶跃响应仍不满足第二次超调量是第一次超调量的1/4的最优性能的条件,但可以通过增大比例环节的系数 K p K_{p} Kp来增大系统的超调量来满足这个条件。令 K p = 120 K_{p}=120 Kp=120、 T i = 1 T_{i}=1 Ti=1、 T d = 0.1 T_{d}=0.1 Td=0.1,此时系统的阶跃响应如图所示。 此时第二次超调量是第一次的0.249倍,基本上符合了最优性能的指标,但不难看出系统达到稳态的时间约为5秒,试凑出来的PID控制器的系数仍有不足。 完整的Matlab代码如下: Kp = 120; Ti = 1; Td = 0.1; sys = tf([Kp*Td*Ti Kp*Ti Kp],[Ti 9*Ti (27+Kp*Td)*Ti (27+Kp)*Ti Kp]); [y,t] = step(sys); %得到阶跃响应达到稳定前的幅值和时间的行向量y,t [Y,T] = max(y); %找出阶跃响应达到稳定前的幅值y中最大值Y和所在时间T %求出阶跃响应达到稳定的时间SteadyStateTime,即求t的长度 SteadyStateTime = length(t); %求出阶跃响应达到稳定的幅值SteadyStateOutput,即t= SteadyStateTime时,y的值 SteadyStateOutput = y(SteadyStateTime); [pks,locs] = findpeaks(y,t);%找出y和t中的峰值pks和达到峰值时间locs figure(1);plot(t,y);hold on if isempty(pks) == 0%如果找到了pks和 locs就进入if条件 PO1 = 100*(pks(1)-SteadyStateOutput)/SteadyStateOutput;%计算第一个超调量 PO2 = 100*(pks(2)-SteadyStateOutput)/SteadyStateOutput;%计算第二个超调量 plot(locs(1),pks(1),'o');hold on;%在图中用“o”标出第一个超调量 plot(locs(2),pks(2),'*');hold on;%在图中用“*”标出第二个超调量 POP = PO2/PO1;%计算两个超调量的比值 POtext = strcat("第二次超调量是第一次超调量的 ",num2str(POP)," 倍"); text(locs(2)*0.8,pks(1)*0.5,POtext); legend(' Step Response',' First Peak',' Second Peak'); else legend(' Step Response'); end ylabel('Amplitude');xlabel('Time(s)'); grid; 扩展:衰减曲线法Trial-and-Error Method是令 T i = 0 T_{i}=0 Ti=0、 T d = 0 T_{d}=0 Td=0,先调节 K p K_{p} Kp、再调节 T i T_{i} Ti、最后调节 T d T_{d} Td。但如果是直接从 K p K_{p} Kp、 K i K_{i} Ki、 K d K_{d} Kd来进行调节,又该如何操作呢?有一种方法叫“衰减曲线法”可以很好地帮助我们解决这个问题。(出自:《基于衰减曲线法整定PID调节器参数仿真技术研究》) 首先将 T i = ∞ T_{i}=\infty Ti=∞、 T d = 0 T_{d}=0 Td=0,逐渐增大 K p K_{p} Kp,使得系统的阶跃响应曲线出现第二次超调量是第一次超调量的0.25倍,此时令 δ s = K p \delta_{s}=K_{p} δs=Kp,计算此时两次超调量的时间间隔 T s T_{s} Ts,代入下表即可计算出最终的PID控制器的系数。 衰减率控制器比例系数 K p K_{p} Kp积分系数 K i K_{i} Ki微分系数 K d K_{d} Kd0.75P δ s \delta_{s} δs000.75PI 1.2 δ s 1.2\delta_{s} 1.2δs 2.4 δ s T s 2.4\delta_{s}T_{s} 2.4δsTs00.75PID 0.8 δ s 0.8\delta_{s} 0.8δs 2.66 δ s T s \frac{2.66\delta_{s}}{T_{s}} Ts2.66δs 0.08 δ s T s 0.08\delta_{s}T_{s} 0.08δsTs但在第一步将 T i = ∞ T_{i}=\infty Ti=∞是不符合Matlab语法规定的,所以可以改变系统的PID控制器的传递函数为: θ o θ i = K p ( 1 + 1 T i s + T d s ) = K p + K i s + K d s \frac{\theta_{o}}{\theta_{i}}=K_{p}(1+\frac{1}{T_{i}s}+T_{d}s)=K_{p}+\frac{K_{i}}{s}+K_{d} s θiθo=Kp(1+Tis1+Tds)=Kp+sKi+Kds 也就实现了直接从 K p K_{p} Kp、 K i K_{i} Ki、 K d K_{d} Kd来进行调节的方法。这里使用Simulink来进行仿真更方便。 先令 K i = 0 K_{i}=0 Ki=0、 K d = 0 K_{d}=0 Kd=0,借助Matlab的Simulink工具搭建可如下所示的模块图。 只需修改 K p K_{p} Kp的Gain使得在Scpoe模快显示的系统的阶跃响应曲线出现第二次超调量是第一次超调量的0.25倍即可,发现当 K p = 81.8 K_{p}=81.8 Kp=81.8时满足此条件,所以有 δ s = 81.8 \delta_{s}=81.8 δs=81.8。 在找出第二次超调量是第一次超调量的0.25倍的 δ s \delta_{s} δs值后,仍需求出其第一次超调量与第二次超调量的时间间隔 T s T_{s} Ts。在搭建完成此Simulink模块图后,进行一些简单的设置可以更方便地展示两次超调的时间间隔: 1、双击Step模块,调节Step time(阶跃发生时间)设置为0s(默认为1s)。 2、设置菜单栏上的SIMULATE工具中的Stop time为10,如图所示。 3、为了方便用Matlab代码来分析,需要将Scope模块的数据导出到Matlab的工作区,双击Scope模块,将Logging功能设置为如图所示,启用仿真功能后会输出一个命名为“out”的结构体到Matlab的工作区。
启用Simulink的仿真功能和运行上述Matlab代码后可得下图。 得到了第二次超调量是第一次超调量的0.25倍的
δ
s
=
81.8
\delta_{s}=81.8
δs=81.8、两次超调时间间隔
T
s
=
1.67
T_{s}=1.67
Ts=1.67。根据上表的经验公式计算其最终的PID控制器系数: 比例系数:
K
p
=
0.8
δ
s
=
0.8
∗
81.8
=
65.44
K_{p}=0.8\delta_{s}=0.8 * 81.8 =65.44
Kp=0.8δs=0.8∗81.8=65.44 积分系数:
K
i
=
2.66
δ
s
T
s
=
2.66
∗
81.8
1.67
≈
130.30
K_{i}=\frac{2.66\delta_{s}}{T_{s}}=\frac{2.66 * 81.8}{1.67} \approx 130.30
Ki=Ts2.66δs=1.672.66∗81.8≈130.30 微分系数:
K
d
=
0.08
δ
s
T
s
=
0.08
∗
81.8
∗
1.67
≈
10.93
K_{d}=0.08\delta_{s}T_{s}=0.08 * 81.8 * 1.67 \approx 10.93
Kd=0.08δsTs=0.08∗81.8∗1.67≈10.93 在采用这个方法前,需判断没有PID控制器影响下的系统开环传函的阶跃响应是否有超调,如果没有则可以使用这个方法,本文的开头以通过实验证明过这点。 通过下面代码可以得到没有PID控制器影响下的系统的阶跃响应曲线,并显示了P、L和PID的参数。 num = 1; den = [1 9 27 27]; G = tf(num,den); sys = feedback(G,1); [y,t] = step(sys); [P,idx] = max(diff(y)./diff(t));%求出响应曲线最大的斜率P及其下标idx Tangent_b = y(idx) - P * t(idx);%最大斜率P的切线Tangent的常数项 L = - Tangent_b/P;%该切线在横轴上的交点 t0 = P * max(t) + Tangent_b; plot(t,y,'b',[L max(t)],[0 t0],'r--') axis([0,max(t),0,max(y)+0.008]); legend('Step Response','Tangent'); ylabel('Amplitude');xlabel('Time(s)'); Ptext = strcat("响应曲线最大的斜率P为 ",num2str(P));text(0.8,0.015,Ptext); Ltext = strcat("最大斜率的切线在横轴上的交点L为 ",num2str(L));text(0.8,0.012,Ltext); PIDtext = strcat("Kp、Ti、Td的系数为 ",num2str(1.2/(P*L)),"、",num2str(2*L),"、",num2str(0.5*L));%代入公式,计算Kp、Ti、Td text(0.8,0.009,PIDtext); grid;
采用这种方法的系统的开环传函的阶跃响应为什么不能有超调?假设一个没有PID控制器的系统的开环传递函数为:
G
H
(
s
)
=
1
(
s
+
3
)
3
=
1
s
3
+
20
s
2
+
27
s
+
27
GH(s)=\frac{1}{(s+3)^3} =\frac{1}{s^3+20s^2+27s+27}
GH(s)=(s+3)31=s3+20s2+27s+271 它的阶跃响应曲线如图所示,有超调。 在只考虑系统的比例控制的情况下,求出系统的临界增益 K p ′ K_{p}' Kp′,使得系统响应曲线出现等幅振荡,此时系统的开环传递函数: G H ( s ) = 1 ( s + 3 ) 3 = 1 s 3 + 9 s 2 + 27 s + 27 GH(s)=\frac{1}{(s+3)^3} =\frac{1}{s^3+9s^2+27s+27} GH(s)=(s+3)31=s3+9s2+27s+271 得特征方程为 s 3 + 9 s 2 + 27 s + ( 27 + K p ′ ) = 0 \text{s}^3+9\text{s}^2+27\text{s}+\left( 27+\text{K}_{\text{p}}' \right) =0 s3+9s2+27s+(27+Kp′)=0,可得劳斯表为: s 3 s^3 s3127 s 2 s^2 s2927 + K p ′ K_{p}' Kp′ s 1 s^1 s1 24 − K p ′ 9 24-\frac{K_{p}'}{9} 24−9Kp′0 s 0 s^0 s0 27 + K p ′ 9 27+\frac{K_{p}'}{9} 27+9Kp′0可知系统的临界稳定条件为:
24
−
K
p
′
9
=
0
24-\frac{K_{p}'}{9}=0
24−9Kp′=0,即
K
p
′
=
216
K_{p}'=216
Kp′=216。同时在计算劳斯表时不难发现:这种方法只适用于三阶及其以上的系统,二阶系统是算不出临界增益
K
p
′
K_{p}'
Kp′的。此时系统的阶跃响应曲线如图所示,系统处于等幅振荡,通过Matlab易求出振荡周期
P
′
=
1.21
s
P'=1.21s
P′=1.21s。 假设一个二阶系统的特征方程为: s 2 + α x + β = 0 s^2+\alpha x+\beta=0 s2+αx+β=0 用第二方法整定PID参数时,有 s 2 + α x + ( β + K p ′ ) = 0 s^2+\alpha x+(\beta + K_{p}')=0 s2+αx+(β+Kp′)=0 有劳斯表: s 2 s^2 s21 β + K p ′ \beta + K_{p}' β+Kp′ s 1 s^1 s1 α \alpha α0 s 0 s^0 s0 α ( β + K p ′ ) α \frac{\alpha(\beta + K_{p}')}{\alpha} αα(β+Kp′)劳斯表第一列均大于0,求不出临界稳定时 K p ′ K_{p}' Kp′ PID Tuner PID Tuner是Maltab的Simulink工具中PID控制器模块提供的一个功能,它能够自动地调节PID控制器的参数,使得系统的响应速度和动态特性达到最大平衡,也可以人为地调节系统的响应速度和动态特性。需要在Simulink中搭建一个如图所示的模块图。 |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |