CN103048927B

CN103048927B - 用于精馏系统的模型预测控制方法

Info

Publication number: CN103048927B
Application number: CN201210586805.XA
Authority: CN
Inventors: 谢磊; 冯皓; 苏宏业
Original assignee: Zhejiang University ZJU
Current assignee: Zhejiang University ZJU
Priority date: 2012-12-28
Filing date: 2012-12-28
Publication date: 2015-03-25
Anticipated expiration: 2032-12-28
Also published as: CN103048927A

Abstract

本发明公开了一种用于精馏系统的模型预测控制方法，步骤如下：1)建立目标函数；2)构造状态空间模型；3)利用子空间方法，构造Hankel矩阵以获得输出变量的子空间函数；4)根据子空间函数求得目标函数最小化状态下的最优控制率；5)将最优控制率作用于精馏系统中，获得各回路输入变量和输出变量的方差；6)通过粒子群优化方法寻找目标函数中权衡矩阵的一次最优值；7)以一次最优值为初值，利用梯度下降法获得二次最优值，将二次最优值代入目标函数以获得最优化的模型预测控制方法。本发明仅通过过程的输入输出数据就可以获得过程特性，对权重矩阵进行精确优化的方法，从而将产品质量和控制代价最优化，提升整体的经济效益。

Description

用于精馏系统的模型预测控制方法

技术领域

本发明涉及一种模型预测控制方法，尤其是一种用于精馏系统的模型预测控制方法。

背景技术

在现代的工业生产中，为了提高经济效益，该生产系统中的控制系统不仅要求产品的质量波动最小，还要求对应的控制成本最低，而模型预测控制器则常被用于平衡这两者间的关系。

模型预测控制器自20世纪70年代诞生以来，在理论研究和工程实践方面都得到了蓬勃发展，特别是在流程工业中得到了广泛应用。特别是在复杂的工业过程中有着越来越广泛的应用，其性能监控技术的研究对提高过程控制质量、提高企业生产效率、促进先进控制的推广具有重要意义。模型预测控制器设计的最终目标是提高经济效益，虽然权重矩阵确定了模型预测控制器的目标函数，但这一目标函数并不直接对应经济效益，提高经济效益可以归结到将过程输出和控制器输出的波动约束到一个范围内，工业现场对模型预测控制器的设计首先应该通过过程输出和控制器输出的方差来确定权重矩阵，再完成控制器的设计。

在以往的模型预测控制器设计中，权重矩阵的选取没有明确的标准，对权重矩阵的优化没有一种对应到过程输出以及控制器输出的方法。Huang的书(Huang B,Shah S L.Performance assessment of control loops:theory and applications[M].Springer,1999.)中描述的权衡曲线描述了这三者间的关系，虽然所用到的假设较为简单，模型预测控制器的权重仅为标量，但为解决权重矩阵的优化问题提供了思路。传统上有两种方法解决问题，一种是通过输入输出模型或者状态空间模型直接计算，另一种是通过广义预测控制(GPC)问题来求解。这两种方法对于模型的限制和依赖程度较高，为了获得精确的模型，往往需要对过程施加一些会引起过程波动较大的激励信号，这在工业现场往往是不被允许的。同时在计算过程中，需要求解Diophantine方程或代数Riccati方程，造成计算量较大，在实际应用中显得较为复杂。权衡曲线提供了已知权重矩阵求解过程输出方差和控制器输出方差的思路，但模型预测控制器设计需要的是已知输入输出方差来求解权衡矩阵，此时权衡曲线没有一个显式的函数关系就会造成优化层级上难度大，精度不高，容易出现局部最优等问题。

因此，针对实际情况，目前需要一种能够不用辨识出过程模型，无需过程机理知识，仅通过过程的输入输出数据就可以获得过程特性，并进一步获得模型预测控制器权重矩阵与过程输入输出方差的显式关系，以此对权重矩阵进行精确优化的方法，能够将产品质量和控制代价最优化，提升整体的经济效益。

发明内容

为了克服现有技术的不足，本发明提供了一种无需辨识出过程模型，仅通过过程的输入输出数据就可以获得过程特性，能够将产品质量和控制代价最优化，提升整体的经济效益的用于精馏系统的模型预测控制方法。

本发明的有益效果是：一种用于精馏系统的模型预测控制方法，步骤如下

1）建立模型预测控制器的目标函数；

2）构造精馏系统的状态空间模型；

3）利用子空间方法，根据精馏系统的输入变量和输出变量构造Hankel矩阵，并将Hankel矩阵代入步骤2）中的状态空间模型，并转化为输出变量的子空间函数；

4）根据步骤3）中输出变量的子空间函数获得其预测值，并将该预测值代入步骤1）中的目标函数内，通过对输入变量求导，求得目标函数最小化状态下的最优控制率；

5）将最优控制率作用于精馏系统中，闭环后计算输入变量和输出变量的方差，让扰动从某一时刻进入，即之前时刻的输入变量和扰动均为0，此时获得各回路输入变量和输出变量的方差；

6）令输入变量和输出变量的波动目标为输入变量和输出变量的方差，通过粒子群优化方法寻找目标函数中权衡矩阵的一次最优值；

7）以权衡矩阵的一次最优值为初值，利用梯度下降法对一次最优值进行进一步的优化，获得二次最优值，将二次最优值代入目标函数以获得最优化的模型预测控制方法。

