PR控制器以及准PR控制器原理分析

  工程近似PR控制器传递函数为： G P R ( s ) = ≈ 1 2 [ G P I ( s + j w 0 ) + G P I ( s − j w 0 ) ] = K p + K i s s 2 + w o 2 G_{PR}(s)=\approx \frac{1}{2}[G_{PI}(s+jw_0)+G_{PI}(s-jw_0)]=Kp+\frac{K_is}{s^2+w_o^2} GPR​(s)=≈21​[GPI​(s+jw0​)+GPI​(s−jw0​)]=Kp+s2+wo2​Ki​s​   根据 L { e − a t } = 1 s + a L\{ e^{-at} \}=\frac{1}{s+a} L{e−at}=s+a1​有： F ( s ) = L { c o s ( w t ) } = L { ∫ 0 ∞ e − s t c o s ( w t ) d t } = L { 1 2 ∫ 0 ∞ e − s t ( e j w t + e − j w t ) } ⇒ s s 2 + w 2 F(s)=L\{cos(wt)\}=L\{ \int_{0}^{\infty} {e^{-st}cos(wt)dt}\}=L\{ \frac{1}{2} \int_{0}^{\infty}{e^{-st}(e^{jwt}+e^{-jwt})} \} \Rightarrow \frac{s}{s^2+w^2} F(s)=L{cos(wt)}=L{∫0∞​e−stcos(wt)dt}=L{21​∫0∞​e−st(ejwt+e−jwt)}⇒s2+w2s​   定义： G P R ( s ) = K p + K r s s 2 + w 0 2 G_{PR}(s)=K_p+\frac{K_rs}{s^2+w_0^2} GPR​(s)=Kp​+s2+w02​Kr​s​显然谐振环节是PR控制器的核心，其中， K p 、 K r K_p、K_r Kp​、Kr​分别为比例增益系数和谐振增益系数， w 0 w_0 w0​为谐振频率。   PR控制器的增益函数为： ∣ G P R ( s ) I s = j w 0 ∣ = K p 2 + K r w 0 − w 0 2 + w 0 2 |G_{PR}(s)I_{s=jw_0}|=\sqrt{K_p^2+ \frac{K_rw_0}{-{w_0}^2+w_0^2}} ∣GPR​(s)Is=jw0​​∣=Kp2​+−w0​2+w02​Kr​w0​​ ​   PR控制器在 w 0 w_0 w0​处的增益接近于无限大，在其他频率下增益低，能够有效地抑制扰动信号。可以把PR看作带宽极窄的二阶带通滤波器。

  下图为PR控制器的波特图：(注：横坐标单位是rad/s，可以右键点击“属性”–>“单位”–>“频率”–>“HZ”)   附matlab绘图代码：

  理想的PR控制器是完全可以实现对应频率的交流量实现无静差跟踪的，但是在谐振频率附近的频段带宽过于狭窄，而在 w 0 w_0 w0​处的增益过高，会使得系统的稳定性不够，当交流信号发生些许偏移时，PR控制器就无法精准工作在预设频率上了，虽然可以通过调节 K r K_r Kr​增大带宽，但是会使得增益变化和相位变化明显增大，会造成系统不稳定。由于PR控制器对于电网参数过于敏感，所以通常不在实际中运用。

  为了提高PR控制器抵抗网侧频率干扰的能力，对PR控制器进行改进，改进后的传递函数如下所示： G P R ( s ) = K p + 2 K r w c s s 2 + 2 w c s + w 0 2 G_{PR}(s)=K_p+\frac{2K_rw_cs}{s^2+2w_cs+w_0^2} GPR​(s)=Kp​+s2+2wc​s+w02​2Kr​wc​s​   其中 w c w_c wc​为截止频率，代表控制器跟踪参考信号的响应速度   准PR控制器的增益函数： ∣ G P R ( s ) ∣ s = j w 0 = K p + K r |G_{PR}(s)|_{s=jw_0}=K_p+K_r ∣GPR​(s)∣s=jw0​​=Kp​+Kr​   根据准PR增益函数可知，当输入信号频率为 w 0 w_0 w0​时，增益为 ( K p + K r ) (K_p+K_r) (Kp​+Kr​)，不再像PR控制器那样增益无穷大。   下图为准PR控制器的波特图：   附matlab绘图代码：

  通过准PR控制器的Bode图可知增益幅度符合传递函数所述，同时调节谐振增益系数可以增大谐振频率附近的频段带宽。

  使用控制变量，分别对与不同的Kp、Kr、wc进行比较，熟悉不同参数对于控制器带来的影响，如下图所示：

  附matlab绘图代码：

