[发明专利]炼油加热炉炉膛压力的多变量预测控制PID控制方法

权利要求书

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(k)=y1(k)y2(k)...yq(k),U(k)=u1(k)u2(k)...up(k)]]>

Si=b1i1b1i2...b1ipb2i1b2i2...b2ip............bqi1bqi2...bqip,i=0,1,...,n]]>

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

Yj(k)=Hjθj,θjT=[L1...Lnbj01...bj0pbj11...bj1p...bjn1...bjnp]]]>

Hj=-yj(n)...-yj(1)UT(n+1)...UT(1)-yj(n+1)...-yj(2)UT(n+2)...UT(2)..................-yj(n+N-1)...-yj(N)UT(n+N)...UT(N)]]>

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

θ^j=(HjTHj)-1HjTYj(k)]]>

按照上式，令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)

其中，

Am=-L1-L2...-Ln-1-LnS2...Sn-1SnIq0...000...000Iq...000...00...........................00...Iq00...0000...000...0000...00Ip...00...........................00...000...Ip0]]>

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

C_m＝[I_q00…0000]

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

1.4选取新的状态变量

z(k)=Δxm(k)y(k)]]>

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

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

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

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

其中，A=Am0CmAmIq,B=BmCmBm,]]>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)

此时选取一个新的状态变量ze(k)=Δxm(k)e(k)]]>

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

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

其中，

Ae=Am0CmAmIq;Be=BmCmBm;Ce=0-Iq]]>

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

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

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

其中，

ζ=AeP-1CeAeP-2Ce...Ce,γ=AeP-1Be]]>

△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(k)=E(k)(-(γTQeγE(k)TE(k))-1γTQe(AePze(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，…；当满足如下条件时结束循环：

kpi(k)=kpi(k-1)kii(k)=kii(k-1)......|e1i(k)|≤ϵkdi(k)=kdi(k-1)]]>

kpi(k)=-wi2(k)-2kdi(k)kii=wi1(k)-kpi(k)-kdi(k)......|e1i(k)|>ϵkdi(k)=wi3(k).]]>