进一步优选是，步骤1）中所建立的模型预测控制器的目标函数为

J = Σ_{k = 1}^{N} [{({\hat{y}}_{k} - r_{k})}^{T} Q ({\hat{y}}_{k} - r_{k}) + u_{k}^{T} {Ru}_{k}] - - - (1)

式中，

N为数据长度；

k∈[1,2,…,N]；

r_k为人为设定的目标值，r_k为常数；

u为精馏系统的输入变量；

y为精馏系统的输出变量；

为精馏系统的输出变量的预测值；

Q为输出的权衡矩阵；

R为输入的权衡矩阵。

进一步优选是，步骤2）中精馏系统的状态空间模型为

\begin{matrix} x_{l} = {Ax}_{l - 1} + {Bu}_{l - 1} + {Ee}_{l - 1} \\ y_{l} = {Cy}_{l - 1} + {Du}_{l - 1} + e_{l - 1} \end{matrix} - - - (2)

式中，A、B、C、D、E为状态空间模型矩阵；

x为状态变量；

e为噪声；

l为精馏系统中某时刻的时间；

进一步优选是，步骤3）利用子空间方法，根据精馏系统的输入变量和输出变量构造Hankel矩阵，

输入变量的Hankel矩阵：

u_{p} \overset{Δ}{=} (\begin{matrix} u_{0} & u_{1} & u_{2} & \cdot \cdot \cdot & u_{j - 1} \\ u_{1} & u_{2} & u_{3} & \cdot \cdot \cdot & u_{j} \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \cdot \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ u_{i - 1} & u_{i} & u_{i + 1} & \cdot \cdot \cdot & u_{i + j - 2} \end{matrix}),

u_{f} \overset{Δ}{=} (\begin{matrix} u_{i} & u_{i + 1} & u_{i + 2} & \cdot \cdot \cdot & u_{i + j - 1} \\ u_{i + 1} & u_{i + 2} & u_{i + 3} & \cdot \cdot \cdot & u_{i + j} \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \cdot \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ u_{2 i - 1} & u_{2 i} & u_{2 i + 1} & \cdot \cdot \cdot & u_{2 i + j - 2} \end{matrix})

（3-1）

输出变量的Hankel矩阵：

y_{p} \overset{Δ}{=} (\begin{matrix} y_{0} & y_{1} & y_{2} & \cdot \cdot \cdot & y_{j - 1} \\ y_{1} & y_{2} & y_{3} & \cdot \cdot \cdot & y_{j} \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \cdot \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ y_{i - 1} & y_{i} & y_{i + 1} & \cdot \cdot \cdot & y_{i + j - 2} \end{matrix}),

y_{f} \overset{Δ}{=} (\begin{matrix} y_{i} & y_{i + 1} & y_{i + 2} & \cdot \cdot \cdot & y_{i + j - 1} \\ y_{i + 1} & y_{i + 2} & y_{i + 3} & \cdot \cdot \cdot & y_{i + j} \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \cdot \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ y_{2 i - 1} & y_{2 i} & y_{2 i + 1} & \cdot \cdot \cdot & y_{2 i + j - 2} \end{matrix})

（3-2）

噪声的Hankel矩阵：

e_{p} \overset{Δ}{=} (\begin{matrix} e_{0} & e_{1} & e_{2} & \cdot \cdot \cdot & e_{j - 1} \\ e_{1} & e_{2} & e_{3} & \cdot \cdot \cdot & e_{j} \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \cdot \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ e_{i - 1} & e_{i} & e_{i + 1} & \cdot \cdot \cdot & e_{i + j - 2} \end{matrix}),

e_{f} \overset{Δ}{=} (\begin{matrix} e_{i} & e_{i + 1} & e_{i + 2} & \cdot \cdot \cdot & e_{i + j - 1} \\ e_{i + 1} & e_{i + 2} & e_{i + 3} & \cdot \cdot \cdot & e_{i + j} \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \cdot \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ e_{2 i - 1} & e_{2 i} & e_{2 i + 1} & \cdot \cdot \cdot & e_{2 i + j - 2} \end{matrix})

（3-3）

式中，i为Hankel矩阵的行数，同时为精馏系统工作时的时刻；

j为Hankel矩阵的列数；

定义精馏系统工作时某一时刻为t，则也为Hankel矩阵的行数；

P代表从0到t-1时刻开始的信息；

f代表从t到t+i-1时刻开始的信息；

将输入变量的Hankel矩阵（3-1）、输出变量的Hankel矩阵（3-2）和噪声的Hankel矩阵（3-3）代入状态空间函数(2)，并转换为子空间函数为

y_f＝L_ww_p+L_uu_f+L_ee_f (4)

式中，

w_{p} = (\begin{matrix} y_{p} \\ u_{p} \end{matrix})

L_w为子空间的系数矩阵；

L_u子空间的系数矩阵；

L_e子空间的系数矩阵。

进一步优选是，步骤（4）中，获得子空间函数式（4）的预测值为（4-1），将式（4-1）代入目标函数（1）中，得到

J = {(L_{w} w_{p} + L_{u} u_{f} - r_{k})}^{T} Q (L_{w} w_{p} + L_{u} u_{f} - r_{k}) + u_{f}^{T} {Ru}_{f} - - - (5)

同时，对式（5）求u_f的偏导并置为0，在J最小化的前提下获得最优控制率

