CN105159095A - 多变量过程的蒸馏塔模型预测控制优化pid控制方法 - Google Patents

Abstract

本发明公开了一种多变量过程的蒸馏塔模型预测控制优化PID控制方法，本发明结合扩展非最小状态空间模型预测控制和传统PID控制算法。首先建立多变量过程的蒸馏塔状态空间模型，挖掘出基本的过程特性；然后结合蒸馏塔内多变量过程的状态过程和输出误差，建立扩展非最小状态空间模型。在模型的基础上，基于优化思想，利用模型预测控制方法来整定PID控制器的参数，最后对蒸馏塔内多变量过程实现PID控制，有效提高了多变量系统的控制性能，明显减小系统超调量，改善了控制性能，又保证了控制结构简单；弥补了传统PID控制的不足，又保证了良好的控制性能。

Description

多变量过程的蒸馏塔模型预测控制优化PID控制方法

技术领域

本发明属于自动化技术领域，涉及一种多变量过程的蒸馏塔模型预测控制优化PID控制方法。

背景技术

在实际工业过程中，由于PID控制结构简单，操作方便，已经广泛地应用于实际过程控制系统中。但是，随着工业的发展，人们对产品质量、能源消耗和环境保护的要求越来越高，对工业过程的控制精度要求也日益严格，传统的PID控制方法难以进一步提升控制水平，控制过程变得更加复杂。针对蒸馏塔这一多变量对象，传统的PID控制的参数整定过程比较复杂，可能往往达不到所要求的控制性能指标和控制精度。模型预测控制(MPC)作为先进控制技术的一种，具有对模型精度要求低、在线计算方便、控制效果好等特点，在实际复杂工业过程中已经得到非常成功的应用。基于扩展非最小状态空间模型预测控制(ENMSSMPC)在多变量过程控制中比PID控制拥有更好的控制性能，如果能够在实际过程中将ENMSSMPC和PID控制技术相结合，系统的控制性能将得到进一步的改善，同时又能保证控制结构的形式比较简单。

发明内容

本发明的目的是针对传统PID控制在多变量过程控制的不足之处，提出了一种多变量过程的蒸馏塔模型预测控制优化PID控制方法，该方法结合扩展非最小状态空间模型预测控制和传统PID控制算法，弥补了传统PID控制的不足，又保证了良好的控制性能。

本发明方法首先建立多变量过程的蒸馏塔状态空间模型，挖掘出基本的过程特性；然后结合蒸馏塔内多变量过程的状态过程和输出误差，建立扩展非最小状态空间模型。在模型的基础上，基于优化思想，利用模型预测控制方法来整定PID控制器的参数，最后对蒸馏塔内多变量过程实现PID控制，有效提高了多变量系统的控制性能。

本发明的技术方案是通过数据采集、预测机理、滚动优化等手段，确立了在蒸馏塔内多变量过程的状态空间模型预测控制优化PID控制方法，利用该方法可明显减小系统超调量，改善了控制性能，又保证了控制结构简单。

本发明方法的步骤如下：

步骤(1).建立被控对象的扩展非最小状态空间模型，具体步骤如下：

1-1.通过采集被控对象的实时数据，用最小二乘法建立多输入多输出系统的模型

Y(k)+H₁Y(k-1)+…H_nY(k-n)＝L₁U(k)+L₂U(k-1)+……+L_nU(k-n+1)

其中，

$Y\left(k\right)=\left[\begin{array}{c}{y}_1\left(k\right)\\ y_2\left(k\right)\\ .\\ .\\ .\\ y_p\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_q\left(k\right)\end{array}\right];$

Y(k)为p维输出，U(k)为q维输入，H₁…H_n,L₁,L₂…L_n是系统需要辨识的系数。

通过最小二乘法辨识系统模型结果如下：

其中，

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

1-2.将步骤1-1中辨识得到的模型进一步处理为如下形式：

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

其中，y(k)和u(k)分别是k时刻输出和输入，△为后移算子。

1-3.选取非最小状态空间变量△x(k)，形式如下：

△x(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]

其中，△x(k)的维数是m＝p×(n-1)+q×n。

将步骤1-2中的模型经过转换后，可得状态空间模型：

△x(k+1)＝A△x(k)+B△u(k)

△y(k+1)＝C△x(k+1)

