CN105159097A - 炼油加热炉炉膛压力的多变量预测控制pid控制方法 - Google Patents

Abstract

本发明公开了一种炼油加热炉炉膛压力的多变量预测控制PID控制方法。本发明方法首先基于炼油加热炉炉膛压力对象的实时多输入多输出数据建立炼油加热炉炉膛压力对象的状态空间模型，然后结合状态过程和输出误差建立扩展的非最小状态空间模型。在此模型的基础上，依据预测函数控制的方法来优化PID控制器的参数，最后对被控对象实现PID控制。本发明弥补了传统控制方式的不足，保证具备良好控制性能，同时又具备PID控制的简单形式，在对多变量也具有良好的控制效果。

Description

炼油加热炉炉膛压力的多变量预测控制PID控制方法

技术领域

本发明属于自动化技术领域，涉及一种基于预测函数控制优化的炼油加热炉炉膛压力比例积分微分(PID)控制方法，具体的是一种炼油加热炉炉膛压力的多变量预测控制PID控制方法。

背景技术

化工过程是我国流程工业过程的重要组成部分，其要求是供给合格的工业产品，以满足我国工业的需要。在目前工业控制中，因为传统的PID控制技术具有结构简单、鲁棒性好、易于操作等优点，所以在工业中得到了广泛的应用。炼油加热炉炉膛是一个典型的多输入多输出(MIMO)对象，传统的PID控制理论针对多变量对象的控制参数整定过程比较复杂，计算量比较大，控制性能差，往往不能满足实际化工过程中日益严格的控制精度和产品质量要求，更加先进、控制效果更好的算法仍然有待研究。针对多变量过程的先进控制方法有很多，现已取得了一定的成果。预测函数控制(PFC)是其中的一种先进控制方法，具有实时性好，输入规律明确且能够有效地减少算法的计算量等优点。如果能够将PFC控制性能好和PID控制结构简单的特性进行结合，使实际的多变量控制系统同时具备PFC和PID二者的优良特点，这样既能够保证控制结构的形式简单，又能够获得更好的控制效果。

发明内容

本发明的目的是针对现有的化工过程系统控制技术的不足之处，提供一种炼油加热炉炉膛压力的多变量预测控制PID控制方法。该方法弥补了传统控制方式的不足，保证具备良好控制性能，同时又具备PID控制的简单形式，在对多变量也具有良好的控制效果。

本发明方法首先基于炼油加热炉炉膛压力对象的实时多输入多输出数据建立炼油加热炉炉膛压力对象的状态空间模型，然后结合状态过程和输出误差建立扩展的非最小状态空间模型。在此模型的基础上，依据预测函数控制的方法来优化PID控制器的参数，最后对被控对象实现PID控制。

本发明的技术方案是通过数据采集、模型建立、优化等手段，提供一种炼油加热炉炉膛压力的多变量预测控制PID控制方法，从而在实际过程中对多变量获得更好地控制性能，有效提高系统的控制性能。

本发明方法的步骤包括：

步骤1、建立被控对象的扩展非最小状态空间模型，具体步骤是：1.1通过采集被控对象的实时输入输出数据，利用最小二乘法建立多输入多输出系统模型。形式如下：

Y(k+1)+L₁Y(k)+L₂Y(k-1)+……+L_nY(k-n+1)＝S₁U(k)+S₂U(k-1)+……+S_nU(k-n+1)

其中，Y(k)表示k时刻预测系统模型的q维输出值，U(k)表示k时刻p维输入，L₁,L₂,……，L_n为需要辨识的标量系数，S₁,S₂,……,S_n为待辨识的q×p矩阵。

$Y\left(k\right)=\left[\begin{array}{c}{y}_1\left(k\right)\\ y_2\left(k\right)\\ .\\ .\\ .\\ y_q\left(k\right)\end{array}\right],U\left(k\right)=\left[\begin{array}{c}{u}_1\left(k\right)\\ u_2\left(k\right)\\ .\\ .\\ .\\ u_p\left(k\right)\end{array}\right]$

$S_i=\left[\begin{array}{cccc}{b}_{1i1}& b_{1i2}& \mathrm{...}& b_{1ip}\\ b_{2i1}& b_{2i2}& \mathrm{...}& b_{2ip}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ b_{qi1}& b_{qi2}& \mathrm{...}& b_{qip}\end{array}\right],i=0,1,...,n$

如果将S_i所有参数同时进行辨识，则计算量很大。所以采用一行一行的进行辨识，利用采集的实时过程数据，得到N组样本数据，形式如下：

$Y_j\left(k\right)=H_j{\mathrm{\θ}}_j,{\mathrm{\θ}}_j^T=\[L_1...L_n{b}_{j01}...b_{j0p}{b}_{j11}...b_{j1p}...b_{jn1}...b_{jnp}\]$