{u_{f}}^{opt} = {(R + L_{u}^{T} Q L_{u})}^{- 1} L_{u}^{T} Q (r_{k} - L_{w} w_{p}) - - - (6)

令r_k＝0，且则式(6)可简化为

{u_{f}}^{opt} = - K (L_{g} {\hat{u}}_{p} + L_{h} {\hat{e}}_{p}) - - - (7)

将式(7)代入式(4)中可以获得模型预测控制器作用下的输出

y_{f} = L_{g} \hat{u} + L_{h} \hat{e} + L_{u} u_{f} + L_{e} e_{f}

(8)

= (I - L_{u} K) L_{g} {\hat{u}}_{p} + (I - L_{u} K) L_{h} {\hat{e}}_{p} + {\hat{L}}_{e} {\hat{e}}_{f}

式中，I为单位矩阵；

为u的预测值；

为e的预测值；

L_{g} = [\begin{matrix} g_{1} & g_{2} & \cdot \cdot \cdot & g_{N - 1} & g_{N} \\ g_{2} & g_{3} & \cdot \cdot \cdot & g_{N - 1} & 0 \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ g_{N} & 0 & 0 & \cdot \cdot \cdot & 0 \end{matrix}],

g为人为定义的未知参数；

L_{h} = [\begin{matrix} l_{0} & l_{1} & \cdot \cdot \cdot & l_{N - 1} & l_{N} \\ l_{1} & l_{2} & \cdot \cdot \cdot & l_{N} & 0 \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ l_{N - 1} & 0 & 0 & \cdot \cdot \cdot & 0 \end{matrix}],

l为人为定义的未知参数；

{\hat{L}}_{e} = [\begin{matrix} 0 & 0 & \cdot \cdot \cdot & 0 \\ l_{0} & 0 & \cdot \cdot \cdot & 0 \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ l_{N} & l_{N - 1} & \cdot \cdot \cdot & 0 \end{matrix}],

为L_e的预期值；

{\hat{u}}_{p} = [\begin{matrix} u_{t} \\ u_{t - 1} \\ \cdot \cdot \cdot \\ u_{t - N + 1} \end{matrix}],

为u_p的预期值；

{\hat{e}}_{p} = [\begin{matrix} e_{t + 1} \\ e_{t} \\ \cdot \cdot \cdot \\ e_{t - N + 1} \end{matrix}],

为e_p的预期值；

{\hat{e}}_{f} = [\begin{matrix} e_{t + 2} \\ e_{t + 3} \\ \cdot \cdot \cdot \\ e_{t + N + 1} \end{matrix}],

为e_f的预期值。

进一步优选是，步骤（5）中，计算精馏系统的输入变量u和精馏系统的输出变量y的方差，让扰动从t+1时刻进入，即之前时刻的控制输入和扰动均为0，可得

u_{f} = - K l_{e} e_{t + 1} = (\begin{matrix} ψ_{0} \\ ψ_{1} \\ \cdot \\ \cdot \\ \cdot \\ ψ_{N - 1} \end{matrix}) e_{t + 1},

y_{f} = (I - L_{u} K) l_{e} e_{t + 1} = (\begin{matrix} γ_{0} \\ γ_{1} \\ \cdot \\ \cdot \\ \cdot \\ γ_{N - 1} \end{matrix}) e_{t + 1}

式中，l_e为L_e的第一列，即噪声模型的脉冲响应；

(\begin{matrix} ψ_{0} \\ ψ_{1} \\ \cdot \\ \cdot \\ \cdot \\ ψ_{N - 1} \end{matrix}) = - K l_{e},

(\begin{matrix} γ_{0} \\ γ_{1} \\ \cdot \\ \cdot \\ \cdot \\ γ_{N - 1} \end{matrix}) = (I - L_{u} K) .

此时可以获得各回路输入变量u和输出变量y的方差分别为

Var [u_{t}] = Σ_{i = 0}^{N - 1} ψ_{i} Var [e_{t}] {ψ_{i}}^{T},

Var [y_{t}] = Σ_{i = 0}^{N - 1} γ_{i} Var [e_{t}] {γ_{i}}^{T} .

进一步优选是，步骤（6）中，给定系统输入变量和输出变量的波动目标分别为和通过粒子群优化方法寻找Q和R的最优值Q1和R1。

进一步优选是，步骤（7）中，以Q1和R1为初值，利用梯度下降法对Q和R进行进一步的优化，获得最优值Q2和R2,Q2和R2即为要求取的权重矩阵。

本发明中，一种能够不用辨识出过程模型，无需过程机理知识，仅通过过程的输入输出数据就可以获得过程特性，并进一步获得模型预测控制器权重矩阵与过程输入输出方差的显式关系，以此对权重矩阵进行精确优化的方法，从而将产品质量和控制代价最优化，提升整体的经济效益。

附图说明

图1为本发明实施例苯-甲苯常压筛板精馏塔的示意图。

具体实施方式

下面结合附图对本发明作进一步描述：

1）建立模型预测控制器的目标函数，该目标函数为

J = Σ_{k = 1}^{N} [{({\hat{y}}_{k} - r_{k})}^{T} Q ({\hat{y}}_{k} - r_{k}) + u_{k}^{T} {Ru}_{k}] - - - (1)

N为数据长度；

k∈[1,2,…,N]；

r_k为人为设定的目标值，r_k为常数；

u为精馏系统的输入变量；

y为精馏系统的输出变量；

为精馏系统的输出变量的预测值；