其中，

$A=\left[\begin{array}{ccccccccc}-H_1& -H_1& \mathrm{...}& -H_1& -H_1& -H_1& \mathrm{...}& -H_1& -H_1\\ 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_q& \mathrm{...}& 0& 0\\ .& .& & .& .& .& & .& .\\ .& .& \mathrm{...}& .& .& .& \mathrm{...}& .& .\\ .& .& & .& .& .& & .& .\\ 0& 0& \mathrm{...}& 0& 0& 0& \mathrm{...}& I_q& 0\end{array}\right]$

B＝[L₁00…0I_p00]

C＝[I_q00…0000]

△x(k+1)、△y(k+1)分别是第k+1时刻的状态和输出，△u(k)是第k时刻的输入变量增量值，A、B、C分别对应的是状态矩阵、输入矩阵和输出矩阵。I_p为p维的单位矩阵，I_q为q维的单位矩阵。

1-4定义预期的输出r(k)，则跟踪误差表示为：

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

结合步骤1-3中的状态空间模型和定义的跟踪误差得到：

e(k+1)＝e(k)+CA△x(k)+CB△u(k)-△r(k+1)

其中△u(k)，△r(k+1)分别是经过后移算子后的输入和预期的输出。

为了获得扩展非最小状态空间模型，构造一个新的状态变量如下：

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

进一步扩展状态空间模型，形式如下：

z(k+1)＝A_mz(k)+B_m△u(k)+C_m△r(k+1)

其中，

$A_m=\left[\begin{array}{cc}A& 0\\ CA& I_q\end{array}\right],B_m=\left[\begin{array}{c}B\\ CB\end{array}\right],C_m=\left[\begin{array}{c}0\\ -I_q\end{array}\right]$

0是一个m×q维的零矩阵，I_q是一个q维单位矩阵。

步骤(2).设计多变量过程的PID控制器，具体步骤是：

2-1.为了预测未来时刻的状态，定义如下形式的矩阵：

$z=\left[\begin{array}{c}z(k+1)\\ z(k+2)\\ .\\ .\\ .\\ z(k+P)\end{array}\right],\mathrm{\Δ}R=\left[\begin{array}{c}\mathrm{\Δ}r(k+1)\\ \mathrm{\Δ}r(k+2)\\ .\\ .\\ .\\ \mathrm{\Δ}r(k+3)\end{array}\right]$

r(k+i)＝αⁱy(k)+(1-αⁱ)c(k)

其中，P是预测时域，α是参考轨迹柔化系数矩阵，c(k)是在k时刻的设定值。

将未来时刻的状态写成如下形式：

Z＝Sz(k)+F△u(k)+θ△R

$S=\left[\begin{array}{c}{A}_m\\ {A_m}^2\\ .\\ .\\ .\\ {A_m}^P\end{array}\right],$ $F=\left[\begin{array}{c}{B}_m\\ A_m{B}_m\\ .\\ .\\ .\\ {A_m}^{P-1}{B}_m\end{array}\right];$

2-1.选择性能指标函数J，形式如下：

J(k)＝Z^TQZ+△u(k)^Tγ△u(k)

Q＝blockdiag{Q₁Q₂…Q_P}

其中，Q和γ分别是多变量过程的状态和输入的权重矩阵。

2-3.结合扩展非最小状态空间预测控制(ENMSSMPC)优化思想，设计PID控制器，将控制量写成增量形式，用如下式子表述：

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)＝c(k)-y(k)

e₁(k)＝[e₁₁(k)e₁₁(k)…e_1q(k)]^Τ

其中，k_p(k),k_i(k),k_d(k)分别代表的是在k时刻比例系数矩阵、积分系数矩阵和微分系数矩阵，e₁(k)是在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}\left(k\right)& 0\\ 0& 0& \mathrm{...}& 0& E_p\left(k\right)\\ .& .& .& .& \\ .& .& .& .& \\ .& .& .& .& 0\end{array}\right]}_{3q\×p}$

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

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

w_i(k)＝[w_i1(k)w_i2(k)w_i3(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)

求解上述式子，得到PID控制器的参数向量：

w(k)＝E(k)(-((F^TQF+γ)E(k)^TE(k))^-1F^TQ(Sz(k)+θ△R))

则得到的多变量控制系统的PID控制参数如下：

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)＝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))。