$H_j=\left[\begin{array}{cccccc}-y_j\left(n\right)& \mathrm{...}& -y_j\left(1\right)& U^T\left(n+1\right)& \mathrm{...}& U^T\left(1\right)\\ -y_j\left(n+1\right)& \mathrm{...}& -y_j\left(2\right)& U^T\left(n+2\right)& \mathrm{...}& U^T\left(2\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ -y_j\left(n+N-1\right)& \mathrm{...}& -y_j\left(N\right)& U^T\left(n+N\right)& \mathrm{...}& U^T\left(N\right)\end{array}\right]$

其中，H_j、y(j)分别表示采集的第j组的输入数据和输出值，N表示样本总数。

辨识结果为：

${\widehat{\mathrm{\θ}}}_j={\left({H_j}^T{H}_j\right)}^{-1}{H_j}^T{Y}_j\left(k\right)$

按照上式，令j＝1,2,…，j-1，j+1，…，q可得其他各行的参数估计值，在求其他行的参数时L₁,L₂,…L_n不必再估计，把这些值代入以减少其他各行的计算量。

1.2将1.1步骤中得到的系统模型转换成差分模型形式：

△y(k+1)+L₁△y(k)+L₂△y(k-1)+…+L_n△y(k-n+1)＝S₁△u(k)+S₂△u(k-1)+…+S_n△u(k-n+1)

其中，△是差分算子。

1.3选取如下所示的非最小状态空间变量△x_m(k)：

△x_m(k)^T＝[△y(k)^T,△y(k-1)^T,…,△y(k-n+1)^T,△u(k-1)^T,

△u(k-2)^T,…,△u(k-n+1)^T]

进而将步骤1.2中的差分模型转化为状态空间模型，其形式如下：

△x_m(k+1)＝A_m△x_m(k)+B_m△u(k)

△y(k+1)＝C_m△x_m(k+1)

其中，

$A_m=\left[\begin{array}{ccccccccc}-L_1& -L_2& \mathrm{...}& -L_{n-1}& -L_n& S_2& \mathrm{...}& S_{n-1}& S_n\\ I_q& 0& \mathrm{...}& 0& 0& 0& \mathrm{...}& 0& 0\\ 0& I_q& \mathrm{...}& 0& 0& 0& \mathrm{...}& 0& 0\\ .& .& & .& .& .& & .& .\\ .& .& \mathrm{...}& .& .& .& \mathrm{...}& .& .\\ .& .& & .& .& .& & .& .\\ 0& 0& \mathrm{...}& I_q& 0& 0& \mathrm{...}& 0& 0\\ 0& 0& \mathrm{...}& 0& 0& 0& \mathrm{...}& 0& 0\\ 0& 0& \mathrm{...}& 0& 0& I_p& \mathrm{...}& 0& 0\\ .& .& & .& .& & .& .& .\\ .& .& \mathrm{...}& .& .& \mathrm{...}& .& .& .\\ .& .& & .& .& & .& .& .\\ 0& 0& \mathrm{...}& 0& 0& 0& \mathrm{...}& I_p& 0\end{array}\right]$

B_m＝[S₁0…0I_p0…0]^Τ

C_m＝[I_q00…0000]

△x_m(k)的维数m＝p×(n-1)+q×n。

1.4选取新的状态变量

$z\left(k\right)=\left[\begin{array}{c}\mathrm{\Δ}{x}_m\left(k\right)\\ y\left(k\right)\end{array}\right]$

然后将步骤1.3中得到的状态空间模型转换成包含状态变量和输出值的扩展非最小状态空间模型，形式如下：

z(k+1)＝Az(k)+B△u(k)

y(k+1)＝Cz(k+1)

e(k)＝y(k)-r(k)

其中， $A=\left[\begin{array}{cc}{A}_m& 0\\ C_m{A}_m& I_q\end{array}\right],$ $B=\left[\begin{array}{c}{B}_m\\ C_m{B}_m\end{array}\right],$ C＝[0I_q]；

0是维数为m×q零矩阵，I_q是维数为q单位矩阵，I_p是维数为p单位矩阵，r(k)是k时刻的期望输出值，e(k)为k时刻的实际输出值与期望输出值之间的差值。那么e(k+1)式子将进一步表示成：

e(k+1)＝e(k)+C_mA_m△x_m(k)+C_mB_m△u(k)-△r(k+1)

这时在选取一个新的状态变量 $z_e\left(k\right)=\left[\begin{array}{c}\mathrm{\Δ}{x}_m\left(k\right)\\ e\left(k\right)\end{array}\right]$

将上式转换成包含状态变量和输出误差的扩展非最小状态空间模型，被表示为如下：

z_e(k+1)＝A_ez(k)+B_e△u(k)+C_e△r(k+1)

其中，

$A_e=\left[\begin{array}{cc}{A}_m& 0\\ C_m{A}_m& I_q\end{array}\right];B_e=\left[\begin{array}{c}{B}_m\\ C_m{B}_m\end{array}\right];C_e=\left[\begin{array}{c}0\\ -I_q\end{array}\right]$

0是维数为m×q的零矩阵，I_q是维数为q的单位矩阵。

步骤2、设计被控对象的PID控制器，具体步骤是：

2.1计算k时刻对第k+P时刻的预测输出值，形式如下：

其中，

$\mathrm{\ζ}=\left[\begin{array}{cccc}{A_e}^{P-1}{C}_e& {A_e}^{P-2}{C}_e& \mathrm{...}& C_e\end{array}\right],\mathrm{\γ}={A_e}^{P-1}{B}_e$

△R＝[△r(k+1)△r(k+2)…△r(k+P)]^Τ

r(k+i)＝αⁱy(k)+(1-αⁱ)c(k),i＝1,2,…,P；

P为预测时域，A^P表示P个A矩阵相乘，α是参考轨迹的柔化因子，c(k)是k时刻的设定值。

2.2选取被控对象的目标函数J(k)，形式如下：

minJ(k)＝z_e(k+P)^ΤQ_ez_e(k+P)

其中，Q_e是(2n-1)×(2n-1)权矩阵，min表示求最小值。

2.3根据步骤2.2中的目标函数求解PID控制器的参数，具体方法是：先将控制量u(k)进行变换：

u(k)＝u(k-1)+K_p(k)(e₁(k)-e₁(k-1))+K_i(k)e₁(k)+K_d(k)(e₁(k)-2e₁(k-1)+e₁(k-2))

e₁(k)＝[e₁₁(k)e₁₂(k)…e_1q(k)]^T

其中，K_p(k)、K_i(k)、K_d(k)分别是k时刻PID控制器的比例、微分、积分参数，e₁(k)是k时刻设定值与实际输出值之间的误差。

进而将控制量u(k)简化成矩阵形式：

u(k)＝u(k-1)+E(k)^Tw(k)

其中，

$E\left(k\right)={\left[\begin{array}{ccccc}{E}_1\left(k\right)& 0& 0& \mathrm{...}& 0\\ 0& E_2\left(k\right)& 0& \mathrm{...}& 0\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ 0& \mathrm{...}& 0& E_{p-1}(k-1)& 0\\ 0& 0& \mathrm{...}& 0& E_p\left(k\right)\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\end{array}\right]}_{3q\×p}$

E_i(k)＝[e_1i(k),e_1i(k-1),e_1i(k-2)]^Τ

w_i(k)＝[w_i1(k),w_i2(k),w_i3(k)]

w(k)＝[w₁(k)w₂(k)…w_q(k)]^T

w_i1(k)＝K_pi(k)+K_ii(k)+K_di(k)

w_i2(k)＝-K_pi(k)-2K_di(k)

w_i3(k)＝K_di(k)

结合控制量u(k)的矩阵形式和步骤2.2中的目标函数，求得：

w(k)＝E(k)(-(γ^TQ_eγE(k)^TE(k))^-1γ^TQ_e(A_e ^Pz_e(k)+ζ△R))

进一步得到：

K_pi(k)＝-w_i2(k)-2K_di(k)

K_ii(k)＝w_i1(k)-K_pi(k)-K_di(k)

K_di(k)＝w_i3(k)

2.4得到PID控制器的参数K_p(k)、K_i(k)、K_d(k)后，构成控制量：u(k)＝u(k-1)+K_p(k)(e₁(k)-e₁(k-1))+K_i(k)e₁(k)+K_d(k)(e₁(k)-2e₁(k-1)+e₁(k-2))，再将u(k)作用于被控对象。

2.5在k+l时刻，依照2.1到2.4中的步骤循环求解PID控制器新的参数K_p(k+l)、K_i(k+l)、K_d(k+l)，l＝1,2,3，…。当满足如下条件时结束循环：

$\left\{\begin{array}{c}{k}_{pi}\left(k\right)=k_{pi}\left(k-1\right)\\ k_{ii}\left(k\right)=k_{ii}\left(k-1\right)\mathrm{......}\left|e_{1i}\left(k\right)\right|\≤\mathrm{\ϵ}\\ k_{di}\left(k\right)=k_{di}\left(k-1\right)\end{array}\right.$

$\left\{\begin{array}{c}{k}_{pi}\left(k\right)=-w_{i2}\left(k\right)-2k_{di}\left(k\right)\\ k_{ii}=w_{i1}\left(k\right)-k_{pi}\left(k\right)-k_{di}\left(k\right)\mathrm{......}\left|e_{1i}\left(k\right)\right|>\mathrm{\ϵ}\\ k_{di}\left(k\right)=w_{i3}\left(k\right)\end{array}\right..$

本发明有益效果：

本发明方法首先基于炼油加热炉炉膛压力对象的实时多输入多输出数据建立炼油加热炉炉膛压力对象的状态空间模型，然后结合状态过程和输出误差建立扩展的非最小状态空间模型。在此模型的基础上，依据预测函数控制的方法来优化PID控制器的参数，最后对被控对象实现PID控制，有效地弥补了传统控制方法的不足，可有效提高系统的控制性能。

具体实施方式

以炼油加热炉炉膛压力过程控制为例：

在炼油加热炉炉膛压力控制过程中，调节手段为调节烟道挡板的开度。

步骤1、建立被控对象的扩展非最小状态空间模型，具体步骤是：

1.1通过采集被控对象的实时输入输出数据，利用最小二乘法建立多输入多输出系统模型。形式如下：

Y(k+1)+L₁Y(k)+L₂Y(k-1)+……+L_nY(k-n+1)＝S₁U(k)+S₂U(k-1)+……+S_nU(k-n+1)

其中，Y(k)表示k时刻预测系统模型的q维输出值，U(k)表示k时刻p维输入，L₁,L₂,……，L_n为需要辨识的标量系数，S₁,S₂,……,S_n为待辨识的q×p矩阵。

$Y\left(k\right)=\left[\begin{array}{c}{y}_1\left(k\right)\\ y_2\left(k\right)\\ .\\ .\\ .\\ y_q\left(k\right)\end{array}\right],U\left(k\right)=\left[\begin{array}{c}{u}_1\left(k\right)\\ u_2\left(k\right)\\ .\\ .\\ .\\ u_p\left(k\right)\end{array}\right]$

$S_i=\left[\begin{array}{cccc}{b}_{1i1}& b_{1i2}& \mathrm{...}& b_{1ip}\\ b_{2i1}& b_{2i2}& \mathrm{...}& b_{2ip}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ b_{qi1}& b_{qi2}& \mathrm{...}& b_{qip}\end{array}\right],i=0,1,...,n$

如果将S_i所有参数同时进行辨识，则计算量很大。所以采用一行一行的进行辨识，利用采集的实时过程数据，得到N组样本数据，形式如下：

$Y_j\left(k\right)=H_j{\mathrm{\θ}}_j,{\mathrm{\θ}}_j^T=\[L_1...L_n{b}_{j01}...b_{j0p}{b}_{j11}...b_{j1p}...b_{jn1}...b_{jnp}\]$

$H_j=\left[\begin{array}{cccccc}-y_j\left(n\right)& \mathrm{...}& -y_j\left(1\right)& U^T\left(n+1\right)& \mathrm{...}& U^T\left(1\right)\\ -y_j\left(n+1\right)& \mathrm{...}& -y_j\left(2\right)& U^T\left(n+2\right)& \mathrm{...}& U^T\left(2\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ -y_j\left(n+N-1\right)& \mathrm{...}& -y_j\left(N\right)& U^T\left(n+N\right)& \mathrm{...}& U^T\left(N\right)\end{array}\right]$

其中，H_j、y(j)分别表示采集的第j组的输入数据和输出值，N表示样本总数。

辨识结果为：

${\widehat{\mathrm{\θ}}}_j={\left({H_j}^H{H}_j\right)}^{-1}{H_j}^T{Y}_j\left(k\right)$

按照上式，令j＝1,2,…，j-1，j+1，…，q可得其他各行的参数估计值，在求其他行的参数时L₁,L₂,…L_n不必再估计，把这些值代入以减少其他各行的计算量。

1.2将1.1步骤中得到的系统模型转换成差分模型形式：

△y(k+1)+L₁△y(k)+L₂△y(k-1)+…+L_n△y(k-n+1)＝S₁△u(k)+S₂△u(k-1)+…+S_n△u(k-n+1)

其中，△是差分算子。

1.3选取如下所示的非最小状态空间变量△x_m(k)：

△x_m(k)^T＝[△y(k)^T,△y(k-1)^T,…,△y(k-n+1)^T,△u(k-1)^T,

△u(k-2)^T,…,△u(k-n+1)^T]

进而将步骤1.2中的差分模型转化为状态空间模型，其形式如下：

△x_m(k+1)＝A_m△x_m(k)+B_m△u(k)

△y(k+1)＝C_m△x_m(k+1)

其中，

$A_m=\left[\begin{array}{ccccccccc}-L_1& -L_2& \mathrm{...}& -L_{n-1}& -L_n& S_2& \mathrm{...}& S_{n-1}& S_n\\ I_q& 0& \mathrm{...}& 0& 0& 0& \mathrm{...}& 0& 0\\ 0& I_q& \mathrm{...}& 0& 0& 0& \mathrm{...}& 0& 0\\ .& .& & .& .& .& & .& .\\ .& .& \mathrm{...}& .& .& .& \mathrm{...}& .& .\\ .& .& & .& .& .& & .& .\\ 0& 0& \mathrm{...}& I_q& 0& 0& \mathrm{...}& 0& 0\\ 0& 0& \mathrm{...}& 0& 0& 0& \mathrm{...}& 0& 0\\ 0& 0& \mathrm{...}& 0& 0& I_p& \mathrm{...}& 0& 0\\ .& .& & .& .& & .& .& .\\ .& .& \mathrm{...}& .& .& \mathrm{...}& .& .& .\\ .& .& & .& .& & .& .& .\\ 0& 0& \mathrm{...}& 0& 0& 0& \mathrm{...}& I_p& 0\end{array}\right]$

B_m＝[S₁0…0I_p0…0]^Τ

C_m＝[I_q00…0000]

△x_m(k)的维数m＝p×(n-1)+q×n。

1.4选取新的状态变量

$z\left(k\right)=\left[\begin{array}{c}\mathrm{\Δ}{x}_m\left(k\right)\\ y\left(k\right)\end{array}\right]$

然后将步骤1.3中得到的状态空间模型转换成包含状态变量和输出值的扩展非最小状态空间模型，形式如下：

z(k+1)＝Az(k)+B△u(k)

y(k+1)＝Cz(k+1)

e(k)＝y(k)-r(k)

其中， $A=\left[\begin{array}{cc}{A}_m& 0\\ C_m{A}_m& I_q\end{array}\right],$ $B=\left[\begin{array}{c}{B}_m\\ C_m{B}_m\end{array}\right],$ C＝[0I_q]；

0是维数为m×q零矩阵，I_q是维数为q单位矩阵，I_p是维数为p单位矩阵，r(k)是k时刻的期望输出值，e(k)为k时刻的实际输出值与期望输出值之间的差值。那么e(k+1)式子将进一步表示成：

e(k+1)＝e(k)+C_mA_m△x_m(k)+C_mB_m△u(k)-△r(k+1)

这时在选取一个新的状态变量 $z_e\left(k\right)=\left[\begin{array}{c}\mathrm{\Δ}{x}_m\left(k\right)\\ e\left(k\right)\end{array}\right]$

将上式转换成包含状态变量和输出误差的扩展非最小状态空间模型，被表示为如下：

z_e(k+1)＝A_ez(k)+B_e△u(k)+C_e△r(k+1)

其中，

$A_e=\left[\begin{array}{cc}{A}_m& 0\\ C_m{A}_m& I_q\end{array}\right];B_e=\left[\begin{array}{c}{B}_m\\ C_m{B}_m\end{array}\right];C_e=\left[\begin{array}{c}0\\ -I_q\end{array}\right]$

0是维数为m×q的零矩阵，I_q是维数为q的单位矩阵。

步骤2、设计被控对象的PID控制器，具体步骤是：

2.1计算k时刻对第k+P时刻的预测输出值，形式如下：

其中，

$\mathrm{\ζ}=\left[\begin{array}{cccc}{A_e}^{P-1}{C}_e& {A_e}^{P-2}{C}_e& \mathrm{...}& C_e\end{array}\right],\mathrm{\γ}={A_e}^{P-1}{B}_e$

△R＝[△r(k+1)△r(k+2)…△r(k+P)]^Τ

r(k+i)＝αⁱy(k)+(1-αⁱ)c(k),i＝1,2,…,P；

P为预测时域，A^P表示P个A矩阵相乘，α是参考轨迹的柔化因子，c(k)是k时刻的设定值。

2.2选取被控对象的目标函数J(k)，形式如下：

minJ(k)＝z_e(k+P)^ΤQ_ez_e(k+P)

其中，Q_e是(2n-1)×(2n-1)权矩阵，min表示求最小值。

2.3根据步骤2.2中的目标函数求解PID控制器的参数，具体方法是：先将控制量u(k)进行变换：

u(k)＝u(k-1)+K_p(k)(e₁(k)-e₁(k-1))+K_i(k)e₁(k)+K_d(k)(e₁(k)-2e₁(k-1)+e₁(k-2))

e₁(k)＝[e₁₁(k)e₁₂(k)…e_1q(k)]^T

其中，K_p(k)、K_i(k)、K_d(k)分别是k时刻PID控制器的比例、微分、积分参数，e₁(k)是k时刻设定值与实际输出值之间的误差。

进而将控制量u(k)简化成矩阵形式：

u(k)＝u(k-1)+E(k)^Tw(k)

其中，

$E\left(k\right)={\left[\begin{array}{ccccc}{E}_1\left(k\right)& 0& 0& \mathrm{...}& 0\\ 0& E_2\left(k\right)& 0& \mathrm{...}& 0\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ 0& \mathrm{...}& 0& E_{p-1}(k-1)& 0\\ 0& 0& \mathrm{...}& 0& E_p\left(k\right)\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\end{array}\right]}_{3q\×p}$

E_i(k)＝[e_1i(k),e_1i(k-1),e_1i(k-2)]^Τ

w_i(k)＝[w_i1(k),w_i2(k),w_i3(k)]

w(k)＝[w₁(k)w₂(k)…w_q(k)]^T

w_i1(k)＝K_pi(k)+K_ii(k)+K_di(k)

w_i2(k)＝-K_pi(k)-2K_di(k)

w_i3(k)＝K_di(k)

结合控制量u(k)的矩阵形式和步骤2.2中的目标函数，求得：

w(k)＝E(k)(-(γ^TQ_eγE(k)^TE(k))^-1γ^TQ_e(A_e ^Pz_e(k)+ζ△R))

进一步得到：

K_pi(k)＝-w_i2(k)-2K_di(k)

K_ii(k)＝w_i1(k)-K_pi(k)-K_di(k)

K_di(k)＝w_i3(k)

2.4得到PID控制器的参数K_p(k)、K_i(k)、K_d(k)后，构成控制量：u(k)＝u(k-1)+K_p(k)(e₁(k)-e₁(k-1))+K_i(k)e₁(k)+K_d(k)(e₁(k)-2e₁(k-1)+e₁(k-2))，再将u(k)作用于被控对象。

2.5在k+l时刻，依照2.1到2.4中的步骤循环求解PID控制器新的参数K_p(k+l)、K_i(k+l)、K_d(k+l)，l＝1,2,3，…。当满足如下条件时结束循环：

$\left\{\begin{array}{c}{k}_{pi}\left(k\right)=k_{pi}\left(k-1\right)\\ k_{ii}\left(k\right)=k_{ii}\left(k-1\right)\mathrm{......}\left|e_{1i}\left(k\right)\right|\≤\mathrm{\ϵ}\\ k_{di}\left(k\right)=k_{di}\left(k-1\right)\end{array}\right.$

$\left\{\begin{array}{c}{k}_{pi}\left(k\right)=-w_{i2}\left(k\right)-2k_{di}\left(k\right)\\ k_{ii}=w_{i1}\left(k\right)-k_{pi}\left(k\right)-k_{di}\left(k\right)\mathrm{......}\left|e_{1i}\left(k\right)\right|>\mathrm{\ϵ}\\ k_{di}\left(k\right)=w_{i3}\left(k\right)\end{array}\right..$

Claims (3)

1.炼油加热炉炉膛压力的多变量预测控制PID控制方法，其特征在于包括如下步骤：

步骤1、建立被控对象的扩展非最小状态空间模型；

步骤2、设计被控对象的PID控制器。

2.如权利要求1所述的炼油加热炉炉膛压力的多变量预测控制PID控制方法，其特征在于步骤1所述的建立被控对象的扩展非最小状态空间模型，具体步骤如下：

1.1通过采集被控对象的实时输入输出数据，利用最小二乘法建立多输入多输出系统模型；形式如下：

Y(k+1)+L₁Y(k)+L₂Y(k-1)+……+L_nY(k-n+1)＝

S₁U(k)+S₂U(k-1)+……+S_nU(k-n+1)

其中，Y(k)表示k时刻预测系统模型的q维输出值，U(k)表示k时刻p维输入，L₁,L₂,……，L_n为需要辨识的标量系数，S₁,S₂,……,S_n为待辨识的q×p矩阵；

$Y\left(k\right)=\left[\begin{array}{c}{y}_1\left(k\right)\\ y_2\left(k\right)\\ .\\ .\\ .\\ y_q\left(k\right)\end{array}\right],U\left(k\right)=\left[\begin{array}{c}{u}_1\left(k\right)\\ u_2\left(k\right)\\ .\\ .\\ .\\ u_p\left(k\right)\end{array}\right]$

$S_i=\left[\begin{array}{cccc}{b}_{1i1}& b_{1i2}& ...& b_{1ip}\\ b_{2i1}& b_{2i2}& ...& b_{2ip}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ b_{qi1}& b_{qi2}& ...& b_{qip}\end{array}\right],i=0,1,...,n$

采用一行一行的进行辨识，利用采集的实时过程数据，得到N组样本数据，形式如下：

$Y_j\left(k\right)=H_j{\mathrm{\θ}}_j,{\mathrm{\θ}}_j^T=\[L_1...L_n{b}_{j01}...b_{j0p}{b}_{j11}...b_{j1p}...b_{jn1}...b_{jnp}\]$

$H_j=\left[\begin{array}{cccccc}-y_j\left(n\right)& ...& -y_j\left(1\right)& U^T\left(n+1\right)& ...& U^T\left(1\right)\\ -y_j\left(n+1\right)& ...& -y_j\left(2\right)& U^T\left(n+2\right)& ...& U^T\left(2\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ -y_j\left(n+N-1\right)& ...& -y_j\left(N\right)& U^T\left(n+N\right)& ...& U^T\left(N\right)\end{array}\right]$

其中，H_j、y(j)分别表示采集的第j组的输入数据和输出值，N表示样本总数，辨识结果为：

${\widehat{\mathrm{\θ}}}_j={\left({H_j}^T{H}_j\right)}^{-1}{H_j}^T{Y}_j\left(k\right)$

按照上式，令j＝1,2,…，j-1，j+1，…，q可得其他各行的参数估计值，在求其他行的参数时L₁,L₂,…L_n不必再估计，把这些值代入以减少其他各行的计算量；

1.2将1.1步骤中得到的系统模型转换成差分模型形式：

△y(k+1)+L₁△y(k)+L₂△y(k-1)+…+L_n△y(k-n+1)

＝S₁△u(k)+S₂△u(k-1)+…+S_n△u(k-n+1)

其中，△是差分算子；

1.3选取如下所示的非最小状态空间变量△x_m(k)：

△x_m(k)^T＝[△y(k)^T,△y(k-1)^T,…,△y(k-n+1)^T,△u(k-1)^T,

△u(k-2)^T,…,△u(k-n+1)^T]

进而将步骤1.2中的差分模型转化为状态空间模型，其形式如下：

△x_m(k+1)＝A_m△x_m(k)+B_m△u(k)

△y(k+1)＝C_m△x_m(k+1)

其中，

$A_m=\left[\begin{array}{ccccccccc}-L_1& -L_2& ...& -L_{n-1}& -L_n& S_2& ...& S_{n-1}& S_n\\ I_q& 0& ...& 0& 0& 0& ...& 0& 0\\ 0& I_q& ...& 0& 0& 0& ...& 0& 0\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ ...\\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ ...\\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ 0& 0& ...& I_q& 0& 0& ...& 0& 0\\ 0& 0& ...& 0& 0& 0& ...& 0& 0\\ 0& 0& ...& 0& 0& I_p& ...& 0& 0\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ ...\\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ ...\\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ 0& 0& ...& 0& 0& 0& ...& I_p& 0\end{array}\right]$

B_m＝[S₁0…0I_p0…0]^Τ

C_m＝[I_q00…0000]

△x_m(k)的维数m＝p×(n-1)+q×n；

1.4选取新的状态变量

$z\left(k\right)=\left[\begin{array}{c}{\mathrm{\Δx}}_m\left(k\right)\\ y\left(k\right)\end{array}\right]$

然后将步骤1.3中得到的状态空间模型转换成包含状态变量和输出值的扩展非最小状态空间模型，形式如下：

z(k+1)＝Az(k)+B△u(k)

y(k+1)＝Cz(k+1)

e(k)＝y(k)-r(k)

其中， $A=\left[\begin{array}{cc}{A}_m& 0\\ C_m{A}_m& I_q\end{array}\right],B=\left[\begin{array}{c}{B}_m\\ C_m{B}_m\end{array}\right],$ C＝[0I_q]；

0是维数为m×q零矩阵，I_q是维数为q单位矩阵，I_p是维数为p单位矩阵，r(k)是k时刻的期望输出值，e(k)为k时刻的实际输出值与期望输出值之间的差值；那么e(k+1)式子将进一步表示成：

e(k+1)＝e(k)+C_mA_m△x_m(k)+C_mB_m△u(k)-△r(k+1)

此时选取一个新的状态变量 $z_e\left(k\right)=\left[\begin{array}{c}{\mathrm{\Δx}}_m\left(k\right)\\ e\left(k\right)\end{array}\right]$

将上式转换成包含状态变量和输出误差的扩展非最小状态空间模型，被表示为如下：

z_e(k+1)＝A_ez(k)+B_e△u(k)+C_e△r(k+1)

其中，

$A_e=\left[\begin{array}{cc}{A}_m& 0\\ C_m{A}_m& I_q\end{array}\right];B_e=\left[\begin{array}{c}{B}_m\\ C_m{B}_m\end{array}\right];C_e=\left[\begin{array}{c}0\\ -I_q\end{array}\right]$

0是维数为m×q的零矩阵，I_q是维数为q的单位矩阵。

3.根据权利要求2所述的炼油加热炉炉膛压力的多变量预测控制PID控制方法，其特征在于步骤2所述的设计被控对象的PID控制器，具体步骤如下：

2.1计算k时刻对第k+P时刻的预测输出值，形式如下：

其中，

$\mathrm{\ζ}=\left[\begin{array}{cccc}{A_e}^{P-1}{C}_e& {A_e}^{P-2}{C}_e& ...& C_e\end{array}\right],\mathrm{\γ}={A_e}^{P-1}{B}_e$

△R＝[△r(k+1)△r(k+2)…△r(k+P)]^Τ

r(k+i)＝αⁱy(k)+(1-αⁱ)c(k),i＝1,2,…,P；

P为预测时域，A^P表示P个A矩阵相乘，α是参考轨迹的柔化因子，c(k)是k时刻的设定值；

2.2选取被控对象的目标函数J(k)，形式如下：

minJ(k)＝z_e(k+P)^ΤQ_ez_e(k+P)

其中，Q_e是(2n-1)×(2n-1)权矩阵，min表示求最小值；

2.3根据步骤2.2中的目标函数求解PID控制器的参数，具体方法是：先将控制量u(k)进行变换：

u(k)＝u(k-1)+K_p(k)(e₁(k)-e₁(k-1))+K_i(k)e₁(k)

+K_d(k)(e₁(k)-2e₁(k-1)+e₁(k-2))

e₁(k)＝[e₁₁(k)e₁₂(k)…e_1q(k)]^T

其中，K_p(k)、K_i(k)、K_d(k)分别是k时刻PID控制器的比例、微分、积分参数，e₁(k)是k时刻设定值与实际输出值之间的误差；

进而将控制量u(k)简化成矩阵形式：

u(k)＝u(k-1)+E(k)^Tw(k)

其中，

E_i(k)＝[e_1i(k),e_1i(k-1),e_1i(k-2)]^Τ

w_i(k)＝[w_i1(k),w_i2(k),w_i3(k)]

w(k)＝[w₁(k)w₂(k)…w_q(k)]^T

w_i1(k)＝K_pi(k)+K_ii(k)+K_di(k)

w_i2(k)＝-K_pi(k)-2K_di(k)

w_i3(k)＝K_di(k)

结合控制量u(k)的矩阵形式和步骤2.2中的目标函数，求得：

$w\left(k\right)=E\left(k\right)(-{\left({\mathrm{\γ}}^T{Q}_e\mathrm{\γ}E{\left(k\right)}^{T}E\left(k\right)\right)}^{-1}{\mathrm{\γ}}^T{Q}_e(A_e^P{z}_e\left(k\right)+\mathrm{\ζ}\mathrm{\Δ}R\left)\right)$

进一步得到：

K_pi(k)＝-w_i2(k)-2K_di(k)

K_ii(k)＝w_i1(k)-K_pi(k)-K_di(k)

K_di(k)＝w_i3(k)

2.4得到PID控制器的参数K_p(k)、K_i(k)、K_d(k)后，构成控制量：u(k)＝u(k-1)+K_p(k)(e₁(k)-e₁(k-1))+K_i(k)e₁(k)+K_d(k)(e₁(k)-2e₁(k-1)+e₁(k-2))，再将u(k)作用于被控对象；

2.5在k+l时刻，依照2.1到2.4中的步骤循环求解PID控制器新的参数K_p(k+l)、K_i(k+l)、K_d(k+l)，l＝1,2,3，…；当满足如下条件时结束循环：

$\left\{\begin{array}{c}{k}_{pi}\left(k\right)=k_{pi}(k-1)\\ k_{ii}\left(k\right)=k_{ii}(k-1)......\left|e_{1i}\left(k\right)\right|\≤\mathrm{\ϵ}\\ k_{di}\left(k\right)=k_{di}(k-1)\end{array}\right.$

$\left\{\begin{array}{c}{k}_{pi}\left(k\right)=-w_{i2}\left(k\right)-2k_{di}\left(k\right)\\ k_{ii}=w_{i1}\left(k\right)-k_{pi}\left(k\right)-k_{di}\left(k\right)......\left|e_{1i}\left(k\right)\right|>\mathrm{\ϵ}\\ k_{di}\left(k\right)=w_{i3}\left(k\right)\end{array}\right..$