在进行理论分析时，Matlab实现离散化很方便。   当 K p = 10 、 K r = 100 、 w c = 0.5 ∗ 2 ∗ π 、 w o = 100 π K_p=10、K_r=100、wc=0.5*2*\pi、wo=100\pi Kp​=10、Kr​=100、wc=0.5∗2∗π、wo=100π时，连续时间模型为：   当采样时间 T s = 1 0 − 6 Ts=10^{-6} Ts=10−6时，其他参数不变，离散时间模型为：

  附matlab转换代码：(注：表达式后面不要加“ ; ”，静态检查的警告忽略)

  在DSP或单片机中对于仅改变 K p 、 K r 、 w c 、 w o K_p、K_r、wc、wo Kp​、Kr​、wc、wo离散化传递函数可过程借助编程得到，而不需要借助三方软件计算得到离散化传递函数，故而具备离散化的计算推导能力是很有必要的。   使用Tustin变换（大多数DSP厂家算法库的选择）， s = 2 T s z − 1 z + 1 s=\frac{2}{T_s}\frac{z-1}{z+1} s=Ts​2​z+1z−1​，带入PR传递函数，便可以得到PR控制器的差分方程，再根据差分方程得到离散域表达式（二者只是形式不一样），根据Z域表达式进行代码实现。 差分方程： Y ( z ) X z = a 0 + a 1 z − 1 + a 2 z − 2 + . . . + a k z − k b 0 + b 1 z − 1 + b 2 z − 2 + . . . + b k z − k \frac{Y(z)}{X{z}}=\frac{a_0+a_1z^{-1}+a_2z^{-2}+...+a_kz^{-k}}{b_0+b_1z^{-1}+b_2z^{-2}+...+b_kz^{-k}} XzY(z)​=b0​+b1​z−1+b2​z−2+...+bk​z−ka0​+a1​z−1+a2​z−2+...+ak​z−k​   Z域表达式： b 0 y [ n ] + b 1 y [ n − 1 ] + b 2 y [ n − 2 ] + . . . + b k y [ n − k ] = a 0 x [ n ] + a 1 x [ n − 1 ] + a 2 x [ n − 2 ] + . . . + a k x [ n − k ] b_0y[n]+b_1y[n-1]+b_2y[n-2]+...+b_ky[n-k]=a_0x[n]+a_1x[n-1]+a_2x[n-2]+...+a_kx[n-k] b0​y[n]+b1​y[n−1]+b2​y[n−2]+...+bk​y[n−k]=a0​x[n]+a1​x[n−1]+a2​x[n−2]+...+ak​x[n−k]   PR控制器离散化推导过程： P R ( s ) = K p + 2 K r w c s s 2 + 2 w c s + w 0 2 = s p K + 2 w c s K p + w 0 2 K p + 2 K r w c s s 2 + 2 w c s + w 0 2 = s 2 K p + 2 w c s ( K p + K r ) + w 0 2 K p s 2 + 2 w c s + w 0 2 \begin{aligned} PR(s)&=K_p+\frac{2K_rwcs}{s^2+2wcs+w_0^2}\\ &=\frac{s^K_p+2wcsK_p+w_0^2K_p+2K_rwcs}{s^2+2wcs+w0^2}\\ &=\frac{s^2K_p+2wcs(K_p+K_r)+w_0^2K_p}{s^2+2wcs+w_0^2}\\ \end{aligned} PR(s)​=Kp​+s2+2wcs+w02​2Kr​wcs​=s2+2wcs+w02spK​+2wcsKp​+w02​Kp​+2Kr​wcs​=s2+2wcs+w02​s2Kp​+2wcs(Kp​+Kr​)+w02​Kp​​​   将 s = 2 T s z − 1 z + 1 s=\frac{2}{T_s}\frac{z-1}{z+1} s=Ts​2​z+1z−1​带入上式，得差分方程： P R ( z ) = Y ( z ) X ( z ) = K p ( 2 T s z − 1 z + 1 ) 2 + ( 2 w c K p + 2 w c K r ) ( 2 T s z − 1 z + 1 ) + K p w 0 2 ( 2 T s z − 1 z + 1 ) 2 + 2 w c ( 2 T s z − 1 z + 1 ) + w 0 2 = 4 K p T s 2 ( z 2 − 2 z + 1 ) + 4 w c T s ( K p + K r ) ( z 2 − 1 ) + K p w 0 2 ( z 2 + 2 z + 1 ) 4 T 2 ( z 2 − 2 z + 1 ) + 4 w c T s ( z 2 − 1 ) + w 0 2 ( z 2 + 2 z + 1 ) = ( 4 K p T s 2 + 4 w c T s + K p w 0 2 ) z 2 + ( − K p T s 2 + 2 K p w 0 2 ) z + [ 4 K p T s 2 − 4 w c T s ( K p + K r ) + K p w 0 2 ] ( 4 T s 2 + 4 w c T s + w 0 2 ) z 2 + ( − 8 T s 2 + 2 w 0 2 ) z + ( 4 T s 2 − 4 w c T s + w 0 2 ) = [ 4 K p T s 2 + 4 w c T s ( K p + K r ) + K p w 0 2 ] + ( − 8 K p T s 2 + 2 K p w 0 2 ) z − 1 + [ 4 K p T s 2 − 4 w c T s ( K p + K r ) + K p w 0 2 ] z − 2 ( 4 T s 2 + 4 w c T s + w 0 2 ) + ( − 8 T s 2 + 2 w 0 2 ) z − 1 + ( 4 T s 2 − 4 w c T s + w 0 2 ) z − 2 \begin{aligned} PR(z)=\frac{Y(z)}{X(z)}&=\frac{K_p(\frac{2}{T_s}\frac{z-1}{z+1})^2+(2wcK_p+2wcK_r)(\frac{2}{T_s}\frac{z-1}{z+1})+K_pw_0^2}{(\frac{2}{T_s} \frac{z-1}{z+1})^2+2wc(\frac{2}{T_s}\frac{z-1}{z+1})+w_0^2}\\ &=\frac{\frac{4K_p}{T_s^2}(z^2-2z+1)+\frac{4wc}{T_s}(K_p+K_r)(z^2-1)+K_pw_0^2(z^2+2z+1)}{\frac{4}{T^2}(z^2-2z+1)+\frac{4wc}{T_s}(z^2-1)+w_0^2(z^2+2z+1)}\\ &=\frac{(\frac{4K_p}{T_s^2}+\frac{4wc}{T_s}+K_pw_0^2)z^2+(-\frac{K_p}{T_s^2}+2K_pw_0^2)z+[\frac{4K_p}{T_s^2}-\frac{4wc}{T_s}(K_p+K_r)+K_pw_0^2]}{(\frac{4}{T_s^2}+\frac{4wc}{T_s}+w_0^2)z^2+(-\frac{8}{T_s^2}+2w_0^2)z+(\frac{4}{T_s^2}-\frac{4wc}{T_s}+w_0^2)}\\ &=\frac{[\frac{4K_p}{T_s^2}+\frac{4wc}{T_s}(K_p+K_r)+K_pw_0^2]+(-\frac{8K_p}{T_s^2}+2K_pw_0^2)z^{-1}+[\frac{4K_p}{T_s^2}-\frac{4wc}{T_s}(K_p+K_r)+K_pw_0^2]z^{-2}}{(\frac{4}{T_s^2}+\frac{4wc}{T_s}+w_0^2)+(-\frac{8}{T_s^2}+2w_0^2)z^{-1}+(\frac{4}{T_s^2}-\frac{4wc}{T_s}+w_0^2)z^{-2}}\\ \end{aligned} PR(z)=X(z)Y(z)​​=(Ts​2​z+1z−1​)2+2wc(Ts​2​z+1z−1​)+w02​Kp​(Ts​2​z+1z−1​)2+(2wcKp​+2wcKr​)(Ts​2​z+1z−1​)+Kp​w02​​=T24​(z2−2z+1)+Ts​4wc​(z2−1)+w02​(z2+2z+1)Ts2​4Kp​​(z2−2z+1)+Ts​4wc​(Kp​+Kr​)(z2−1)+Kp​w02​(z2+2z+1)​=(Ts2​4​+Ts​4wc​+w02​)z2+(−Ts2​8​+2w02​)z+(Ts2​4​−Ts​4wc​+w02​)(Ts2​4Kp​​+Ts​4wc​+Kp​w02​)z2+(−Ts2​Kp​​+2Kp​w02​)z+[Ts2​4Kp​​−Ts​4wc​(Kp​+Kr​)+Kp​w02​]​=(Ts2​4​+Ts​4wc​+w02​)+(−Ts2​8​+2w02​)z−1+(Ts2​4​−Ts​4wc​+w02​)z−2[Ts2​4Kp​​+Ts​4wc​(Kp​+Kr​)+Kp​w02​]+(−Ts2​8Kp​​+2Kp​w02​)z−1+[Ts2​4Kp​​−Ts​4wc​(Kp​+Kr​)+Kp​w02​]z−2​​   Z域表达式： b 0 y [ n ] + b 1 y [ n − 1 ] + b 2 y [ n − 2 ] = a 0 x [ n ] + a 1 x [ n − 1 ] + a 2 x [ n − 2 ] b_0y[n]+b_1y[n-1]+b_2y[n-2]=a_0x[n]+a_1x[n-1]+a_2x[n-2] b0​y[n]+b1​y[n−1]+b2​y[n−2]=a0​x[n]+a1​x[n−1]+a2​x[n−2] { a 0 = ( 4 K p T s 2 + 4 w c T s ( K p + K r ) + K p w 0 2 ) a 1 = ( − 8 K p T s 2 + 2 K p w 0 2 ) a 2 = [ 4 K p T s 2 − 4 w c T s ( K p + K r ) + K p w 0 2 ] { b 0 = ( 4 T s 2 + 4 w c T s + w 0 2 ) b 1 = ( − 8 T s 2 + 2 w 0 2 ) b 2 = ( 4 T s 2 − 4 w c T s + w 0 2 ) \begin{cases} a_0=(\frac{4K_p}{T_s^2}+\frac{4wc}{T_s}(K_p+K_r)+K_pw_0^2)\\ a_1=(-\frac{8K_p}{T_s^2}+2K_pw_0^2)\\ a_2=[\frac{4K_p}{T_s^2}-\frac{4wc}{T_s}(K_p+K_r)+K_pw_0^2]\\ \end{cases} \begin{cases} b_0=(\frac{4}{T_s^2}+\frac{4wc}{T_s}+w_0^2)\\ b_1=(-\frac{8}{T_s^2}+2w_0^2)\\ b_2=(\frac{4}{T_s^2}-\frac{4wc}{T_s}+w_0^2)\\ \end{cases} ⎩ ⎨ ⎧​a0​=(Ts2​4Kp​​+Ts​4wc​(Kp​+Kr​)+Kp​w02​)a1​=(−Ts2​8Kp​​+2Kp​w02​)a2​=[Ts2​4Kp​​−Ts​4wc​(Kp​+Kr​)+Kp​w02​]​⎩ ⎨ ⎧​b0​=(Ts2​4​+Ts​4wc​+w02​)b1​=(−Ts2​8​+2w02​)b2​=(Ts2​4​−Ts​4wc​+w02​)​   通常会将差分方程表示成： Y ( z ) X ( z ) = a 0 + a 1 z − 1 + a 2 z − 2 + . . . + a k z − k 1 + b 1 z − 1 + b 2 z − 2 + . . . + b k z − k y [ n ] + b 1 y [ n − 1 ] + b 2 y [ n − 2 ] = a 0 x [ n ] + a 1 x [ n − 1 ] + a 2 x [ n − 2 ] \frac{Y(z)}{X(z)}=\frac{a_0+a_1z^{-1}+a_2z^{-2}+...+a_kz^{-k}}{1+b_1z^{-1}+b_2z^{-2}+...+b_kz^{-k}}\\ y[n]+b_1y[n-1]+b_2y[n-2]=a_0x[n]+a_1x[n-1]+a_2x[n-2]\\ X(z)Y(z)​=1+b1​z−1+b2​z−2+...+bk​z−ka0​+a1​z−1+a2​z−2+...+ak​z−k​y[n]+b1​y[n−1]+b2​y[n−2]=a0​x[n]+a1​x[n−1]+a2​x[n−2]   故而在代码实现中可以看到 a 0 、 a 1 、 a 2 、 b 1 、 b 2 a_0、a_1、a_2、b_1、b_2 a0​、a1​、a2​、b1​、b2​各系数除以 b 0 b_0 b0​。

  比例增益系数 K p K_p Kp​和谐振增益系数 K r K_r Kr​主要影响控制器的增益和相位裕度，截止频率 w c w_c wc​主要影响谐振频率 w o w_o wo​处的带宽，调节 K p K_p Kp​和 K r K_r Kr​可以优化系统的动态性能和稳态性能，调节 w c w_c wc​可以改善系统的抗干扰能力。

PR控制器基于Matlab绘制Bode图 准PR控制器程序