2-5.在k+l时刻，重复步骤2-1至2-4，继续求解新的控制量u(k+1)，并依次循环求解k_p(k+l),k_i(k+l),k_d(k+l)，l＝1,2,3，…。

本发明有益效果如下：

本发明方法首先基于实时多输入多输出数据建立对象的状态空间模型，然后结合多输入多输出对象的状态过程和输出误差建立扩展的非最小状态空间模型。在此模型的基础上，依据ENMSSMPC的方法来优化PID控制器的参数，最后对多变量过程对象实现PID控制，该方法明显地减小了系统超调量，获得了良好的控制效果。

具体实施方式

以蒸馏塔这一多变量过程为例：

这里以蒸馏塔内多个变量控制为例加以描述，蒸馏塔内是一个复杂的情形，同时可能受到多种不同的控制量的影响。

步骤(1).建立被控对象的扩展非最小状态空间模型，具体步骤如下：

1-1.通过采集被控对象的实时数据，用最小二乘法建立多输入多输出系统的模型

Y(k)+H₁Y(k-1)+…H_nY(k-n)＝L₁U(k)+L₂U(k-1)+……+L_nU(k-n+1)

其中，

$Y\left(k\right)=\left[\begin{array}{c}{y}_1\left(k\right)\\ y_2\left(k\right)\\ .\\ .\\ .\\ y_p\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_q\left(k\right)\end{array}\right];$

Y(k)为p维输出，U(k)为q维输入，H₁…H_n,L₁,L₂…L_n是系统需要辨识的系数。

通过最小二乘法辨识系统模型结果如下：

其中，

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

1-2.将步骤1-1中辨识得到的模型进一步处理为如下形式：

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

其中，y(k)和u(k)分别是k时刻输出和输入，△为后移算子。

1-3.选取非最小状态空间变量△x(k)，形式如下：

△x(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]

其中，△x(k)的维数是m＝p×(n-1)+q×n。

将步骤1-2中的模型经过转换后，可得状态空间模型：

△x(k+1)＝A△x(k)+B△u(k)

△y(k+1)＝C△x(k+1)

其中，

$A=\left[\begin{array}{ccccccccc}-H_1& -H_1& \mathrm{...}& -H_1& -H_1& -H_1& \mathrm{...}& -H_1& -H_1\\ 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_q& \mathrm{...}& 0& 0\\ .& .& & .& .& .& & .& .\\ .& .& \mathrm{...}& .& .& .& \mathrm{...}& .& .\\ .& .& & .& .& .& & .& .\\ 0& 0& \mathrm{...}& 0& 0& 0& \mathrm{...}& I_q& 0\end{array}\right]$

B＝[L₁00…0I_p00]

C＝[I_q00…0000]

△x(k+1)、△y(k+1)分别是第k+1时刻的状态和输出，△u(k)是第k时刻的输入变量增量值，A、B、C分别对应的是状态矩阵、输入矩阵和输出矩阵。I_p为p维的单位矩阵，I_q为q维的单位矩阵。

1-4定义预期的输出r(k)，则跟踪误差表示为：

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

结合步骤1-3中的状态空间模型和定义的跟踪误差得到：

e(k+1)＝e(k)+CA△x(k)+CB△u(k)-△r(k+1)

其中△u(k)，△r(k+1)分别是经过后移算子后的输入和预期的输出。

为了获得扩展非最小状态空间模型，构造一个新的状态变量如下：

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

进一步扩展状态空间模型，形式如下：

z(k+1)＝A_mz(k)+B_m△u(k)+C_m△r(k+1)

其中，

$A_m=\left[\begin{array}{cc}A& 0\\ CA& I_q\end{array}\right],B_m=\left[\begin{array}{c}B\\ CB\end{array}\right],C_m=\left[\begin{array}{c}0\\ -I_q\end{array}\right]$

0是一个m×q维的零矩阵，I_q是一个q维单位矩阵。

步骤(2).设计多变量过程的PID控制器，具体步骤是：

2-1.为了预测未来时刻的状态，定义如下形式的矩阵：