Q为输出的权衡矩阵；

R为输入的权衡矩阵。

2）构造精馏系统的状态空间模型，该状态空间模型为

x_l＝Ax_l-1+Bu_l-1+Ee_l-1

(2)

y_l＝Cy_l-1+Du_l-1+e_l-1

式中，A、B、C、D、E为状态空间模型矩阵；

x为状态变量；

e为噪声；

l为精馏系统中某时刻的时间。

3）利用子空间方法，根据精馏系统的输入变量、输出变量和噪声构造Hankel矩阵，

输入变量的Hankel矩阵：

u_{p} \overset{Δ}{=} (\begin{matrix} u_{0} & u_{1} & u_{2} & \cdot \cdot \cdot & u_{j - 1} \\ u_{1} & u_{2} & u_{3} & \cdot \cdot \cdot & u_{j} \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \cdot \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ u_{i - 1} & u_{i} & u_{i + 1} & \cdot \cdot \cdot & u_{i + j - 2} \end{matrix}),

u_{f} \overset{Δ}{=} (\begin{matrix} u_{i} & u_{i + 1} & u_{i + 2} & \cdot \cdot \cdot & u_{i + j - 1} \\ u_{i + 1} & u_{i + 2} & u_{i + 3} & \cdot \cdot \cdot & u_{i + j} \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \cdot \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ u_{2 i - 1} & u_{2 i} & u_{2 i + 1} & \cdot \cdot \cdot & u_{2 i + j - 2} \end{matrix})

（3-1）；

输出变量的Hankel矩阵：

y_{p} \overset{Δ}{=} (\begin{matrix} y_{0} & y_{1} & y_{2} & \cdot \cdot \cdot & y_{j - 1} \\ y_{1} & y_{2} & y_{3} & \cdot \cdot \cdot & y_{j} \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \cdot \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ y_{i - 1} & y_{i} & y_{i + 1} & \cdot \cdot \cdot & y_{i + j - 2} \end{matrix}),

y_{f} \overset{Δ}{=} (\begin{matrix} y_{i} & y_{i + 1} & y_{i + 2} & \cdot \cdot \cdot & y_{i + j - 1} \\ y_{i + 1} & y_{i + 2} & y_{i + 3} & \cdot \cdot \cdot & y_{i + j} \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \cdot \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ y_{2 i - 1} & y_{2 i} & y_{2 i + 1} & \cdot \cdot \cdot & y_{2 i + j - 2} \end{matrix})

（3-2）；

噪声的Hankel矩阵：

e_{p} \overset{Δ}{=} (\begin{matrix} e_{0} & e_{1} & e_{2} & \cdot \cdot \cdot & e_{j - 1} \\ e_{1} & e_{2} & e_{3} & \cdot \cdot \cdot & e_{j} \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \cdot \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ e_{i - 1} & e_{i} & e_{i + 1} & \cdot \cdot \cdot & e_{i + j - 2} \end{matrix}),

e_{f} \overset{Δ}{=} (\begin{matrix} e_{i} & e_{i + 1} & e_{i + 2} & \cdot \cdot \cdot & e_{i + j - 1} \\ e_{i + 1} & e_{i + 2} & e_{i + 3} & \cdot \cdot \cdot & e_{i + j} \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \cdot \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ e_{2 i - 1} & e_{2 i} & e_{2 i + 1} & \cdot \cdot \cdot & e_{2 i + j - 2} \end{matrix})

（3-3）；

式中，i为Hankel矩阵的行数，同时为精馏系统工作时的时刻；

j为Hankel矩阵的列数；

定义精馏系统工作时某一时刻为t，则t时刻也即Hankel矩阵的行数；

P代表从0到t-1时刻开始的信息；

f代表从t到t+i-1时刻开始的信息。

此处，i没有物理意义，指的是人为设定的Hankel矩阵的行数。如果i是无穷大，那么就和推导中说的一样逼近了真值，但不可能取无穷大，所以，一般需要i远大于系统的阶次。

然后，将输入变量的Hankel矩阵（3-1）、输出变量的Hankel矩阵（3-2）和噪声的Hankel矩阵（3-3）代入状态空间函数(2)，并转换为子空间函数为

y_f＝L_ww_p+L_uu_f+L_ee_f (4)

式中，

w_{p} = (\begin{matrix} y_{p} \\ u_{p} \end{matrix})

L_w为子空间的系数矩阵；

L_u子空间的系数矩阵；

L_e子空间的系数矩阵。

步骤（4）中，获得子空间函数式（4）的预测值为（4-1），将式（4-1）代入目标函数（1）中，得到

J = {(L_{w} w_{p} + L_{u} u_{f} - r_{k})}^{T} Q (L_{w} w_{p} + L_{u} u_{f} - r_{k}) + u_{f}^{T} {Ru}_{f} - - - (5)

式中，

{u_{f}}^{opt} = {(R + L_{u}^{T} Q L_{u})}^{- 1} L_{u}^{T} Q (r_{k} - L_{w} w_{p}) - - - (6)

令r_k＝0，且则式(6)可简化为

{u_{f}}^{opt} = - K (L_{g} {\hat{u}}_{p} + L_{h} {\hat{e}}_{p}) - - - (7)

将式(7)代入式(4)中可以获得模型预测控制器作用下的输出

y_{f} = L_{g} \hat{u} + L_{h} \hat{e} + L_{u} u_{f} + L_{e} e_{f}

