-
摘要: 针对存在临界点的A类被控对象及不存在临界点的B类被控对象, 分别采用其$-180^\circ$和$-120^\circ$相位点的频率和增益提出了PID (Proportional-integral-derivative) 控制器参数的优化整定方法. 基于Tchebyshev多项式和分数阶积分器求取被控对象$-180^\circ$或$-120^\circ$相位点的频率和增益, 建立其积分滞后模型. 采用负载扰动下跟踪误差平方和(Sum of squares of tracking errors, SSE)最小作为优化指标, 使闭环系统具有强的鲁棒性的最大灵敏度和最大补灵敏度为约束方程, 针对两类被控对象, 分别建立了基于$-180^\circ$和$-120^\circ$相位点频率和增益的PID控制器比例、积分与微分三个参数的优化整定规则. 通过与其他常用PID控制方法的仿真与物理对比实验, 表明所提方法的优越性.
-
关键词:
- PID控制 /
- Tchebyshev多项式 /
- 积分滞后模型 /
- 跟踪误差平方和 /
- 优化整定规则
Abstract: In this paper, a tuning method of PID (proportional-integral-derivative) controller parameters is proposed for class A controlled object with critical point and class B controlled object without critical point using the frequency and gain of the point at which its phase is $-180^\circ$ and $-120^\circ$ respectively. The frequency and gain of the controlled object at which its phase is $-180^\circ$ or $-120^\circ$ are obtained based on Tchebyshev polynomial and fractional order integrator, and the gain-integrator-delay model is established for the controlled object. Taking the minimum sum of squares of tracking error (SSE) under load disturbance as the optimization objective, and the maximum sensitivity and maximum complementary sensitivity as the constraint equations that give the closed-loop system strong robustness, a tuning rule of the proportion, integral and differential parameters of PID controller is established for two classes of plants based on the frequency and gain at which the phase is $-180^\circ$ and $-120^\circ$ respectively. Through simulation and physical comparison experiments with other common PID control methods, the superiority of the proposed method is demonstrated. -
比例−积分−微分(Proportional-integral-derivative, PID) 控制技术因其结构简单, 使用方便而广泛应用于实际工业过程, 目前约95%的控制回路采用PID控制技术[1-3]. PID控制器的性能取决于比例、积分和微分三个参数. 文献[4]指出: “··· 没有一种整定方法可以说是通用的或最优的, 迄今为止, 我们仍缺乏足够的知识来确定一个通用且最优的PID控制器参数的整定方法.” 因此PID控制器参数整定方法一直是控制领域研究的重要方向.
目前工业领域广泛采用的PID控制器参数整定技术是Ziegler-Nichols (Z-N)法[5]. Z-N法采用可以反映被控对象动态特性的临界点, 即其Nyquist曲线穿越负实轴的$ -180^\circ $相位点的增益$ K_{180} $和频率$ \omega_{180} $来设计PID控制器参数. Z-N法使闭环系统对阶跃负载扰动的响应以1/4的震荡比衰减, 得到基于临界点增益和频率的简单的PID参数整定规则. Z-N法将PID控制系统的开环传递函数在$ \omega_{180} $点的频率响应调节为$ -0.6-0.28{\rm{j}} $, 在保证控制系统稳定的条件下, 获得好的抗扰性. Z-N法仅适用于时间常数较大的慢动态工业过程[6-7]. 文献[8-9]以给定的PID控制系统开环传递函数在$ \omega_{180} $点的频率响应的幅值和相角为目标, 从而提出了基于临界点增益和频率的PID控制器参数整定规则. 文献[10-11] 采用被控对象在临界点的频率和增益建立积分滞后模型, 文献[12] 采用被控对象的稳态增益与$ -180^\circ $相位点的增益和频率建立一阶惯性滞后或二阶欠阻尼滞后模型, 文献[13]采用被控对象的稳态增益、$ -180^\circ $相位点的增益和频率及该点Nyquist曲线的切角建立二阶滞后模型, 以期望的控制系统传递函数为目标, 提出了基于临界点增益和频率的PID控制器参数整定方法. 为提高PID控制系统的性能和鲁棒稳定性, 文献[14] 针对基于临界点频率和增益的被控对象模型, 以PID控制器的积分增益最大为优化指标, 以敏感度为约束, 提出了基于临界点频率和增益的PID控制器参数整定规则. 但积分增益最大的优化指标只适用于平稳的被控对象, 当负载扰动频繁波动时, 不能减小跟踪误差[9].
上述PID参数整定方法要求被控对象必须存在临界点. 然而低阶或相对阶次较小、延时很小的被控对象不存在临界点. 文献[15-16]采用反映其动态特性的$ -120^\circ $相位点的增益$ K_{120} $和频率$ \omega_{120} $, 以给定的PID控制系统开环传递函数在$ \omega_{120} $点的频率响应的幅值和相角为目标, 提出了PID参数整定规则. 针对离散被控对象, 文献[17-18]提出离散的Z-N法. 文献[19]以给定的PID控制系统开环传递函数在$ \omega_{180} $点的频率响应的幅值与相角为目标, 提出基于临界点的频率和增益的离散PID控制器参数整定方法.
上述PID控制器参数整定方法的$ -180^\circ $或$ -120^\circ $相位点的增益和频率参数采用实验方法近似求取. 当被控对象模型已知时, 文献[18, 20] 提出临界点增益和频率的数学计算方法. 文献[18] 引入比例控制与被控对象组成闭环回路, 通过临界稳定条件下闭环特征方程的根计算使闭环系统临界稳定的比例增益和对应的频率, 得到临界点的增益$ K_{180} $和频率$ \omega_{180} $. 该方法只适合低阶的被控对象. 文献[20]建立被控对象的相位等于$ -180^\circ $的方程, 采用Newton-Raphson算法迭代求解该方程, 计算临界点的增益$ K_{180} $和频率$ \omega_{180} $. 迭代算法时间长且存在计算误差, 相位点增益和频率的求取误差影响PID控制系统的性能.
本文通过引入Tchebyshev多项式和分数阶积分器, 提出被控对象$ -180^\circ $和$ -120^\circ $相位点的频率和增益的求取方法. 为了获得基于相位点频率与增益的PID参数优化整定规则, 本文基于被控对象指定相位点频率和增益的积分滞后模型, 以负载扰动下跟踪误差的平方和(sum of squares of tracking errors, SSE)最小为优化目标, 以最大灵敏度和最大补灵敏度为约束, 提出实现上述优化性能指标的基于被控对象相位点频率和增益的PID控制器参数整定规则. 采用存在临界点和不存在临界点的被控对象的仿真和物理实验, 表明所提方法的优越性.
1. 被控对象模型
$$ A(z^{-1})y(k) = z^{-d}B(z^{-1})u(k) $$ (1) 其中, $ u(k) $为被控对象的输入, $ y(k) $为被控对象的输出, $ d $为系统延时, $ k = t/T_0 = 0,1,2,{\cdots} $为离散的采样时间, $ T_0 $为采样周期, $ A(z^{-1}) $和$ B(z^{-1}) $为关于$ z^{-1} $的多项式, 即为:
$$ \begin{split} A(z^{-1}) =\; &1+\sum_{i = 1}^na_iz^{-i} = 1+a_1z^{-1}+a_2z^{-2} \;+ \\ &{\cdots}+a_nz^{-n} \end{split} $$ (2) $$ \begin{split} B(z^{-1}) = \;&\sum_{i = 1}^nb_iz^{-i} = b_1z^{-1}+b_2z^{-2} \;+ \\ &{\cdots}+b_nz^{-n} \end{split} $$ (3) 其中, $ n $表示系统阶次. 以上被控对象的传递函数为:
$$ G_P (z)= \frac{B(z)}{A(z)} = \frac{b_1z^{n-1}+b_2z^{n-2}+{\cdots}+b_n}{z^d(z^n+a_1z^{n-1}+{\cdots}+a_n)} $$ (4) 本文考虑的被控对象包括两类: 一类是频率响应经过$ -180^\circ $相位线(负实轴)的被控对象, 称为A类被控对象, 这类被控对象存在临界点, 包括大多数存在延时的工业被控对象; 一类是频率响应不经过$ -180^\circ $相位线但经过$ -120^\circ $相位线的被控对象, 称为B类被控对象, 这类被控对象不存在临界点. B类被控对象主要包括稳定的二阶最小相位系统和很多相对阶次小于3的系统[15].
对于A类被控对象, 采用其$ -180^\circ $相位点的数字频率$ \theta_{180} $和增益$ K_{180} $表征被控对象的动态特性, 该点被称为临界点[14]; 对于B类被控对象, 采用其$ -120^\circ $相位点的数字频率$ \theta_{120} $和增益$ K_{120} $描述被控对象的动态特性[15-16]. 为统一描述两类被控对象在指定相位点的频率和增益的求取方法及PID控制器参数的优化, 定义$ \theta_{\varphi} $和$ K_{\varphi} $表示被控对象在指定相位点的频率和增益:
$$ \begin{split} &\theta_{\varphi} = \underset{0< \theta < \pi} {\min \theta }:\angle G_P({\rm{e}}^{{\rm{j}}\theta}) = -\varphi \\& K_{\varphi} = \lvert G_P({\rm{e}}^{{\rm{j}}\theta_{\varphi}}) \rvert \end{split} $$ (5) 若$ G_P(z) $为A类被控对象, $ \theta_{\varphi} $和$ K_{\varphi} $表示其$ -180^\circ $相位点的频率和增益, 即
$$ \varphi = \pi, \theta_{\varphi} = \theta_{180},K_{\varphi} = K_{180} $$ (6) 若$ G_P(z) $为B类被控对象, $ \theta_{\varphi} $和$ K_{\varphi} $表示其$ -120^\circ $相位点的频率和增益, 即
$$ \varphi = \frac{2\pi}{3}, \theta_{\varphi} = \theta_{120},K_{\varphi} = K_{120} $$ (7) 1.1 ${\boldsymbol{-180^\circ }}$和${\boldsymbol{-120^\circ}}$相位点的频率和增益的求取方法
本文引入Tchebyshev多项式[21]和分数阶积分器[22]提出两类被控对象的$ -180^\circ $和$ -120^\circ $相位点的数字频率和增益的求取方法如下. 引入分数阶积分器$ F(z) $, 串联B类被控对象传递函数(4)使其频率响应穿越$ -180^\circ $相位. 对于A类被控对象, $ F(z) = 1 $. 这样可采用统一的方法求取两类对象的频率$ \theta_{\varphi} $和增益$ K_{\varphi} $. 引入Tchebyshev多项式, 将被控对象传递函数关于超越函数$ {\rm{e}}^{{\rm{j}}\theta} $的高次幂的复杂的相位方程转化为关于实变量$ u $的方程, 通过求解该方程即可得到被控对象传递函数的$ -180^\circ $和$ -120^\circ $相位点的数字频率和增益.
定义传递函数:
$$ H(z) = G_P(z)F(z) $$ (8) 若$ G_P(z) $为A类被控对象, 令$ F(z) = 1 $; 若$ G_P(z) $为B类被控对象, 令$ F(z) $为如下的分数阶积分器:
$$ F(z) = \left(\frac{2}{T_0}\frac{z-1}{z+1}\right)^{-\frac{2}{3}} $$ (9) 由文献[15, 23]知, 以上分数阶积分器的相角始终为$ -60^\circ $, 即:
$$ \angle F(z) = -60^\circ $$ (10) 由于分数阶积分器(9)无法直接求解, 采用连分式展开(Continued fraction expansion, CFE) 技术[22]将其近似为整数阶系统计算频率响应:
$$ \begin{split} \widetilde{F}(z)=\;& CFE\{F(z)\} = \\ &\frac{P(z)}{\Pi(z)} = \frac{p_0z^m+p_1z^{m-1}+{\cdots}+p_m}{\alpha_0z^m+\alpha_1z^{m-1}+{\cdots}+\alpha_m} \end{split} $$ (11) 其中, $ CFE\{x\} $表示对$ x $的连续分式展开, $ P(z) $和$ \Pi(z) $为关于$ z $的多项式, $ m $表示近似的阶次, 本文令$ m = 5 $, 由文献[22]知$ P(z) $和$ \Pi(z) $为:
$$ P(z) = p_0z^5+p_1z^{4}+p_2z^3+p_3z^2+p_4z+p_5 $$ (12) 其中,
$$ \begin{split}& p_0 = 25.63\left(\frac{2}{T_0}\right)^{-\frac{2}{3}};\;p_1 = 17.09\left(\frac{2}{T_0}\right)^{-\frac{2}{3}};\\ &p_2 = -23.41\left(\frac{2}{T_0}\right)^{-\frac{2}{3}};\;p_3 = -12.45\left(\frac{2}{T_0}\right)^{-\frac{2}{3}}; \\ &p_4 = 3.83\left(\frac{2}{T_0}\right)^{-\frac{2}{3}};\;p_5 = \left(\frac{2}{T_0}\right)^{-\frac{2}{3}} \\ & \Pi(z) = \alpha_0z^5+\alpha_1z^{4}+\alpha_2z^3+\alpha_3z^2+\alpha_4z+\alpha_5 \end{split} $$ (13) 其中,
$$ \alpha_0 = 25.63;\;\alpha_1 = -17.09;\;\alpha_2 = -23.41; $$ $$ \alpha_3 = 12.45;\;\alpha_4 = 3.83;\;\alpha_5 = -1 $$ 采用近似分数阶积分器(11), $ H(z) $可写为:
$$ \begin{split} H(z)=\;& G_P(z)\tilde{F}(z) = \frac{\beta(z)}{\chi(z)} = \\ &\frac{\beta_1z^{m+n-1}+\beta_2z^{m+n-2}+{\cdots}+\beta_{m+n}}{z^{m+n+d}+\chi_1z^{m+n+d-1}+{\cdots}+\chi_{m+n}z^{d}} \end{split} $$ (14) 通过引入$ F(z) $, 对两类被控对象, 构造的传递函数$ H(z) $必与负实轴相交.
令
$$ u = -{\rm{cos}}\theta $$ (15) 则单位圆上$ z $可以写为[21]:
$$ z = {\rm{e}}^{{\rm{j}}\theta} = {\rm{cos}}\theta+{\rm{j}}{\rm{sin}}\theta = -u+{\rm{j}}\sqrt{1-u^2} $$ (16) 随着$ \theta $从0到$ \pi $, $u $在−1到1之间运动.
由于
$$ z^\kappa = {\rm{e}}^{{\rm{j}}\kappa\theta} = {\rm{cos}}\kappa\theta+{\rm{j}}{\rm{sin}}\kappa\theta $$ (17) 定义:
$$ c_\kappa: = {\rm{cos}}\kappa\theta,s_\kappa: = \frac{{\rm{sin}}\kappa\theta}{{\rm{sin}}\theta} $$ (18) 由文献[21, 24]知, $ c_\kappa $和$ s_\kappa $为关于$ u $的第一类和第二类Tchebyshev多项式. 则$ z^\kappa $等于
$$ z^\kappa = c_\kappa(u)+{\rm{j}}\sqrt{1-u^2} s_\kappa(u) $$ (19) 采用如下的递推公式可求得$ c_\kappa(u) $和$ s_\kappa(u) $:
$$ \begin{split}& s_\kappa(u) = -\frac{1}{\kappa}\frac{\partial c_\kappa(u)}{\partial u},\\ &c_{\kappa+1}(u) = -uc_\kappa(u)-(1-u^2)s_\kappa(u), \kappa = 1,2,\cdots \end{split} $$ (20) 将式(19)代入式(14)可得$ \beta (u) $和$ \chi (u) $:
$$ \begin{align} \beta (u) = R_\beta(u)+{\rm{j}}\sqrt{1-u^2}T_\beta(u) \end{align} $$ (21) $$ \begin{align} \chi (u) = R_\chi(u)+{\rm{j}}\sqrt{1-u^2}T_\chi(u) \end{align} $$ (22) 式中, $ R_\beta(u) $, $ T_\beta(u) $, $ R_\chi(u) $, $ T_\chi(u) $分别为
$$ \begin{split} R_\beta(u) =\; &\beta_1c_{m+n-1}(u)+\cdots+\beta_{m+n-1}c_1(u)\;+\\ &\beta_{m+n} \end{split} $$ (23) $$ T_\beta(u) = \beta_1s_{m+n-1}(u)+\cdots+\beta_{m+n-1}s_1(u) $$ (24) $$ \begin{split} R_\chi(u) =\; &c_{m+n+d}(u)+\cdots+\chi_{m+n-1}c_{d+1}(u)\;+\\ &\chi_{m+n}c_{d}(u) \end{split} $$ (25) $$ \begin{split} T_\chi(u) =\; &s_{m+n+d}(u)+\cdots+ \chi_{m+n-1}s_{d+1}(u)\;+ \\ &\chi_{m+n}s_{d}(u) \end{split} $$ (26) $ H(z) $可以写为:
$$ \begin{split} H(z)|&_{z = -u+{\rm{j}}\sqrt{1-u^2}} = \frac{\beta(z)\overline{\chi}(z)}{\chi(z)\overline{\chi}(z)}\Bigg|_{z = -u+{\rm{j}}\sqrt{1-u^2}}\; = \\ &H(u) = {\rm{Re}}(H(u))+{\rm{j}}{\rm{Im}}(H(u)) \\[-10pt] \end{split} $$ (27) 其中, $ \overline{\chi}(z) $为$ \chi(z) $的复共轭. $ {\rm{Re}}(\cdot) $和$ {\rm{Im}}(\cdot) $分别表示实部和虚部, $ {\rm{Re}}(H(u)) $和$ {\rm{Im}}(H(u)) $为
$$ \begin{align} {\rm{Re}}(H(u))& = \frac{R_\beta(u)R_\chi(u)+(1-u^2)T_\beta(u)T_\chi(u)}{R^2_\chi(u)+(1-u^2)T^2_\chi(u)} \end{align} $$ (28) $$ \begin{align} {\rm{Im}}(H(u))& = \sqrt{1-u^2}\frac{T_\beta(u)R_\chi(u)-R_\beta(u)T_\chi(u)}{{R^2_\chi(u)}+(1-u^2){T^2_\chi(u)}} \end{align} $$ (29) 由于所求$ K_{\varphi} $和$ \theta_{\varphi} $为$ H(z) $的Nyquist曲线与负实轴相交的第一个交点的增益和频率, 定义
$$ u_{\varphi} = -{\rm{cos}}\theta_\varphi $$ (30) 则$ K_{\varphi} $和$u_{\varphi}$满足如下的方程:
$$ \begin{align} {\rm{Re}}(H(u_\varphi))& = -K_{\varphi} \lvert \tilde{F}(u_{\varphi}) \rvert \end{align} $$ (31) $$ \begin{align} {\rm{Im}}(H(u_\varphi))& = 0 \end{align} $$ (32) 由式(29)知
$$ T_\beta(u_\varphi)R_\chi(u_\varphi)-R_\beta(u_\varphi)T_\chi(u_\varphi) = 0 $$ (33) 或者
$$ \sqrt{1-u_\varphi^2} = 0 $$ (34) 由于$ u $的取值范围为 $ [-1,1] $, 且$ u = \pm1 $对应于$ \theta = 0 $和$ \theta = \pi $, 而$ \theta_\varphi\neq0 $且$ \theta_\varphi\neq\pi $, 因此$ u_\varphi $满足
$$ T_\beta(u_\varphi)R_\chi(u_\varphi)-R_\beta(u_\varphi)T_\chi(u_\varphi) = 0 $$ (35) $$ R_\beta(u_\varphi)R_\chi(u_\varphi)+(1-u_\varphi^2)T_\beta(u_\varphi)T_\chi(u_\varphi)<0 $$ (36) $$ -1<u_\varphi<1 $$ (37) 由式(20)和(23) ~ (26)知, 由于传递函数$ G_P(z) $和$ \tilde{F}(z) $系数可知, 式(35)和(36) 所有系数已知. 在$ H(z) $的Nyquist曲线上, 可能存在多个虚部为0的点, 而$ \theta_\varphi $为其中使$ {\rm{Re}}(H)<0 $ 的最小频率. 由式(15)知, $ u $为关于$ \theta $的单调递增函数, 因此$ u_\varphi $为$ (-1,1) $区间内使式(36)小于0的式(35)的最小解.
由式(30)可得
$$ \theta_\varphi = {\rm{arccos}}(-u_\varphi) $$ (38) 将$ u_\varphi $代入式(28)可得$ {\rm{Re}}(H(u_\varphi)) $, 当被控对象为A类, 由式(31)得
$$ K_{180} = -{\rm{Re}}(H(u_\varphi)) $$ (39) 当被控对象为B类, 将$ u_\varphi $代入式(19), 并由式(11) ~ (13)求得$ \tilde{F}(u_\varphi) $为:
$$ \tilde{F}(u_\varphi) = {\rm{Re}}(\tilde{F}(u_\varphi))+{\rm{j}}{\rm{Im}}(\tilde{F}(u_\varphi)) $$ (40) 其中, $ {\rm{Re}}(\tilde{F}(u_\varphi)) $和$ {\rm{Im}}(\tilde{F}(u_\varphi)) $为
$$ \begin{aligned} {\rm{Re}}(\tilde{F}&(u_\varphi)) = \\ &\frac{R_P(u_\varphi)R_\Pi(u_\varphi)+(1-u_\varphi^2)T_P(u_\varphi)T_\Pi(u_\varphi)}{{R^2_\Pi(u_\varphi)}+(1-u_\varphi^2){T^2_\Pi(u_\varphi)}} \nonumber \end{aligned} $$ $$ \begin{aligned}{\rm{Im}}(\tilde{F}(u_\varphi))& = \\ \sqrt{1-u_\varphi^2} & \frac{T_P(u_\varphi)R_\Pi(u_\varphi)-R_P(u_\varphi)T_\Pi(u_\varphi)}{{R^2_\Pi(u_\varphi)}+(1-u_\varphi^2){T^2_\Pi(u_\varphi)}} \nonumber \end{aligned} $$ 其中,
$$ \begin{aligned} R_P(u_\varphi) =\; &p_0c_5(u_\varphi)+ p_1c_4(u_\varphi)+p_2c_3(u_\varphi)\;+ \\ &p_3c_2(u_\varphi)+p_4c_1(u_\varphi)+p_5 \nonumber \end{aligned} $$ $$ \begin{aligned} T_P(u_\varphi) = \;&p_0s_5(u_\varphi)+ p_1s_4(u_\varphi)+p_2s_3(u_\varphi)\;+ \\ &p_3s_2(u_\varphi)+p_4s_1(u_\varphi) \nonumber \end{aligned} $$ $$ \begin{aligned} R_\Pi(u_\varphi) =\; &\alpha_0c_5(u_\varphi)+ \alpha_1c_4(u_\varphi)+\alpha_2c_3(u_\varphi)\;+ \\ &\alpha_3c_2(u_\varphi)+\alpha_4c_1(u_\varphi)+\alpha_5 \nonumber \end{aligned} $$ $$ \begin{aligned} T_\Pi(u_\varphi) =\; &\alpha_0s_5(u_\varphi)+ \alpha_1s_4(u_\varphi)+\alpha_2s_3(u_\varphi)\;+ \\ &\alpha_3s_2(u_\varphi)+\alpha_4s_1(u_\varphi) \nonumber \end{aligned} $$ 由式(31)得
$$ K_{120} = -\frac{{\rm{Re}}(H(u_\varphi))}{\lvert \tilde{F}(u_{\varphi}) \rvert} $$ (41) 其中,
$$ \lvert \tilde{F}(u_{\varphi}) \rvert = \sqrt{{\rm{Re}}^2(\tilde{F}(u_{\varphi}))+{\rm{Im}}^2(\tilde{F}(u_{\varphi}))} $$ (42) 1.2 基于指定相位点频率和增益的被控对象模型
由1.1节得到的两类被控对象在$ -180^\circ $相位点或$ -120^\circ $相位点的频率和增益建立积分滞后模型. 为方便PID控制器参数优化指标的设计, 将两类被控对象的离散积分滞后模型统一描述为[11, 25-26]:
$$ \tilde{G}_P(z) = \frac{K_{\varphi}\lvert {\rm{e}}^{{\rm{j}}\theta_\varphi}-1 \rvert}{z^{\frac{\pi}{\eta\theta_{{\varphi}}}}(z-1)} $$ (43) 其中, $ K_{\varphi} $和$ \theta_\varphi $由式(5) ~ (7)给出, 分别为$ -180^\circ $或$ -120^\circ $相位点的增益和频率. $ \eta $决定了两类被控对象的积分滞后模型的延时项的相位滞后, 若$ G_P $为A类被控对象, $ \eta = 2 $; 若$ G_P $为B类被控对象, $ \eta = 6 $.
对于积分滞后模型(43), 若$ G_P(z) $为A类被控对象, 在$ \theta_{180} $点, 有
$$ \begin{split} &\angle\tilde{G}_P({\rm{e}}^{{\rm{j}}\theta_{180}}) = \angle {\rm{e}}^{-{\rm{j}}\frac{\pi}{2}}-\frac{\pi}{2} = -\pi \\& \lvert \tilde{G}_P({\rm{e}}^{{\rm{j}}\theta_{180}}) \rvert = K_{180} \end{split} $$ (44) 即积分滞后模型(43)与被控对象的传递函数(4) 在$ \theta_{180} $点的频率响应一致. 同理, 若$ G_P(z) $为B类被控对象, 在$ \theta_{120} $点, 有
$$ \begin{split} & \angle\tilde{G}_P({\rm{e}}^{{\rm{j}}\theta_{120}})= \angle {\rm{e}}^{-{\rm{j}}\frac{\pi}{6}}-\frac{\pi}{2} = -\frac{2\pi}{3} \\ &\lvert \tilde{G}_P({\rm{e}}^{{\rm{j}}\theta_{120}}) \rvert = K_{120} \end{split} $$ (45) 此时积分滞后模型(43)与被控对象的传递函数(4)在$ \theta_{120} $点的频率响应一致.
采用基于指定相位点频率和增益的积分滞后模型(43)的目的是设计优化的PID控制器参数, 式(43)不需要物理实现, 因此该模型不需要必须为有理真分式.
2. PID控制器参数的优化整定方法
2.1 PID控制器参数的优化整定策略
采用被控对象的积分滞后模型(43), 以负载扰动下跟踪误差的平方和最小作为优化指标, 闭环系统的最大灵敏度在1.2到2之间[14, 27], 最大补灵敏度在1到1.5之间[14, 27], 及PID控制器的微分时间常数$ T_D $等于积分时间常数$ T_I $的四分之一[28]作为约束方程, 求取基于被控对象指定相位点频率和增益的PID控制器参数整定规则.
首先利用负载扰动的低频特性, 跟踪误差对负载扰动的响应主要取决于PID控制器[27], 得到负载扰动下跟踪误差平方和的近似表达式.
为了使PID控制器参数的选择与被控对象指定相位点的频率$ \theta_{\varphi} $和增益$ K_{\varphi} $相联系, 引入参数$ \rho_K $和$ \rho_T $使$ K_P = \rho_K/K_{\varphi} $, $ T_I = 2\pi\rho_T T_0/\theta_{\varphi} $. 采用序列二次规划算法求解实现上述优化指标的基于$ \theta_{\varphi} $的$ \rho_K $和$ \rho_T $, 从而建立基于被控对象指定相位点频率和增益的PID控制器参数优化的整定规则.
2.2 PID参数整定的优化性能指标
PID控制器为:
$$ \begin{split} u(k) =\; &u(k-1)+K_P[e(k)-e(k-1)]+\frac{K_PT_0}{T_I}e(k)\;+\\ &\frac{K_PT_D}{T_0}[e(k)-2e(k-1)+e(k-2)]\\[-15pt]\end{split} $$ (46) 其中, $ K_P $为PID控制器的比例增益, $ T_I $和$ T_D $为积分时间常数和微分时间常数, $ T_0 $为采样周期. 上述PID控制器的传递函数为:
$$ \begin{aligned} G_c(z) = K_P \frac{(1+\frac{T_0}{T_I}+\frac{T_D}{T_0})z^2+(-1-2\frac{T_D}{T_0})z + \frac{T_D}{T_0}}{z(z-1)} \end{aligned} $$ (47) 在负载扰动$ \upsilon(k) $下, 被控对象模型(43)在PID控制下的闭环控制系统如图1所示. 其中$ y_{sp}(k) $为设定值, $ u(k) $为被控对象输入, $ y(k) $为被控过程输出, $ e(k) $为跟踪误差, $ e(k) = y_{sp}(k)-y(k) $.
定义$ e_{\upsilon}(k) $为扰动信号$ \upsilon(k) $下的跟踪误差, 即$ y_{sp}(k) = 0 $时, 跟踪误差$ e(k) $对扰动信号$ \upsilon(k) $的响应. 扰动$ \upsilon(k) $到跟踪误差$ e(k) $的传递函数为
$$ G_\upsilon(z) = -\frac{1}{G_c(z)}\frac{G_c(z)\tilde{G}_P(z)}{1+G_c(z)\tilde{G}_P(z)} $$ (48) 由于PID控制器包含积分, 因此上式中,
$$ \underset {z \rightarrow 1}{{\rm{lim}}} \frac{G_c(z)\tilde{G}_P(z)}{1+G_c(z)\tilde{G}_P(z)} = 1 $$ (49) 由于负载扰动$ \upsilon(k) $通常是低频的[27], 由式(48)和(49)知, 跟踪误差$ e(k) $对负载扰动$ \upsilon(k) $的响应主要取决PID控制器. 因此扰动$ \upsilon(k) $到跟踪误差$ e(k) $的传递函数可近似表示为[27]:
$$ G_\upsilon(z)\approx-\frac{1}{G_c(z)} $$ (50) 在采用指定相位点频率和增益的PID控制器参数整定方法[5, 8-12, 14-18]中, PID控制器参数$ K_P $、$ T_I $和$ T_D $可写为关于指定相位点增益$ K_{\varphi} $的倒数和周期$ T_{\varphi}\;(T_{\varphi} = 2\pi T_0/\theta_\varphi) $的线性表达式. 因此本文为了使PID控制器参数的选择与被控对象在指定相位点的频率$ \theta_\varphi $和增益$ K_{\varphi} $相联系, 引入无量纲参数$ \rho_K $和$ \rho_T $, 并利用约束方程$ T_D = T_I/4 $, 将PID控制器参数表示为:
$$ K_P = \frac{\rho_K}{K_{\varphi}} $$ (51) $$ T_I = \rho_T T_{\varphi} = \frac{2\pi\rho_T}{\theta_\varphi}T_0 $$ (52) $$ T_D = \frac{1}{4}\rho_T T_{\varphi} = \frac{ \pi \rho_T}{2\theta_\varphi}T_0 $$ (53) 令$h = \theta_\varphi/(\rho_T\pi)$, 由式(47)知, $ G_c(z) $为:
$$ G_c(z) = \frac{\rho_K[(h^2+2h+1)z^2+(-2h-2)z+1]}{2hK_\varphi z(z-1)} $$ (54) 将其代入式(50), 得:
$$ \begin{split} G_\upsilon(z)\approx \;& -\frac{2K_\varphi hz(z-1)}{\rho_K[(z-1)^2+2hz(z-1)+h^2z^2]} = \\ &-\frac{2K_\varphi hz(z-1)}{\rho_K(1+h)^2(z-\frac{1}{1+h})^2}\\[-15pt] \end{split} $$ (55) 由文献[29]知当扰动$ \upsilon(k) $为单位阶跃信号, 则$ e_\upsilon(k) $为:
$$ \begin{split} e_\upsilon(k)=\;& Z^{-1}\left\{\frac{z}{z-1}G_\upsilon(z)\right\} \approx \\ &Z^{-1}\left\{-\frac{2K_\varphi hz^2}{\rho_K(1+h)^2(z-\frac{1}{1+h})^2}\right\} = \\ &-\frac{2K_\varphi h(k+1)}{\rho_K(1+h)^2(1+h)^k} \\[-15pt]\end{split} $$ (56) 其中$ Z^{-1} $为$z $逆变换. 则闭环系统在单位阶跃扰动下的跟踪误差的平方和为:
$$ \begin{split} SSE_\upsilon = \;&\sum_{k = 0}^\infty [e_\upsilon(k)]^2\approx \\ &\frac{4K_{\varphi}^2}{\rho_K^2}\sum_{k = 0}^\infty\left[\frac{h(k+1)}{(1+h)^2(1+h)^k}\right]^2 = \\ &\frac{4K_{\varphi}^2 (h^2+2h+2)}{\rho_K^2h(h+2)^3} = \\ &\frac{4K_{\varphi}^2 \left[(\frac{\theta_\varphi}{\pi\rho_T})^2+2(\frac{\theta_\varphi}{\pi\rho_T})+2\right]}{\rho_K^2(\frac{\theta_\varphi}{\pi\rho_T})(\frac{\theta_\varphi}{\pi\rho_T}+2)^3} \end{split} $$ (57) 上式中, $ K_{\varphi} $只决定$ SSE_\upsilon $指标的值, 不影响$ \rho_K $和$ \rho_T $的优化解, 因此选择优化指标为:
$$ J_1= \frac{SSE_\upsilon}{K_{\varphi}^2} = \frac{4\left[(\frac{\theta_\varphi}{\pi\rho_T})^2+2(\frac{\theta_\varphi}{\pi\rho_T})+2\right]}{\rho_K^2(\frac{\theta_\varphi}{\pi\rho_T})(\frac{\theta_\varphi}{\pi\rho_T}+2)^3} $$ (58) 以$ \rho_K $和$ \rho_T $为决策变量, 式(58)为优化指标的PID参数整定的优化模型为:
$$ \underset {\rho_K,\rho_T\in {\bf{R}} ^+}{{\rm{min}}}{J_1 = \frac{4\left[(\frac{\theta_\varphi}{\pi\rho_T})^2+2(\frac{\theta_\varphi}{\pi\rho_T})+2\right]}{\rho_K^2(\frac{\theta_\varphi}{\pi\rho_T})(\frac{\theta_\varphi}{\pi\rho{}_T}+2)^3}} $$ (59) $$ {\rm{s.t}}.\quad M_s\leq M_s^{ub} \qquad \qquad \qquad \quad\;\;\;\quad$$ (60) $$\;\;\;\;M_t \leq M_t^{ub} \quad \;\; \qquad \qquad\qquad $$ (61) $$ \begin{align} &T_D = \frac{T_I}{4} \end{align} $$ (62) 式中$ M_s $和$ M_t $分别为闭环系统的最大灵敏度和最大补灵敏度, 由文献[14]知$ M_s $和$ M_t $分别为:
$$ \begin{align} M_s = \left\|\frac{1}{1+\tilde{G}_{ol}({\rm{e}}^{{\rm{j}}\theta})}\right\| _\infty \end{align} $$ (63) $$ \begin{align} M_t = \left\|\frac{\tilde{G}_{ol}({\rm{e}}^{{\rm{j}}\theta})}{1+\tilde{G}_{ol}({\rm{e}}^{{\rm{j}}\theta})}\right\| _\infty \end{align} $$ (64) 其中$ \parallel \cdot \parallel_\infty $表示 “·” 的$ H_\infty $范数, $ \tilde{G}_{ol} $为图1所示的PID闭环控制系统的开环传递函数, 由式(43)和(54)知$ \tilde{G}_{ol}(z) $为
$$ \begin{split} \tilde{G}_{ol}(z)=\;& G_c(z) \tilde{G}_{P}(z) = \\ &\rho_K\lvert {\rm{e}}^{{\rm{j}}\theta_\varphi}-1 \rvert\;\times \\ &\frac{[(h^2+2h+1)z^2+(-2h-2)z+1]}{2hz(z-1)^2z^{\frac{\pi}{\eta\theta_{{\varphi}}}}} = \\ &\rho_K\lvert {\rm{e}}^{{\rm{j}}\theta_\varphi}-1 \rvert \;\times\\ &\frac{[(\frac{\theta_\varphi^2}{\pi^2\rho_T^2}+\frac{2\theta_\varphi}{\pi\rho_T}+1)z^2+(-\frac{2\theta_\varphi}{\pi\rho_T}-2)z+1]}{\frac{2\theta_\varphi}{\pi\rho_T}z(z-1)^2z^{\frac{\pi}{\eta\theta_{{\varphi}}}}} \end{split} $$ (65) $ M_s^{ub} $ 和 $ M_t^{ub} $ 为 $ M_s $ 和 $ M_t $ 的上界, 设 $ M_s^{ub} = 1.7 $, $ M_t^{ub} = 1.5 $.
2.3 优化PID控制器参数的整定规则
由于式(58)和(65)中所采用的被控对象的数字频率$ \theta_\varphi $的取值范围为$ (0,\pi) $, 为了求出使$ J_1 $最小的$ \rho_K $和$ \rho_T $, 将$ \theta_\varphi $按间隔0.02网格化, 对每个网格点, 使用序列二次规划(Sequential quadratic programming, SQP) 算法[30] 求性能指标(59)最小的最优解, 图2和图3分别表示A类被控对象与B类被控对象$ \rho_K $和$ \rho_T $的最优解曲线. 采用$ \theta_\varphi $和相应的$ \rho_K $与$ \rho_T $的最优解数据建立以$ \theta_\varphi $为输入, 相应的$ \rho_K $和$ \rho_T $的最优解为输出的多项式. 采用基于AIC (Akaike＇s information criterion)准则的多项式阶次确定方法[31]和最小二乘参数设计方法[32], 得A类被控对象和B类被控对象使性能指标(59)最小的$ \rho_K $和$ \rho_T $如式(66) ~ (69)所示.
图 2 A类被控对象$ \rho_K $和$ \rho_T $的最优解曲线Fig. 2 Optimal solution curve of $ \rho_K $ and $ \rho_T $ of class A controlled object
图 3 B类被控对象$ \rho_K $和$ \rho_T $的最优解曲线Fig. 3 Optimal solution curve of $ \rho_K $ and $ \rho_T $ of class B controlled objectA类被控对象:
$$ \rho_K = -0.02\theta_{180}^3+0.15\theta_{180}^2-0.34\theta_{180}+0.39 $$ (66) $$ \rho_T = 0.45\theta_{180}+0.65 $$ (67) B类被控对象:
$$ \rho_K = -0.04\theta_{120}^3+0.28\theta_{120}^2-0.65\theta_{120}+0.67 $$ (68) $$ \rho_T = 0.39\theta_{120}+0.25 $$ (69) 由式(51) ~ (53)知, PID控制器参数的整定规则为:
$$ K_P = \frac{\rho_K}{K_{\varphi}} $$ (70) $$ T_I = \frac{2\pi\rho_T}{\theta_\varphi}T_0 $$ (71) $$ T_D = \frac{T_I}{4} = \frac{\pi\rho_T}{2\theta_\varphi}T_0 $$ (72) A类被控对象$ K_{\varphi} = K_{180} $, $ \theta_\varphi = \theta_{180} $; B类被控对象$ K_{\varphi} = K_{120}, \theta_\varphi = \theta_{120} $.
PID控制器参数整定算法如下:
1) 采用第1.1节的相位点频率$ \theta_\varphi $和增益$ K_\varphi $的求取方法, 针对A类对象求得$ \theta_{180} $和$ K_{180} $, 针对B类对象求得$ \theta_{120} $和$ K_{120} $.
2) 针对A类对象由式(66)和(67)求$ \rho_K $和$ \rho_T $, 针对B类对象由式(68)和(69)求$ \rho_K $和$ \rho_T $, 由式(70) ~ (72)求取PID控制器参数$ K_P $、$ T_I $和$ T_D $.
注 1. 由文献[33]的定理3.5可知, 若实际的被控对象传递函数(4)在单位圆外没有极点, 且PID控制器设计模型(43)和实际的被控对象传递函数(4) 的频率响应满足
$$ \left| \frac{G_P({\rm{e}}^{{\rm{j}}\theta})-\tilde{G}_P({\rm{e}}^{{\rm{j}}\theta})}{\tilde{G}_P({\rm{e}}^{{\rm{j}}\theta})} \right|<\left| \frac{1}{T} \right| = \left| \frac{1+\tilde{G}_{ol}({\rm{e}}^{{\rm{j}}\theta})}{\tilde{G}_{ol}({\rm{e}}^{{\rm{j}}\theta})} \right| $$ (73) 则采用优化指标(59) ~ (62)所提出的PID控制器参数整定方法(66) ~ (72)所设计的PID控制器可以使实际的被控对象闭环稳定.
注 2. 上述PID参数优化整定规则(式(66) ~ (72))的关键是求取被控对象的相位点参数$ K_\varphi $和$ \theta_\varphi $. 本文给出了被控对象模型已知时的计算方法, 对于难以建立数学模型的被控对象, 可采用阶跃实验法获得[5].
3. 仿真实验
采用PID控制器参数的离散Z-N法[18]、不存在临界点的被控对象的PID整定方法[15]、Astrom的PID参数优化方法[14]及抗干扰PID整定方法[7]与本文的PID优化整定方法进行仿真对比实验.
例 1. 采用文献[18]的仿真模型:
$$ G_{P1}(z) = \frac{0.0329z^{-1}+0.0269z^{-2}}{1-1.4891z^{-1}+0.5488z^{-2}} $$ (74) 系统采样周期为$ T_0 = 2\;{\rm{s}} $. 设定值$ y_{sp}(t) $为:
$$ y_{sp}(t) = \begin{cases} 1,\quad\;\;\;\delta T<t\leq\delta T + \dfrac{T}{2} \\ -1,\quad \delta T + \dfrac{T}{2}<t\leq(\delta+1)T \end{cases} $$ (75) 其中, $ \delta = 0,1 $, $ T $为方波信号的周期, 设$T = 800\;{\rm{s}}$. 负载扰动$ \upsilon(t) $为:
$$ \upsilon(t) = \begin{cases} 0.2,&200\;{\rm{s}}<t\leq400\;{\rm{s}} \\ -0.5,&400\;{\rm{s}}<t\leq800\;{\rm{s}} \\ 0.2{\rm{cos}}(t),&800\;{\rm{s}}<t \end{cases} $$ (76) 采用文献[18]的方法求得被控对象(74)的$ -180^\circ $相位点的增益$ K_{180} $和周期$ T_{180} $($ T_{180} = 2\pi T_0/\theta_{180} $) 为0.0596和11.6027 (对应$ \theta_{180} = 1.0827 $), PID控制器参数$ K_P $、$ T_I $和$ T_D $为:
$$ \begin{split}& K_P = \frac{0.6}{K_{180}} = 10.0671 \\ &T_I = 0.5T_{180} = 5.8014 \\& T_D = 0.125T_{180} = 1.4503 \end{split} $$ (77) 采用本文方法求得被控对象(74)的$ -180^\circ $相位点的增益$ K_{180} $和数字频率$ \theta_{180} $为0.0596和1.0827, $ \rho_K $和$ \rho_T $为0.1698和1.1320, PID控制器参数$ K_P $、$ T_I $和$ T_D $为:
$$ \begin{split}& K_P = \frac{\rho_K}{{K}_{180}} = 2.8490 \\& T_I = \frac{2\pi\rho_T}{\theta_{180}}\; T_0= 13.1319 \\& T_D =\frac{ T_I}{4} = 3.2830 \end{split} $$ (78) 采用文献[18]方法与本文方法的设定值$ y_{sp}(t) $、被控对象输出$ y(t) $及输入$ u(t) $的曲线分别如图4和图5所示. 采用下列误差绝对值累积和(Sum of absolute error, SAE)[34]和均方误差(Mean-square error, MSE)[35]对两种方法的评价结果如表1所示.
图 4 Z-N法的设定值$y_{sp}(t)$与输出$y(t)$及输入$u(t)$的曲线Fig. 4 Curves of set value $y_{sp}(t)$, output $y(t)$ and input $u(t)$ with Z-N method
图 5 本文方法的设定值$y_{sp}(t)$与输出$y(t)$及输入$u(t)$的曲线Fig. 5 Curves of set value $y_{sp}(t)$, output $y(t)$ and input $u(t)$ with the proposed method表 1 本文方法与Z-N法下误差$ e(t)$的性能评价Table 1 Performance evaluation of error $ e(t)$ with the proposed method and Z-N method整定方法 性能指标 SAE MSE 本文 3835.3081 0.0274 Z-N法 8004.8341 0.0479 $$ \begin{align} SAE = &\sum\limits_{k = 1}^N \lvert y_{sp}(k)-y(k) \rvert \end{align} $$ (79) $$ \begin{align} MSE = &\frac{1}{N}\sum\limits_{k = 1}^N [ y_{sp}(k)-y(k)]^2 \end{align} $$ (80) 其中$ N = 800 $.
由图4、5和表1可以看出, 虽然采用本文方法求取的模型(74)的$ -180^\circ $相位点增益和周期与文献[18]方法相同, 但控制效果明显优于Z-N法, SAE和MSE分别降低了52.09%和42.80%. 这是因为Z-N法下闭环系统具有过大的最大灵敏度$ M_s $和最大补灵敏度$ M_t $, 分别为4.81和4.36, 远超出$ M_s $和$ M_t $的理想范围[14, 27]. 本文方法的$ M_s $和$ M_t $为1.42和1.00, 因此在设定值和扰动频繁变化下本文方法使闭环系统保持好的动态性能.
例 2. 采用文献[15]的B类被控对象的仿真模型:
$$ {G_{P2}(s) = \frac{4.7\times10^5s^2+1.4\times10^7s+4.67\times10^{10}}{s^4+178s^3+1.09\times10^6s^2+4.21\times10^7s+4.67\times10^{10}} } $$ (81) 图6为以上被控对象的Nyquist曲线, 可以看到, 其Nyquist曲线不穿越负实轴. 设定值$ y_{sp}(t) $为:
图 6 被控对象(81)的Nyquist曲线$(K_{120} = 2.26,$ $\omega_{120} = 224.25)$Fig. 6 Nyquist curve of controlled plant (81) $(K_{120} = 2.26,$ $\omega_{120} = 224.25)$$$ y_{sp}(t) = \begin{cases} 1,\quad\;\;\;\delta T<t\leq\delta T +\dfrac{T}{2} \\ -1,\quad \delta T + \dfrac{T}{2}<t\leq(\delta+1)T \end{cases} $$ (82) 设$ \delta = 0,1 $, $T = 10\;{\rm{s}}$. 负载扰动$ \upsilon(t) $为:
$$ \upsilon(t) = \begin{cases} -0.2,&2\;{\rm{s}}<t\leq5\;{\rm{s}} \\ -0.5,&5\;{\rm{s}}<t\leq8\;{\rm{s}} \\ 0.2{\rm{cos}}(5t),&8\;{\rm{s}}<t\\ \end{cases} $$ (83) 采用文献[15]的方法求得被控对象(81)的$ -120^\circ $相位点的增益$ K_{120} $和频率$ \omega_{120} $为2.3581和232, PID控制器参数$ K_P $、$ T_I $和$ T_D $为:
$$ \begin{split} &K_P = \frac{{\rm{cos}}^2(10^\circ)}{K_{120}} = 0.41 \\ &T_I =\frac{ 1}{\omega_{120}{\rm{tan}}(10^\circ) }= 0.024 \\ &T_D =\frac{ {\rm{tan}}(10^\circ)}{\omega_{120} }= 7.6\times10^{-4} \end{split} $$ (84) 设置采样周期$T_0=0.001\;{\rm{s}}$, 将被控对象(81)离散化. 采用本文方法求得其$ -120^\circ $相位点的增益$ K_{120} $和频率$ \omega_{120} $为2.2856和223.9, 数字频率$ \theta_{120} $$( \theta_{120} = \omega_{120}T_0) $ 为0.2239, $ \rho_K $和$ \rho_T $为0.5425和0.3335, PID控制器参数$ K_P $、$ T_I $和$ T_D $为:
$$ \begin{split}& K_P = \frac{\rho_K}{K_{120}} = 0.2374 \\ &T_I = \frac{2\pi\rho_TT_0}{\theta_{120}} = 0.0094 \\ &T_D = \frac{T_I}{4 }= 0.0024 \end{split} $$ (85) 采用文献[15]方法与本文方法的设定值$ y_{sp}(t) $、被控对象输出$ y(t) $及输入$ u(t) $的曲线分别如图7和图8所示. 采用式(79)的SAE和式(80)的MSE对两种方法的评价结果如表2所示, 其中$N = \;20\,000$.
图 8 本文方法的设定值$y_{sp}(t)$与输出$y(t)$及输入$u(t)$的曲线Fig. 8 Curves of set value $y_{sp}(t)$, output $y(t)$ and input $u(t)$ with the proposed method由图7、8和表2可以看出, 本文方法的控制效果优于文献[15]的方法, SAE和MSE分别降低了35.98%和40.79%. 这是因为文献[15]方法$ -120^\circ $相位点的增益$ K_{120} $和频率$ \omega_{120} $是采用继电反馈实验近似求取, 与本文方法相比, 存在较大的误差. 由图6知被控对象(81)的$ -120^\circ $相位点的增益$ K_{120} $和频率$ \omega_{120} $为$ K_{120} = 2.26 $, $ \omega_{120} = 224.25 $, 文献[15]方法的误差为4.16%和3.34%, 本文方法的误差为1.12% 和0.16%. 此外文献[15]基于给定的PID控制系统开环传递函数的相角裕度的整定方法无法有效保证闭环系统的抗扰性[36]. 本文通过使负载扰动下跟踪误差平方和最小, 使得闭环系统具有强的抗扰能力.
例 3. 采用文献[14]的仿真模型:
$$ \begin{aligned} G_{P3} (s)= \frac{1}{(s+1)(0.1s+1)(0.01s+1)(0.001s+1)} \end{aligned} $$ (86) 设定值$ y_{sp}(t) $为
$$ y_{sp}(t) = \begin{cases} 1,&\delta T<t\leq\delta T + \dfrac{T}{2} \\ -1,& \delta T +\dfrac{T}{2}<t\leq(\delta+1)T \end{cases} $$ (87) 设$ \delta = 0,1 $, $T = 10\;{\rm{s}}$. 负载扰动$ \upsilon(t) $为:
$$ \upsilon(t) = \begin{cases} -0.5,&2\;{\rm{s}}<t\leq5\;{\rm{s}} \\ {\rm{cos}}(20t),&5\;{\rm{s}}<t\leq10\;{\rm{s}} \\ 2{\rm{sin}}(t){\rm{cos}}(20t),&10\;{\rm{s}}<t \end{cases} $$ (88) 采用文献[14]的方法得到被控对象(86)的稳态增益$ K_0 $为1, $ -180^\circ $相位点的增益$ K_{180} $和周期$ T_{180} $为0.0091和0.199, PID控制器参数$ K_P $、$ T_I $和$ T_D $为:
$$ \begin{split} &K_P =\frac{ 0.3-0.1\lambda^4}{K_{180} }= 32.8670 \\ &T_I = \frac{0.6T_{180}}{1+2\lambda} = 0.1173 \\& T_D = \frac{0.15(1-\lambda)T_{180}}{1-0.95\lambda} = 0.0298 \end{split} $$ (89) 其中$ \lambda = K_{180}/K_0 $.
设置采样周期$T_0=0.01\;{\rm{s}}$, 将被控对象(86)离散化. 采用本文方法求得其$ -180^\circ $相位点的增益$ K_{180} $和周期$ T_{180} $为0.0135和0.2003, 数字频率$ \theta_{180} $ $ (\theta_{180}= 2\pi T_0/T_{180}) $ 为0.3136, $ \rho_K $和$ \rho_T $为0.2989和0.7861, PID控制器参数$ K_P $、$ T_I $和$ T_D $为
$$ \begin{split}& K_P = \frac{\rho_K}{K_{180}} = 22.1407 \\ &T_I = \frac{2\pi\rho_TT_0}{\theta_{180} }= 0.1574 \\ &T_D = \frac{T_I}{4} = 0.0394 \end{split} $$ (90) 采用文献[14]方法与本文方法的设定值$ y_{sp}(t) $、被控对象输出$ y(t) $及输入$ u(t) $的曲线分别如图9和图10所示. 采用式(79)的SAE和式(80)的MSE对两种方法的评价结果如表3所示, 其中$ N = 2\,000 $.
图 10 本文方法的设定值$y_{sp}(t)$与输出$y(t)$及输入$u(t)$的曲线Fig. 10 Curves of set value $y_{sp}(t)$, output $y(t)$ and input $u(t)$ with the proposed method由图9、10和表3可以看出, 本文方法的控制效果优于文献[14]的方法, SAE和MSE分别降低了19.21%和11.75%. 这是因为Astrom的PID参数优化方法使负载扰动下跟踪误差的积分最小的优化指标只适用于平稳的被控过程, 当负载扰动频繁变化, 不能降低跟踪误差[9, 27]. 本文的跟踪误差平方和最小的优化指标可以在负载扰动频繁变化时有效降低跟踪误差, 使闭环系统具有好的跟踪性能.
例 4. 采用文献[7]的仿真模型:
$$ \begin{aligned} G_{P4} (s)= \frac{-1.4s+1}{(s+1)^3} \end{aligned} $$ (91) 设定值$ y_{sp}(t) $为
$$ y_{sp}(t) = \begin{cases} 2,&\delta T<t\leq\delta T + \dfrac{T}{2} \\ -2,& \delta T + \dfrac{T}{2}<t\leq(\delta+1)T \end{cases} $$ (92) 设$ \delta = 0,1 $, $T = 800\;{\rm{s}}$. 负载扰动$ \upsilon(t) $为:
$$ \upsilon(t) = \begin{cases} 0.2,&100\;{\rm{s}}<t\leq500\;{\rm{s}} \\ 0.2{\rm{cos}}(0.5t),&500\;{\rm{s}}<t\leq800\;{\rm{s}} \\ 0.3{\rm{sin}}(0.1t){\rm{cos}}(0.5t),&800\;{\rm{s}}<t \end{cases} $$ (93) 采用文献[7]的方法求得$ a_0/b_0 = 1 $, 扰动观测带宽$ \omega_{0} $为0.356, 控制带宽$ \omega_c $为0.58, 调节参数$ \alpha $为0.5, 得到PID控制器参数$ K_P $、$ T_I $和$ T_D $为:
$$ \begin{split}& K_P =\frac{ a_0\omega_{0}(1+\alpha)}{b_0 \omega_c}= 0.92 \\& T_I =\frac{ \alpha+1}{\omega_c }= 2.59 \\ &T_D = \frac{\alpha}{(1+\alpha)\omega_c} = 0.57 \end{split} $$ (94) 设置采样周期$T_0 = 0.1\;{\rm{s}}$, 将被控对象(91)离散化. 采用本文方法求得其$ -180^\circ $相位点的增益$ K_{180} $和数字频率$ \theta_{180} $为0.6608和0.0899, $ \rho_K $和$ \rho_T $为0.3622和0.6854, PID控制器参数$ K_P $、$ T_I $和$ T_D $为:
$$ \begin{split} &K_P = \frac{\rho_K}{K_{180}} = 0.5481 \\ &T_I =\frac{2\pi\rho_TT_0}{\theta_{180}} = 4.7879 \\ &T_D = \frac{T_I}{4} = 1.1970 \end{split} $$ (95) 采用文献[7]方法与本文方法的设定值$ y_{sp}(t) $、被控对象输出$ y(t) $及输入$ u(t) $的曲线分别如图11和图12所示. 采用式(79)的SAE和式(80)的MSE对两种方法的评价结果如表4所示, 其中$ N = 16\,000 $.
图 12 本文方法的设定值$y_{sp}(t)$与输出$y(t)$及输入$u(t)$的曲线Fig. 12 Curves of set value $y_{sp}(t)$, output $y(t)$ and input $u(t)$ with the proposed method由图11、12和表4可以看出, 本文方法的控制效果优于文献[7]的方法, SAE和MSE分别降低了16.50%和5.64%. 这是因为文献[7]设计的扰动观测器只能估计平稳变化的扰动. 在扰动频繁变化下, 扰动观测器的观测误差影响PID控制器的抗扰性能.
4. 物理实验
为验证所提方法的有效性, 以空气流量实验系统为对象进行物理实验. 实验系统如图13所示, 由空气流量系统, 监控计算机, 西门子S7-300 PLC控制系统, 检测仪表和执行机构组成. 以冷风机频率为控制输入, 冷风流量为输出, 令采样周期$ T_0 = 1\;{\rm{s}} $, 采集被控对象的输入输出数据, 辨识其ARX模型如下:
$$ \begin{split} y(k) = \;&0.720y(k-1)+0.872u(k-2)\;+ \\ &0.871u(k-3) \end{split} $$ (96) 设定值为:
$$ y_{sp}(t) = \begin{cases} 2,&\delta T<t\leq\delta T + \dfrac{T}{2} \\ -2,& \delta T + \dfrac{T}{2}<t\leq(\delta+1)T \end{cases} $$ (97) 设$ \delta = 0,1,2,3 $, $T = 5\;{\rm{min}}$.
针对实际物理系统, 将本文方法与PID控制器参数的离散Z-N法[18]及Astrom的PID参数优化方法[14]进行比较. 采用三种方法的PID控制器参数如表5所示. 三种方法的设定值$ y_{sp}(t) $、被控对象的输出(冷风流量)$ y(t) $及控制量(冷风机频率)$ u(t) $的曲线如图14 ~ 16所示. 可以看到, 本文方法下, 被控对象的输出对设定值的跟踪性能最好, 跟踪精度最高. 采用式(79)的SAE和式(80)的MSE对三种方法的评价结果如表6所示, 其中$N = 1\;200$. 本文方法相对于Z-N法, SAE和MSE分别降低76.33%和70.38%, 相对于Astrom方法SAE和MSE分别降低18.35%和9.76%, 验证了所提方法的优越性.
表 5 三种整定方法下PID控制器参数Table 5 PID controller parameters with three tuning methods整定方法 PID参数 $K_P$ $T_I$ $T_D$ Z-N法 0.3158 3.3412 0.8353 Astrom方法 0.1377 2.6910 1.1111 本文 0.0974 7.1364 1.7841 
图 14 Z-N法的设定值$y_{sp}(t)$与冷风流量$y(t)$及冷风机频率$u(t)$的曲线Fig. 14 Curves of set value $y_{sp}(t)$, cold air flow $y(t)$ and air cooler frequency $u(t)$ with Z-N method
图 15 Astrom优化方法的设定值$y_{sp}(t)$与冷风流量$y(t)$及冷风机频率$u(t)$的曲线Fig. 15 Curves of set value $y_{sp}(t)$, cold air flow $y(t)$ and air cooler frequency $u(t)$ with the optimization method proposed by Astrom
图 16 本文方法的设定值$y_{sp}(t)$与冷风流量$y(t)$及冷风机频率$u(t)$的曲线Fig. 16 Curves of set value $y_{sp}(t)$, cold air flow $y(t)$ and air cooler frequency $u(t)$ with the proposed method表 6 三种方法下误差$ e(t)$的性能评价Table 6 Performance evaluation of error $ e(t)$ with the three methods整定方法 性能指标 SAE MSE Z-N法 37624.0129 1736.8307 Astrom方法 10905.8089 570.0473 本文 8904.2220 514.4361 5. 结论
本文针对存在临界点的A类被控对象和不存在临界点的B类被控对象, 提出了基于$ -180^\circ $和$ -120^\circ $相位点的数字频率和增益的一种新的PID控制器参数优化整定方法. 该方法包括基于分数阶积分器和Tchebyshev多项式求取被控对象$ -180^\circ $或$ -120^\circ $相位点的频率和增益的方法和以阶跃负载扰动下跟踪误差平方和最小为优化指标, 以闭环系统的最大灵敏度和最大补灵敏度为约束的PID控制器参数优化整定规则. 仿真和物理实验表明所提PID控制器参数整定方法的优越性.
-
图 2 A类被控对象$ \rho_K $和$ \rho_T $的最优解曲线
Fig. 2 Optimal solution curve of $ \rho_K $ and $ \rho_T $ of class A controlled object
图 3 B类被控对象$ \rho_K $和$ \rho_T $的最优解曲线
Fig. 3 Optimal solution curve of $ \rho_K $ and $ \rho_T $ of class B controlled object
图 4 Z-N法的设定值$y_{sp}(t)$与输出$y(t)$及输入$u(t)$的曲线
Fig. 4 Curves of set value $y_{sp}(t)$, output $y(t)$ and input $u(t)$ with Z-N method
图 5 本文方法的设定值$y_{sp}(t)$与输出$y(t)$及输入$u(t)$的曲线
Fig. 5 Curves of set value $y_{sp}(t)$, output $y(t)$ and input $u(t)$ with the proposed method
图 6 被控对象(81)的Nyquist曲线$(K_{120} = 2.26,$ $\omega_{120} = 224.25)$
Fig. 6 Nyquist curve of controlled plant (81) $(K_{120} = 2.26,$ $\omega_{120} = 224.25)$
图 8 本文方法的设定值$y_{sp}(t)$与输出$y(t)$及输入$u(t)$的曲线
Fig. 8 Curves of set value $y_{sp}(t)$, output $y(t)$ and input $u(t)$ with the proposed method
图 10 本文方法的设定值$y_{sp}(t)$与输出$y(t)$及输入$u(t)$的曲线
Fig. 10 Curves of set value $y_{sp}(t)$, output $y(t)$ and input $u(t)$ with the proposed method
图 12 本文方法的设定值$y_{sp}(t)$与输出$y(t)$及输入$u(t)$的曲线
Fig. 12 Curves of set value $y_{sp}(t)$, output $y(t)$ and input $u(t)$ with the proposed method
图 14 Z-N法的设定值$y_{sp}(t)$与冷风流量$y(t)$及冷风机频率$u(t)$的曲线
Fig. 14 Curves of set value $y_{sp}(t)$, cold air flow $y(t)$ and air cooler frequency $u(t)$ with Z-N method
图 15 Astrom优化方法的设定值$y_{sp}(t)$与冷风流量$y(t)$及冷风机频率$u(t)$的曲线
Fig. 15 Curves of set value $y_{sp}(t)$, cold air flow $y(t)$ and air cooler frequency $u(t)$ with the optimization method proposed by Astrom
图 16 本文方法的设定值$y_{sp}(t)$与冷风流量$y(t)$及冷风机频率$u(t)$的曲线
Fig. 16 Curves of set value $y_{sp}(t)$, cold air flow $y(t)$ and air cooler frequency $u(t)$ with the proposed method
表 1 本文方法与Z-N法下误差$ e(t)$的性能评价
Table 1 Performance evaluation of error $ e(t)$ with the proposed method and Z-N method
整定方法 性能指标 SAE MSE 本文 3835.3081 0.0274 Z-N法 8004.8341 0.0479 
下载: 导出CSV
表 2 本文方法与文献[15]方法下误差$e(t)$的性能评价
Table 2 Performance evaluation of error $e(t)$ with the proposed method and the method in reference [15]
整定方法 性能指标 SAE MSE 本文 635.7247 0.0164 文献[15] 993.0357 0.0277
下载: 导出CSV
表 3 本文方法与文献[14]方法下误差$ e(t)$的性能评价
Table 3 Performance evaluation of error $ e(t)$ with the proposed method and the method in reference [14]
整定方法 性能指标 SAE MSE 本文 142.6175 0.0473 文献[14] 176.5287 0.0536
下载: 导出CSV
表 4 本文方法与文献[7]方法下误差$ e(t)$的性能评价
Table 4 Performance evaluation of error $ e(t)$ with the proposed method and the method in reference [7]
整定方法 性能指标 SAE MSE 本文 28889.3906 0.1806 文献[7] 34598.9258 0.1914
下载: 导出CSV
表 5 三种整定方法下PID控制器参数
Table 5 PID controller parameters with three tuning methods
整定方法 PID参数 $K_P$ $T_I$ $T_D$ Z-N法 0.3158 3.3412 0.8353 Astrom方法 0.1377 2.6910 1.1111 本文 0.0974 7.1364 1.7841
下载: 导出CSV
表 6 三种方法下误差$ e(t)$的性能评价
Table 6 Performance evaluation of error $ e(t)$ with the three methods
整定方法 性能指标 SAE MSE Z-N法 37624.0129 1736.8307 Astrom方法 10905.8089 570.0473 本文 8904.2220 514.4361
下载: 导出CSV
-
[1] Guo L. Feedback and uncertainty: Some basic problems and results. Annual Reviews in Control, 2020, 49: 27-36 doi: 10.1016/j.arcontrol.2020.04.001 [2] Samad T. A survey on industry impact and challenges thereof[Technical Activities]. IEEE Control Systems Magazine, 2017, 37(1): 17-18 doi: 10.1109/MCS.2016.2621438 [3] Borase R P, Maghade D K, Sondkar S Y, Pawar S N. A review of PID control, tuning methods and applications. International Journal of Dynamics and Control, 2021, 9(2): 818-827 doi: 10.1007/s40435-020-00665-4 [4] Somefun O A, Akingbade K, Dahunsi F. The dilemma of PID tuning. Annual Reviews in Control, 2022, 52: 65-74 [5] Ziegler J G, Nichols N B. Optimum settings for automatic controllers. Transactions of the ASME, 1942, 64: 759-768 [6] Åström K J, Hang C C, Persson P, Ho W K. Towards intelligent PID control. Automatica, 1992, 28(1): 1-9 doi: 10.1016/0005-1098(92)90002-W [7] Nie Z Y, Li Z Y, Wang Q G, Gao Z Q, Luo J L. A unifying Ziegler-nichols tuning method based on active disturbance rejection. International Journal of Robust and Nonlinear Control, 2022, 32(18): 9525-9541 doi: 10.1002/rnc.5848 [8] Åström K J, Hägglund T. Automatic tuning of simple regulators with specifications on phase and amplitude margins. Automatica, 1984, 20(5): 645-651 doi: 10.1016/0005-1098(84)90014-1 [9] Åström K J, Hägglund T. PID Controllers: Theory, Design, and Tuning (Second edition). Research Triangle Park: ISA, 1995. [10] Poulin E, Pomerleau A. PI settings for integrating processes based on ultimate cycle information. IEEE Transactions on Control Systems Technology, 1999, 7(4): 509-511 doi: 10.1109/87.772167 [11] Ossareh H R, Wisotzki S, Seeds J B, Jankovic M. An internal model control-based approach for characterization and controller tuning of turbocharged gasoline engines. IEEE Transactions on Control Systems Technology, 2021, 29(2): 866-875 doi: 10.1109/TCST.2019.2955917 [12] Huang H P, Jeng J C, Luo K Y. Auto-tune system using single-run relay feedback test and model-based controller design. Journal of Process Control, 2005, 15(6): 713-727 doi: 10.1016/j.jprocont.2004.11.004 [13] Mataušek M R, Šekara T B. PID controller frequency-domain tuning for stable, integrating and unstable processes, including dead-time. Journal of Process Control, 2011, 21(1): 17-27 doi: 10.1016/j.jprocont.2010.09.007 [14] Åström K J, Hägglund T. Advanced PID Control. Pittsburgh: ISA, 2006. [15] Bazanella A S, Pereira L F A, Parraga A. A new method for PID tuning including plants without ultimate frequency. IEEE Transactions on Control Systems Technology, 2017, 25(2): 637-644 doi: 10.1109/TCST.2016.2557723 [16] Lorenzini C, Bazanella A S, Pereira L F A, da Silva G R G. The generalized forced oscillation method for tuning PID controllers. ISA Transactions, 2019, 87: 68-87 doi: 10.1016/j.isatra.2018.11.014 [17] Takahashi Y, Chan C S, Auslander D M. Parametereinstellung bei linearen DDC-algorithmen. at - Automatisierungstechnik, 1971, 19(JG): 237-284 doi: 10.1524/auto.1971.19.jg.237 [18] Bobál V, Böhm J, Fessl J, Macháček J. Digital Self-tuning Controllers: Algorithms, Implementation and Applications. London: Springer, 2005. [19] Bazanella A S, Parraga A. Tuning PID controllers from sampled-data relay feedback experiments. IFAC-PapersOnLine, 2018, 51(4): 125-130 doi: 10.1016/j.ifacol.2018.06.106 [20] Rad A B, Lo W L, Tsang K M. Self-tuning PID controller using Newton-raphson search method. IEEE Transactions on Industrial Electronics, 1997, 44(5): 717-725 doi: 10.1109/41.633479 [21] Díaz-Rodríguez I D, Han S J, Bhattacharyya S P. Analytical Design of PID Controllers. Switzerland: Springer, 2019. [22] Chen Y Q, Vinagre B M, Podlubny I. Continued fraction expansion approaches to discretizing fractional order derivatives—an expository review. Nonlinear Dynamics, 2004, 38(1-4): 155-170 doi: 10.1007/s11071-004-3752-x [23] Zhang Q, Song B Y, Zhao H D, Zhang J S. Discretization of fractional order differentiator and integrator with different fractional orders. Intelligent Control and Automation, 2017, 8(2): 75-85 doi: 10.4236/ica.2017.82006 [24] 崔凤午, 梁怀学. 棣莫弗公式在推导倍角公式中的应用. 松辽学刊(自然科学版), 1997(4): 64-65Cui Feng-Wu, Liang Huai-Xue. The application of Demoivre Formula in the derivation of double Angle formula. Songliao Journal (Natural Science Edition), 1997(4): 64-65 [25] Huba M, Vrančič D. Introduction to the discrete time $ PID^m_n$ control for the IPDT plant. IFAC-PapersOnLine, 2018, 51(6): 119-124 doi: 10.1016/j.ifacol.2018.07.140 [26] Vítečcková M, Viteček A. 2DOF analog and digital controller tunning for integrating plants with time delay. Acta Mechanica Slovaca, 2009, 13(4): 48-53 doi: 10.2478/v10147-010-0036-y [27] Veinović S, Stojić D, Joksimović D. Optimized four-parameter PID controller for AVR systems with respect to robustness. International Journal of Electrical Power & Energy Systems, 2022, 135: Article No. 107529 [28] Wallén A, Årström K J, Hägglun T. Loop-shaping design of PID controllers with constant Ti/Td ratio. Asian Journal of Control, 2002, 4(4): 403-409 [29] 胡寿松. 自动控制原理. 第5版. 北京: 科学出版社, 2007.Hu Shou-Song. Automatic Control Theory (Fifth edition). Beijing: Science Press, 2007. [30] Solodov M V. On the sequential quadratically constrained quadratic programming methods. Mathematics of Operations Research, 2004, 29(1): 64-79 doi: 10.1287/moor.1030.0069 [31] Qi M, Zhang G P. An investigation of model selection criteria for neural network time series forecasting. European Journal of Operational Research, 2001, 132(3): 666-680 doi: 10.1016/S0377-2217(00)00171-5 [32] Haberstich C, Nouy A, Perrin G. Boosted optimal weighted least-squares. Mathematics of Computation, 2022, 91(335): 1281−1315 [33] Åström K J, Wittenmark B. Computer-controlled Systems: Theory and Design (Second edition). Englewood Cliffs: Prentice-Hall, 1990. [34] Hägglund T. A control-loop performance monitor. Control Engineering Practice, 1995, 3(11): 1543-1551 doi: 10.1016/0967-0661(95)00164-P [35] Ghanassi M, Champagne B, Kabal P. On the steady-state mean squared error of the fixed-point LMS algorithm. Signal Processing, 2007, 87(12): 3226-3233 doi: 10.1016/j.sigpro.2007.05.029 [36] Euzébio T A M, Barros P R. Optimal integral gain for smooth PI control. IFAC Proceedings Volumes, 2013, 46(11): 529-533 doi: 10.3182/20130703-3-FR-4038.00125 期刊类型引用(12)
1. 何立新，曹辰宇，张峥，雷晓辉，李翔. 引黄济青明渠段输水控制系统的MIL测试系统设计与实现. 南水北调与水利科技(中英文). 2025(01): 1-9+58 . 
百度学术2. 陈家煜，刘丽桑，张友渊，陈炯晖，王晨曦. 基于改进PID算法的钢铁行车系统建模与控制. 福建理工大学学报. 2024(01): 47-53 . 百度学术3. 黄鑫，谢董玉. 基于模糊控制的永磁同步电机协同控制. 化工自动化及仪表. 2024(02): 168-172 . 百度学术4. 丁勇军，王建宏，罗熙，张金龙. 紧格式无模型自适应控制在四旋翼飞行器中的应用. 电光与控制. 2024(06): 87-93 . 百度学术5. 舒奕彬，李立君，张振翮，戚浩，刘姜毅. 基于改进海马优化算法的PID参数优化. 机床与液压. 2024(13): 189-194 . 百度学术6. 周毓金，孙文彩，戴山力，张磊. 船用蒸汽排放系统动态特性分析及控制研究. 舰船电子工程. 2024(05): 124-132 . 百度学术7. 杨婷，安娜，王昊，祝小虎，马海军，翟杰，句瑞婷. 基于PID自控系统的原油稳定工艺关键设备参数动态模拟. 油气田地面工程. 2024(07): 29-34 . 百度学术8. 安世硕，高海，李维军，唐静，张文才. 基于改进鲸鱼算法的收卷张力PID控制的性能优化. 辽宁石油化工大学学报. 2024(04): 91-96 . 百度学术9. 张文林. 高效节能泵站的反步控制研究. 自动化应用. 2024(14): 56-58+61 . 百度学术10. 郑德红，殷利平. 基于扰动观测器的输入饱和非线性系统运行优化控制. 科学技术与工程. 2024(35): 15082-15089 . 百度学术11. 夏知胜，陈诚，杨爱斌. 基于遗传算法优化PID的固定翼无人机俯仰控制设计. 科技创新与应用. 2023(34): 41-44 . 百度学术12. 张宇，赵胜雪，邓炜航. 玉米电控排种系统的设计与试验. 中国农机装备. 2023(12): 45-49 . 百度学术其他类型引用(49)
-
下载:
计量
- 文章访问数: 2263
- HTML全文浏览量: 1001
- PDF下载量: 740
- 被引次数: 61