$z=\left[\begin{array}{c}z(k+1)\\ z(k+2)\\ .\\ .\\ .\\ z(k+P)\end{array}\right],\mathrm{\Δ}R=\left[\begin{array}{c}\mathrm{\Δ}r(k+1)\\ \mathrm{\Δ}r(k+2)\\ .\\ .\\ .\\ \mathrm{\Δ}r(k+3)\end{array}\right]$

r(k+i)＝αⁱy(k)+(1-αⁱ)c(k)

其中，P是预测时域，α是参考轨迹柔化系数矩阵，c(k)是在k时刻的设定值。

将未来时刻的状态写成如下形式：

Z＝Sz(k)+F△u(k)+θ△R

$S=\left[\begin{array}{c}{A}_m\\ {A_m}^2\\ .\\ .\\ .\\ {A_m}^P\end{array}\right],$ $F=\left[\begin{array}{c}{B}_m\\ A_m{B}_m\\ .\\ .\\ .\\ {A_m}^{P-1}{B}_m\end{array}\right];$

2-1.选择性能指标函数J，形式如下：

J(k)＝Z^TQZ+△u(k)^Tγ△u(k)

Q＝blockdiag{Q₁Q₂…Q_P}

其中，Q和γ分别是多变量过程的状态和输入的权重矩阵。

2-3.结合扩展非最小状态空间预测控制(ENMSSMPC)优化思想，设计PID控制器，将控制量写成增量形式，用如下式子表述：

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)＝c(k)-y(k)

e₁(k)＝[e₁₁(k)e₁₁(k)…e_1q(k)]^Τ

其中，k_p(k),k_i(k),k_d(k)分别代表的是在k时刻比例系数矩阵、积分系数矩阵和微分系数矩阵，e₁(k)是在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}\left(k\right)& 0\\ 0& 0& \mathrm{...}& 0& E_p\left(k\right)\\ .& .& .& .& \\ .& .& .& .& \\ .& .& .& .& 0\end{array}\right]}_{3q\×p}$

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

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

w_i(k)＝[w_i1(k)w_i2(k)w_i3(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)

求解上述式子，得到PID控制器的参数向量：

w(k)＝E(k)(-((F^TQF+γ)E(k)^TE(k))^-1F^TQ(Sz(k)+θ△R))

则得到的多变量控制系统的PID控制参数如下：

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)＝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))。

2-5.在k+l时刻，重复步骤2-1至2-4，继续求解新的控制量u(k+1)，并依次循环求解k_p(k+l),k_i(k+l),k_d(k+l)，l＝1,2,3，…。

Claims (2)

1.多变量过程的蒸馏塔模型预测控制优化PID控制方法，其特征在于包括如下步骤：

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

步骤(2).设计多变量过程的PID控制器；

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

1-1.通过采集被控对象的实时数据，用最小二乘法建立多输入多输出系统的模型

Y(k)+H₁Y(k-1)+…H_nY(k-n)＝L₁U(k)+L₂U(k-1)+……+L_nU(k-n+1)

其中，

$Y\left(k\right)=\left[\begin{array}{c}{y}_1\left(k\right)\\ y_2\left(k\right)\\ \·\\ \·\\ \·\\ y_p\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_q\left(k\right)\end{array}\right];$

Y(k)为p维输出，U(k)为q维输入，H₁…H_n,L₁,L₂…L_n是系统需要辨识的系数；

通过最小二乘法辨识系统模型结果如下：

其中，

$Y_j=\left[\begin{array}{c}{y}_j(n+1)\\ y_j(n+2)\\ \begin{array}{c}\·\\ \·\\ \·\end{array}\\ y_j(n+N)\end{array}\right]H_j=\left[\begin{array}{cccccc}-y_j\left(n\right)& ...& -y_j\left(1\right)& U^T(n+1)& ...& U^T\left(1\right)\\ -y_j(n+1)& ...& -y_j\left(2\right)& U^T(n+2)& ...& U^T\left(2\right)\\ \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}\\ -y_j(n+N-1)& ...& -y_j\left(N\right)& U^T(n+N)& ...& U^T\left(N\right)\end{array}\right]$

1-2.将步骤1-1中辨识得到的模型进一步处理为如下形式：

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

其中，y(k)和u(k)分别是k时刻输出和输入，△为后移算子；

1-3.选取非最小状态空间变量△x(k)，形式如下：

△x(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]

其中，△x(k)的维数是m＝p×(n-1)+q×n；