(8)

= (I - L_{u} K) L_{g} {\hat{u}}_{p} + (I - L_{u} K) L_{h} {\hat{e}}_{p} + {\hat{L}}_{e} {\hat{e}}_{f}

式中，I为单位矩阵；

为u的预测值；

为e的预测值；

L_{g} = [\begin{matrix} g_{1} & g_{2} & \cdot \cdot \cdot & g_{N - 1} & g_{N} \\ g_{2} & g_{3} & \cdot \cdot \cdot & g_{N - 1} & 0 \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ g_{N} & 0 & 0 & \cdot \cdot \cdot & 0 \end{matrix}],

g为人为定义的未知参数；

L_{h} = [\begin{matrix} l_{0} & l_{1} & \cdot \cdot \cdot & l_{N - 1} & l_{N} \\ l_{1} & l_{2} & \cdot \cdot \cdot & l_{N} & 0 \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ l_{N - 1} & 0 & 0 & \cdot \cdot \cdot & 0 \end{matrix}],

l为人为定义的未知参数；

{\hat{L}}_{e} = [\begin{matrix} 0 & 0 & \cdot \cdot \cdot & 0 \\ l_{0} & 0 & \cdot \cdot \cdot & 0 \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ l_{N} & l_{N - 1} & \cdot \cdot \cdot & 0 \end{matrix}],

为L_e的预期值；

{\hat{u}}_{p} = [\begin{matrix} u_{t} \\ u_{t - 1} \\ \cdot \cdot \cdot \\ u_{t - N + 1} \end{matrix}],

为u_p的预期值；

{\hat{e}}_{p} = [\begin{matrix} e_{t + 1} \\ e_{t} \\ \cdot \cdot \cdot \\ e_{t - N + 1} \end{matrix}],

为e_p的预期值；

{\hat{e}}_{f} = [\begin{matrix} e_{t + 2} \\ e_{t + 3} \\ \cdot \cdot \cdot \\ e_{t + N + 1} \end{matrix}],

为e_f的预期值。

5）将最优控制率作用于精馏系统中，闭环后，计算精馏系统的输入变量u和精馏系统的输出变量y的方差，让扰动从t+1时刻进入，即之前时刻的控制输入和扰动均为0，可得

u_{f} = - K l_{e} e_{t + 1} = (\begin{matrix} ψ_{0} \\ ψ_{1} \\ \cdot \\ \cdot \\ \cdot \\ ψ_{N - 1} \end{matrix}) e_{t + 1},

y_{f} = (I - L_{u} K) l_{e} e_{t + 1} = (\begin{matrix} γ_{0} \\ γ_{1} \\ \cdot \\ \cdot \\ \cdot \\ γ_{N - 1} \end{matrix}) e_{t + 1}

式中，l_e为L_e的第一列，即噪声模型的脉冲响应,

(\begin{matrix} ψ_{0} \\ ψ_{1} \\ \cdot \\ \cdot \\ \cdot \\ ψ_{N - 1} \end{matrix}) = - K l_{e},

(\begin{matrix} γ_{0} \\ γ_{1} \\ \cdot \\ \cdot \\ \cdot \\ γ_{N - 1} \end{matrix}) = (I - L_{u} K) .

其中，

(\begin{matrix} ψ_{0} \\ ψ_{1} \\ \cdot \\ \cdot \\ \cdot \\ ψ_{N - 1} \end{matrix})

是-Kl_e化简后所得到的矩阵，同理，

(\begin{matrix} γ_{0} \\ γ_{1} \\ \cdot \\ \cdot \\ \cdot \\ γ_{N - 1} \end{matrix})

是(I-L_uK)化简后所得到的矩阵。

此时可以获得各回路输入变量u和输出变量y的方差分别为

Var [u_{t}] = Σ_{i = 0}^{N - 1} ψ_{i} Var [e_{t}] {ψ_{i}}^{T},

Var [y_{t}] = Σ_{i = 0}^{N - 1} γ_{i} Var [e_{t}] {γ_{i}}^{T} .

6）给定系统输入变量和输出变量的波动目标分别为和通过粒子群优化方法寻找Q和R的最优值Q1和R1。

7）以Q1和R1为初值，利用梯度下降法对Q和R进行进一步的优化，获得最优值Q2和R2,Q2和R2即为要求取的权重矩阵。

下面针对国内某厂采用苯-甲苯常压筛板精馏塔过程中模型预测控制器的权重设定为例，对本发明具体实施方式做进一步描述。

如图1所示，苯-甲苯常压筛板精馏塔的示意图。

精馏塔1、再沸器2、冷凝器3、气液分离器4、回流泵5；

精馏塔1共41块塔板,每块塔板上液体流动的时间常数为0.06分钟，进料塔板为从下往上数第21块塔板，板上液体的相对挥发度为1.5，原料液为饱和液体，其中苯的摩尔分数ZF=50%，进料流量F=8kmol/min，回流量L=16kmol/min，再沸器2蒸汽流量V=20kmol/min，馏出液产量D=4kmol/min，釜液产量B=4kmol/min。

该过程有两个输出量和两个输入量，输出量为为釜液浓度XB和馏出液浓度XD，即y_t=[XD,XB]^T(T表示矩阵的转置)，输入量为回流量L和再沸器蒸汽流量V，即u_t=[L,V]^T。根据经济效益，需要设定两者的方差为

Var[y_t]=[1.46,1.46]^T，Var[u_t]=[0.11,0.017]^T。

根据本发明所示的步骤，首先获得一段过程状况良好时的数据，分别是系统的输出变量y_t和输入变量u_t。

1）建立模型预测控制器的目标函数为

J = Σ_{k = 1}^{N} [{({\hat{y}}_{k} - r_{k})}^{T} Q ({\hat{y}}_{k} - r_{k}) + u_{k}^{T} {Ru}_{k}];

2）构造精馏系统的状态空间模型为

x_l＝Ax_l-1+Bu_l-1+Ee_l-1

；

y_l＝Cy_l-1+Du_l-1+e_l-1

3）利用子空间方法，根据精馏系统的输入变量u_t和输出变量y_t构造Hankel矩阵。在该实施例中，选取k=20，以下选取一个数据长度N=4，k=4的情况给出更详细的计算过程，但由于实际运算中要求N和k均较大，且N>>k，以下方法只为解释计算过程，不能获得正确的结果。当采样获得

y = [\begin{matrix} 0.3273 & 0.2587 & - 7.5755 & 3.9271 \\ - 0.5883 & 1.9229 & - 0.6819 & 0.7693 \end{matrix}],

u = [\begin{matrix} - 0.4326 & - 1.6656 & 0.1253 & 0.2877 \\ - 1.1465 & 1.1909 & 1.1892 & - 0.0376 \end{matrix}]

构造Hankel矩阵获得

u_{p} = (\begin{matrix} - 0.4326 \\ - 1.1465 \\ - 1.6656 \\ 1.1909 \end{matrix}),

u_{f} \overset{Δ}{=} (\begin{matrix} 0.1253 \\ 1.1892 \\ 0.2877 \\ - 0.0376 \end{matrix}),

y_{p} = (\begin{matrix} 0.3273 \\ - 0.5883 \\ 0.2587 \\ 1.9229 \end{matrix}),

y_{p} \overset{Δ}{=} (\begin{matrix} - 7.5755 \\ - 0.6819 \\ 3.9271 \\ 0.7693 \end{matrix})

将Hankel矩阵代入步骤2）中的状态空间模型，并转化为输出变量的子空间函数y_f＝L_ww_p+L_uu_f+L_ee_f；

4）根据步骤3）中输出变量的子空间函数获得其预测值，并将该预测值代入步骤1）中的目标函数内，通过对输入变量求导，求得目标函数最小化状态下的最优控制率

{u_{f}}^{opt} = {(R + L_{u}^{T} Q L_{u})}^{- 1} L_{u}^{T} Q (r_{k} - L_{w} w_{p});

5）将最优控制率作用于精馏系统中，闭环后计算输入变量和输出变量的方差，让扰动从某一时刻进入，即之前时刻的输入变量和扰动均为0，获得

u_{f} = - K l_{e} e_{t + 1} = (\begin{matrix} ψ_{0} \\ ψ_{1} \\ \cdot \\ \cdot \\ \cdot \\ ψ_{N - 1} \end{matrix}) e_{t + 1},

y_{f} = (I - L_{u} K) l_{e} e_{t + 1} = (\begin{matrix} γ_{0} \\ γ_{1} \\ \cdot \\ \cdot \\ \cdot \\ γ_{N - 1} \end{matrix}) e_{t + 1}

其中，

Lu = [\begin{matrix} - 0.0831 & - 0.7884 & - 0.1907 & 0.0250 \\ - 0.0075 & - 0.0710 & - 0.0172 & 0.0022 \\ 0.0431 & 0.4087 & 0.0989 & - 0.0129 \\ 0.0084 & 0.0801 & 0.0194 & - 0.0025 \end{matrix}],

le = [\begin{matrix} 0.9412 & 0.2353 \\ 0.2353 & 0.0588 \\ - 0.9412 & - 0.2353 \\ - 0.2353 & - 0.0588 \end{matrix}]

K = {(R + L_{u}^{T} Q L_{u})}^{- 1} L_{u}^{T} Q

此时求得系统输出和输入的方差关于权重矩阵的显示表达：

Var(y)＝f₁(Q,R)；Var(u)＝f₂(Q,R)。

6）令输入变量和输出变量的波动目标为输入变量和输出变量的方差，通过粒子群优化方法寻找目标函数中权衡矩阵的一次最优值，在该实施例中，粒子群优化算法的粒子个数为60，最大迭代次数为2000次，结果如表1所示；

7）以权衡矩阵的一次最优值Q1和R1为初值，利用梯度下降法对一次最优值进行进一步的优化，获得二次最优值Q2和R2，结果如表1所示，将二次最优值代入目标函数即可获得最优化的模型预测控制方法，运用这一结果来设计模型预测控制器可以让输入输出方差达到设定值。

表1

在经过两次优化后，优化结果已经达到目标值，运用本发明的方法，可以使模型预测控制器输入和输出的方差达到预定的目标，提升工业过程的经济效益。

Claims

1.一种用于精馏系统的模型预测控制方法，其特征是：步骤如下

1)建立模型预测控制器的目标函数；

2)构造精馏系统的状态空间模型；

3)利用子空间方法，根据精馏系统的输入变量和输出变量构造Hankel矩阵，并将Hankel矩阵代入步骤2)中的状态空间模型，并转化为输出变量的子空间函数；