将步骤1-2中的模型经过转换后，可得状态空间模型：

△x(k+1)＝A△x(k)+B△u(k)

△y(k+1)＝C△x(k+1)

其中，

$A=\left[\begin{array}{ccccccccc}-H_1& -H_1& ...& -H_1& -H_1& -H_1& ...& -H_1& -H_1\\ 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}\\ 0& 0& ...& I_q& 0& 0& ...& 0& 0\\ 0& 0& ...& 0& 0& 0& ...& 0& 0\\ 0& 0& ...& 0& 0& I_q& ...& 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}\\ 0& 0& ...& 0& 0& 0& ...& I_q& 0\end{array}\right]$

B＝[L₁00…0I_p00]

C＝[I_q00…0000]

△x(k+1)、△y(k+1)分别是第k+1时刻的状态和输出，△u(k)是第k时刻的输入变量增量值，A、B、C分别对应的是状态矩阵、输入矩阵和输出矩阵；I_p为p维的单位矩阵，I_q为q维的单位矩阵；

1-4定义预期的输出r(k)，则跟踪误差表示为：

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

结合步骤1-3中的状态空间模型和定义的跟踪误差得到：

e(k+1)＝e(k)+CA△x(k)+CB△u(k)-△r(k+1)

其中△u(k)，△r(k+1)分别是经过后移算子后的输入和预期的输出；

为了获得扩展非最小状态空间模型，构造一个新的状态变量如下：

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

进一步扩展状态空间模型，形式如下：

z(k+1)＝A_mz(k)+B_m△u(k)+C_m△r(k+1)

其中，

$A_m=\left[\begin{array}{cc}A& 0\\ CA& I_q\end{array}\right],B_m=\left[\begin{array}{c}B\\ CB\end{array}\right],C_m=\left[\begin{array}{c}0\\ -I_q\end{array}\right]$

0是一个m×q维的零矩阵，I_q是一个q维单位矩阵。

2.根据权利要求1所述的多变量过程的蒸馏塔模型预测控制优化PID控制方法，其特征在于步骤(2)所述的设计多变量过程的PID控制器，具体步骤如下：

2-1.为了预测未来时刻的状态，定义如下形式的矩阵：

$z=\left[\begin{array}{c}z(k+1)\\ z(k+2)\\ \·\\ \·\\ \·\\ z(k+P)\end{array}\right],\mathrm{\Δ}R=\left[\begin{array}{c}\mathrm{\Δ}r(k+1)\\ \mathrm{\Δ}r(k+2)\\ \·\\ \·\\ \·\\ \mathrm{\Δ}r(k+3)\end{array}\right]$

r(k+i)＝αⁱy(k)+(1-αⁱ)c(k)

其中，P是预测时域，α是参考轨迹柔化系数矩阵，c(k)是在k时刻的设定值；

将未来时刻的状态写成如下形式：

Z＝Sz(k)+F△u(k)+θ△R

2-1.选择性能指标函数J，形式如下：

J(k)＝Z^TQZ+△u(k)^Tγ△u(k)

Q＝blockdiag{Q₁Q₂…Q_P}

其中，Q和γ分别是多变量过程的状态和输入的权重矩阵；

2-3.结合扩展非最小状态空间预测控制优化思想，设计PID控制器，将控制量写成增量形式，用如下式子表述：

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)＝c(k)-y(k)

e₁(k)＝[e₁₁(k)e₁₁(k)…e_1q(k)]^Τ

其中，k_p(k),k_i(k),k_d(k)分别代表的是在k时刻比例系数矩阵、积分系数矩阵和微分系数矩阵，e₁(k)是在k时刻设定值矩阵和实际输出矩阵的差值；

进一步将控制量简化为如下形式：

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

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

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

w_i(k)＝[w_i1(k)w_i2(k)w_i3(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)

求解上述式子，得到PID控制器的参数向量：

w(k)＝E(k)(-((F^TQF+γ)E(k)^TE(k))^-1F^TQ(Sz(k)+θ△R))

则得到的多变量控制系统的PID控制参数如下：

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)＝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))；

2-5.在k+l时刻，重复步骤2-1至2-4，继续求解新的控制量u(k+1)，并依次循环求解k_p(k+l),k_i(k+l),k_d(k+l)，l＝1,2,3，…。