4)根据步骤3)中输出变量的子空间函数获得其预测值，并将该预测值代入步骤1)中的目标函数内，通过对输入变量求导，求得目标函数最小化状态下的最优控制率；

5)将最优控制率作用于精馏系统中，闭环后计算输入变量和输出变量的方差，让扰动从某一时刻进入，即之前时刻的输入变量和扰动均为0，此时获得各回路输入变量和输出变量的方差；

6)令输入变量和输出变量的波动目标为输入变量和输出变量的方差，通过粒子群优化方法寻找目标函数中权衡矩阵的一次最优值；

7)以权衡矩阵的一次最优值为初值，利用梯度下降法对一次最优值进行进一步的优化，获得二次最优值，将二次最优值代入目标函数以获得最优化的模型预测控制方法。

2.根据权利要求1所述的用于精馏系统的模型预测控制方法，其特征是：步骤1)中所建立的模型预测控制器的目标函数为

J = Σ_{k = 1}^{N} [{({\hat{y}}_{k} - r_{k})}^{T} Q ({\hat{y}}_{k} - r_{k}) + u_{k}^{T} {Ru}_{k}] - - - (1)

式中，

N为数据长度；

k∈[1,2,...,N]；

r_k为人为设定的目标值，r_k为常数；

u为精馏系统的输入变量；

y为精馏系统的输出变量；

为精馏系统的输出变量的预测值；

Q为输出的权衡矩阵；

R为输入的权衡矩阵。

3.根据权利要求2所述的用于精馏系统的模型预测控制方法，其特征是：步骤2)中精馏系统的状态空间模型为

x_l＝Ax_l-1+Bu_l-1+Ee_l-1

(2)

y_l＝Cy_l-1+Du_l-1+e_l-1

式中，A、B、C、D、E为状态空间模型矩阵；

x为状态变量；

e为噪声；

l为精馏系统中某时刻的时间。

4.根据权利要求3所述的用于精馏系统的模型预测控制方法，其特征是：步骤3)利用子空间方法，根据精馏系统的输入变量和输出变量构造Hankel矩阵，

输入变量的Hankel矩阵：

u_{p} \overset{Δ}{=} (\begin{matrix} u_{0} & u_{1} & u_{2} & . . . & u_{j - 1} \\ u_{1} & u_{2} & u_{3} & . . . & u_{j} \\ . & . & . & . \\ . & . & . & . . . & . \\ . & . & . & . \\ u_{i - 1} & u_{i} & u_{i + 1} & . . . & u_{i + j - 2} \end{matrix}), u_{f} \overset{Δ}{=} (\begin{matrix} u_{i} & u_{i + 1} & u_{i + 2} & . . . & u_{i + j - 1} \\ u_{i + 1} & u_{i + 2} & u_{i + 3} & . . . & u_{i + j} \\ . & . & . & . \\ . & . & . & . . . & . \\ . & . & . & . \\ u_{2 i - 1} & u_{2 i} & u_{2 i + 1} & . . . & u_{2 i + j - 2} \end{matrix}) - - - (3 - 1)

输出变量的Hankel矩阵：

y_{p} \overset{Δ}{=} (\begin{matrix} y_{0} & y_{1} & y_{2} & . . . & y_{j - 1} \\ y_{1} & y_{2} & y_{3} & . . . & y_{j} \\ . & . & . & . \\ . & . & . & . . . & . \\ . & . & . & . \\ y_{i - 1} & y_{i} & y_{i + 1} & . . . & y_{i + j - 2} \end{matrix}), y_{f} \overset{Δ}{=} (\begin{matrix} y_{i} & y_{i + 1} & y_{i + 2} & . . . & y_{i + j - 1} \\ y_{i + 1} & y_{i + 2} & y_{i + 3} & . . . & y_{i + j} \\ . & . & . & . \\ . & . & . & . . . & . \\ . & . & . & . \\ y_{2 i - 1} & y_{2 i} & y_{2 i + 1} & . . . & y_{2 i + j - 2} \end{matrix}) - - - (3 - 2)

噪声的Hankel矩阵：

e_{p} \overset{Δ}{=} (\begin{matrix} e_{0} & e_{1} & e_{2} & . . . & e_{j - 1} \\ e_{1} & e_{2} & e_{3} & . . . & e_{j} \\ . & . & . & . \\ . & . & . & . . . & . \\ . & . & . & . \\ e_{i - 1} & e_{i} & e_{i + 1} & . . . & e_{i + j - 2} \end{matrix}), e_{f} \overset{Δ}{=} (\begin{matrix} e_{i} & e_{i + 1} & e_{i + 2} & . . . & e_{i + j - 1} \\ e_{i + 1} & e_{i + 2} & e_{i + 3} & . . . & e_{i + j} \\ . & . & . & . \\ . & . & . & . . . & . \\ . & . & . & . \\ e_{2 i - 1} & e_{2 i} & e_{2 i + 1} & . . . & e_{2 i + j - 2} \end{matrix}) - - - (3 - 3)

式中，i为Hankel矩阵的行数，同时为精馏系统工作时的时刻；

j为Hankel矩阵的列数；

定义精馏系统工作时某一时刻为t，则t也为Hankel矩阵的行数；

p代表从0到t-1时刻开始的信息；

f代表从t到t+i-1时刻开始的信息；

将输入变量的Hankel矩阵(3-1)、输出变量的Hankel矩阵(3-2)和噪声的Hankel矩阵(3-3)代入状态空间函数(2)，并转换为子空间函数为

y_f＝L_ww_p+L_uu_f+L_ee_f (4)

式中，

w_{p} = (\begin{matrix} y_{p} \\ u_{p} \end{matrix})

L_w为子空间的系数矩阵；

L_u为子空间的系数矩阵；

L_e为子空间的系数矩阵。

5.根据权利要求4所述的用于精馏系统的模型预测控制方法，其特征是：步骤(4)中，获得子空间函数式(4)的预测值为将式(4-1)代入目标函数(1)中，得到

J = {(L_{w} w_{p} + L_{u} u_{f} - r_{k})}^{T} Q (L_{w} w_{p} + L_{u} u_{f} - r_{k}) + u_{f}^{T} {Ru}_{f} - - - (5)

同时，对式(5)求u_f的偏导并置为0，在J最小化的前提下获得最优控制率

{u_{f}}^{opt} = {(R + L_{u}^{T} Q L_{u})}^{- 1} L_{u}^{T} Q (r_{k} - L_{w} w_{p}) - - - (6)

令r_k＝0，且则式(6)可简化为

\hat{e} {u_{f}}^{opt} = - K (L_{g} {\hat{u}}_{p} + L_{h} {\hat{e}}_{p}) - - - (7)

将式(7)代入式(4)中可以获得模型预测控制器作用下的输出

\begin{matrix} y_{f} = L_{g} \hat{u} + L_{h} \hat{e} + L_{u} u_{f} + L_{e} e_{f} \\ = (I - L_{u} K) L_{g} {\hat{u}}_{p} + (I - L_{u} K) L_{h} {\hat{e}}_{p} + {\hat{L}}_{e} {\hat{e}}_{f} \end{matrix} - - - (8)

式中，I为单位矩阵；

为u的预测值；

为e的预测值；

L_{g} = [\begin{matrix} g_{1} & g_{2} & . . . & g_{N - 1} & g_{N} \\ g_{2} & g_{3} & . . . & g_{N - 1} & 0 \\ . . . & . . . & . . . & . . . & . . . \\ g_{N} & 0 & 0 & . . . & 0 \end{matrix}],

g为人为定义的未知参数；

L_{h} = [\begin{matrix} l_{0} & l_{1} & . . . & l_{N - 1} & l_{N} \\ l_{1} & l_{2} & . . . & l_{N} & 0 \\ . . . & . . . & . . . & . . . & . . . \\ l_{N - 1} & 0 & 0 & . . . & 0 \end{matrix}],

l为人为定义的未知参数；

{\hat{L}}_{e} = [\begin{matrix} 0 & 0 & . . . & 0 \\ l_{0} & 0 & . . . & 0 \\ . . . & . . . & . . . & . . . \\ l_{N} & l_{N - 1} & . . . & 0 \end{matrix}],

为L_e的预期值；

{\hat{u}}_{p} = [\begin{matrix} u_{t} \\ u_{t - 1} \\ . . . \\ u_{t - N + 1} \end{matrix}],

为u_p的预期值；

{\hat{e}}_{p} = [\begin{matrix} e_{t + 1} \\ e_{t} \\ . . . \\ e_{t - N + 1} \end{matrix}],

为e_p的预期值；

{\hat{e}}_{f} = [\begin{matrix} e_{t + 2} \\ e_{t + 3} \\ . . . \\ e_{t + N + 1} \end{matrix}],

为e_f的预期值。

6.根据权利要求5所述的用于精馏系统的模型预测控制方法，其特征是：步骤(5)中，计算精馏系统的输入变量u和精馏系统的输出变量y的方差，让扰动从t+1时刻进入，即之前时刻的控制输入和扰动均为0，可得

u_{f} = - {Kl}_{e} e_{t + 1} = (\begin{matrix} ψ_{0} \\ ψ_{1} \\ . \\ . \\ . \\ ψ_{N - 1} \end{matrix}) e_{t + 1}, y_{f} = (I - L_{u} K) l_{e} e_{t + 1} = (\begin{matrix} γ_{0} \\ γ_{1} \\ . \\ . \\ . \\ γ_{N - 1} \end{matrix}) e_{t + 1}

式中，l_e为L_e的第一列，即噪声模型的脉冲响应；

(\begin{matrix} ψ_{0} \\ ψ_{1} \\ . \\ . \\ . \\ ψ_{N - 1} \end{matrix}) = - {Kl}_{e}, (\begin{matrix} γ_{0} \\ γ_{1} \\ . \\ . \\ . \\ γ_{N - 1} \end{matrix}) = (I - L_{u} K);

此时可以获得各回路输入变量u和输出变量y的方差分别为

Var [u_{t}] = Σ_{i = 0}^{N - 1} ψ_{i} Var [e_{t}] {ψ_{i}}^{T}, Var [y_{t}] = Σ_{i = 0}^{N - 1} γ_{i} Var [e_{t}] {γ_{i}}^{T} .

7.根据权利要求6所述的用于精馏系统的模型预测控制方法，其特征是：步骤(6)中，给定系统输入变量和输出变量的波动目标分别为和通过粒子群优化方法寻找Q和R的最优值Q1和R1。

8.根据权利要求7所述的用于精馏系统的模型预测控制方法，其特征是：步骤(7)中，以Q1和R1为初值，利用梯度下降法对Q和R进行进一步的优化，获得最优值Q2和R2,Q2和R2即为要求取的权重矩阵。