CN104606912A - 基于温度波特性的内部热耦合精馏在线观测器 - Google Patents

Abstract

本发明涉及一种基于温度波特性的内部热耦合精馏在线观测器，包括内部热耦合精馏塔、智能仪表、控制站、数据存储装置、上位机观测器及上位机显示界面，智能仪表与内部热耦合精馏塔连接，控制站与内部热耦合精馏塔连接，数据存储装置与智能仪表和控制站连接，上位机观测器与数据存储装置连接，上位机显示界面与上位机观测器连接，上位机观测器包括浓度梯度描述模块、温度梯度描述模块、温度波静态描述模块、温度波动态描述模块、误差计算模块及未来时刻状态预测模块。本发明采用基于内部热耦合精馏过程的温度波特性的非线性预测模型，能够高效率、高精度地跟踪内部热耦合的复杂动态特性，在线运行效率高，抗干扰能力强，实现了高效而准确的动态观测效果。

Description

基于温度波特性的内部热耦合精馏在线观测器

技术领域

本发明属于精馏节能过程中的动态建模领域，具体地说，涉及一种基于温度波特性的内部热耦合精馏在线观测器。

背景技术

精馏过程能耗占国民经济总能耗的20％，占石油化工行业的67％，是石油、化工、冶金、煤化等行业广泛使用的单元操作，与我国国民经济的诸多支柱产业息息相关。然而，精馏过程能源利用率极低，仅为5％-10％，严重制约了经济的发展。

内部热耦合精馏技术充分利用精馏段与提馏段之间的热交换，比常规精馏节能30％以上。然而，内部热耦合精馏过程的热耦合导致该过程具有显著的非线性动态特性，使得目前应用广泛的基于线性化辨识模型或者数据驱动模型的观测器，难以准确的描述该过程的复杂动态特性，往往满足不了实际生产对预测精度的需求，使得内部热耦合精馏过程的生产控制成为国际上的一个难点。因此，基于高效、准确的非线性建模方法的观测器设计，是实现该过程产品浓度准确预测的前提，也是提高内部热耦合精馏过程的产品品质的保障。

发明内容

本发明针对现有内部热耦合精馏观测器存在的测量精度不高、在线运行效率低下、对噪声敏感度高等上述不足，提供了一种基于温度波特性的内部热耦合精馏在线观测器，该观测器能够实现在线测量和实时误差校正，具有观测精度高、运行速度快、抗噪声能力强等优点。

本发明的技术方案是：一种基于温度波特性的内部热耦合精馏在线观测器，包括内部热耦合精馏塔、智能仪表、控制站、数据存储装置、上位机观测器及上位机显示界面，所述智能仪表与所述内部热耦合精馏塔连接，用于进行数据采集；所述控制站与所述内部热耦合精馏塔连接，用于实现对内部热耦合精馏塔的控制；所述数据存储装置与所述智能仪表和所述控制站连接，用于实现数据存储；所述上位机观测器与所述数据存储装置连接，用于实现在线观测并预测未来时刻内部热耦合精馏塔的状态，所述上位机显示界面与所述上位机观测器连接，用于显示内部热耦合精馏塔的状态；所述上位机观测器包括用于观测浓度及浓度梯度的浓度梯度描述模块、用于观测温度梯度的温度梯度描述模块、用于观测温度波静态的温度波静态描述模块、用于观测温度波动态的温度波动态描述模块、用于观测计算误差的误差计 算模块及用于观测未来时刻状态的未来时刻状态预测模块；其中，

(一)所述浓度梯度描述模块观测浓度及浓度梯度的步骤为：通过智能仪表中的温度检测元件、压力检测元件、流量检测元件采集相应的温度、压力、流量参数，传输至数据存储装置，再由数据存储装置传输至所述浓度梯度描述模块，通过所述浓度梯度描述模块确定出浓度及浓度梯度，所述浓度梯度描述模块包括各塔板的浓度观测和各塔板的浓度梯度观测两部分；

1)各塔板的浓度观测，根据公式(1)、(2)获得各塔板当前时刻的浓度值，并将结果传输至数据存储装置及上位机显示界面，公式(1)、(2)的表达式如下：

$X_i\left(k\right)=\frac{P_r\left(k\right)\mathrm{\α}{e}^{\frac{T_i\left(k\right)+c}{b}-a}-1}{\mathrm{\α}-1},i=\mathrm{1,2},......,f-1---\left(1\right)$

$X_i\left(k\right)=\frac{P_s\left(k\right)\mathrm{\α}{e}^{\frac{T_i\left(k\right)+c}{b}-a}-1}{\mathrm{\α}-1},i=f,f+1,......,n---\left(2\right)$

式中，k为当前采样时刻，P_r(k)为k采样时刻的精馏段压强、P_s(k)为k采样时刻的提馏段压强，T_i(k)为k采样时刻第i块塔板的温度，i表示塔板编号(i＝1,2,...,f,f+1,...,n，其中，1为塔顶编号，f为进料板编号，n为塔底编号)，P_r(k)、P_s(k)及T_i(k)由智能仪表测得，α为相对挥发度，a、b、c为安东尼常数，X_i(k)为k采样时刻第i块板塔的液相轻组分的浓度测量值；

2)各塔板的浓度梯度观测，根据公式(3)、(4)、(5)、(6)获得各塔板当前时刻的浓度梯度值，并将结果传输至温度梯度描述模块、数据存储装置及上位机显示界面，公式(3)、(4)、(5)、(6)的表达式如下：

$\frac{d{X}_1\left(k\right)}{\mathrm{dt}}=\frac{1}{H}[V_2\left(k\right)Y_2\left(k\right)-V_1\left(k\right)Y_1\left(k\right)-L_1\left(k\right)X_1\left(k\right)]---\left(3\right)$

$\frac{d{X}_i\left(k\right)}{\mathrm{dt}}=\frac{1}{H}[V_{i+1}\left(k\right)Y_{i+1}\left(k\right)-V_i\left(k\right)Y_i\left(k\right)+L_{i-1}\left(k\right)X_{i-1}\left(k\right)-L_i\left(k\right)X_i\left(k\right)]---\left(4\right)$

(i＝2,...,n-1并且i≠f)

$\frac{d{X}_f\left(k\right)}{\mathrm{dt}}=\frac{1}{H}[V_{f+1}\left(k\right)Y_{f+1}\left(k\right)-V_f\left(k\right)Y_f\left(k\right)+L_{f-1}\left(k\right)X_{f-1}\left(k\right)-L_f\left(k\right)X_f\left(k\right)+F\left(k\right)Z_f\left(k\right)]---\left(5\right)$

$\frac{d{X}_n\left(k\right)}{\mathrm{dt}}=\frac{1}{H}[-V_n\left(k\right)Y_n\left(k\right)+L_{n-1}\left(k\right)X_{n-1}\left(k\right)-L_n\left(k\right)X_n\left(k\right)]---\left(6\right)$

式中，H为持液量，V_i(k)为k采样时刻第i块板塔的气相流率，L_i(k)为k采样时刻第i块板塔的液相流率，X_i(k)为k采样时刻第i块板塔的液相轻组分的浓度测量值，Y_i(k)为k采样时 刻第i块板塔的气相轻组分浓度，为k采样时刻第i块板塔的液相轻组分的浓度梯度值，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号，n为塔底编号)，F(k)为k采样时刻的进料流量，Z_f(k)为k采样时刻的进料组分；所述Y_i(k)由公式(7)得到，公式(7)的表达式如下：

Y_i(k)＝αX_i(k)/[(α-1)X_i(k)+1] i＝1,2,...,f,f+1,...,n    (7) 

式中，α为相对挥发度，X_i(k)为k采样时刻第i块板塔的液相轻组分的浓度测量值；

(二)所述温度梯度描述模块观测温度梯度的步骤为：通过数据存储装置提取智能仪表中的温度检测元件、压力检测元件收集的温度、压力参数，以及所述浓度梯度描述模块得到的浓度梯度信息，通过所述温度梯度描述模块确定温度梯度，所述温度梯度描述模块根据公式(8)、(9)获得温度梯度，公式(8)、(9)的表达式如下：

$\frac{d{T}_i\left(k\right)}{\mathrm{dt}}=-\frac{d{X}_i\left(k\right)}{\mathrm{dt}}\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_r\left(k\right)\mathrm{\αb}{e}^{\frac{b}{T_i\left(k\right)+c}-a}},i=\mathrm{1,2},......,f-1---\left(8\right)$

$\frac{d{T}_i\left(k\right)}{\mathrm{dt}}=-\frac{d{X}_i\left(k\right)}{\mathrm{dt}}\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_s\left(k\right)\mathrm{\αb}{e}^{\frac{b}{T_i\left(k\right)+c}-a}},i=f,f+1,......,n---\left(9\right)$

式中，为k采样时刻第i块板塔的液相轻组分的浓度梯度值，α为相对挥发度，a、b、c为安东尼常数，T_i(k)为k采样时刻第i块塔板的温度，P_r(k)为k采样时刻的精馏段压强、P_s(k)为k采样时刻的提馏段压强，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号，n为塔底编号)，为k采样时刻第i块板塔的温度梯度值；

(三)所述温度波静态描述模块观测温度波静态步骤为：根据内部热耦合精馏塔的温度波特性，通过温度检测元件采集各塔板当前时刻的温度测量值，并将各温度值在坐标轴中对应的点连成连续的光滑曲线，获得当前时刻精馏段与提馏段的温度波形，进而得到温度波静态描述函数的常系数T_r1、T_r2、T_s1、T_s2、γ_r、γ_s的值，以及精馏段、提馏段温度波的拐点初始值S_r(0)、S_s(0)，通过温度波静态描述函数公式(10)、(11)获得各塔板的温度预测初始值公式(10)、(11)的表达式如下：

${\hat{T}}_i\left(k\right)=\frac{T_{r1}{e}^{-{\mathrm{\γ}}_r(i-S_r\left(k\right))}+T_{r2}}{1+e^{-{\mathrm{\γ}}_r(i-S_r\left(k\right))}},i=\mathrm{1,2},...,f-1---\left(10\right)$

${\hat{T}}_i\left(k\right)=\frac{T_{s1}{e}^{-{\mathrm{\γ}}_s(i-S_s\left(k\right))}+T_{s2}}{1+e^{-{\mathrm{\γ}}_s(i-S_s\left(k\right))}},i=f,f+1,......,n---\left(11\right)$

式中，k为当前采样时刻，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号,n为塔底编号)，为k采样时刻第i块板塔的温度预测值，S_r(k)、S_s(k)分别为k采样时刻内部热耦合精馏塔精馏段、提馏段温度波的拐点，T_r1、T_r2、T_s1、T_s2、γ_r、γ_s为温度波静态描述函数的常系数，T_r1、T_r2、T_s1、T_s2分别表示精馏段与提馏段温度波两端的渐进浓度，γ_r、γ_s分别表征精馏段与提馏段温度波拐点处的斜率大小；

(四)所述温度波动态描述模块观测温度波动态的步骤为：根据所述温度梯度描述模块和所述温度波静态描述模块所获得的信息，确定各塔板温度波在未来时刻的变化趋势，该变化趋势用拐点的移动速度来表示，具体由公式(12)、(13)获得，公式(12)、(13)的表达式如下：

$\frac{d{S}_r}{\mathrm{dt}}\left(k\right)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}\left[\frac{d{T}_i\left(k\right)}{\mathrm{dt}}\frac{T_{r2}-T_{r1}}{{\mathrm{\γ}}_r(T_i\left(k\right)-T_{r2})(T_i\left(k\right)-T_{r1})}\right]---\left(12\right)$

$\frac{d{S}_s}{\mathrm{dt}}\left(k\right)=\underset{i=f}{\overset{n}{\mathrm{\Σ}}}\left[\frac{d{T}_i\left(k\right)}{\mathrm{dt}}\frac{T_{s2}-T_{s1}}{{\mathrm{\γ}}_s(T_i\left(k\right)-T_{s2})(T_i\left(k\right)-T_{s1})}\right]---\left(13\right)$

式中，分别为k采样时刻精馏段和提馏段温度波的拐点移动速度，为k采样时刻第i块板塔的温度梯度值，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号,n为塔底编号)，T_i(k)为k采样时刻第i块塔板的温度，T_r1、T_r2、T_s1、T_s2、γ_r、γ_s为温度波静态描述函数的常系数；

(五)所述误差计算模块用观测计算误差的步骤为：根据公式(14)计算k采用时刻的温度预测误差，公式(14)的表达式如下：

$e_i\left(k\right)=T_i\left(k\right)-{\hat{T}}_i\left(k\right)---\left(14\right)$

式中，e_i(k)为k采样时刻的温度预测误差，T_i(k)为k采样时刻第i块塔板的温度测量值，为k采样时刻第i块塔板的温度预测值；

(六)所述未来时刻状态模块其观测未来时刻状态的步骤为：根据当前时刻的参数信息预测未来时刻的状态，包括以下三部分：

1)根据公式(15)、(16)、(17)、(18)、(19)(20)获得下一时刻拐点位置以及温度和组分浓度，公式(15)、(16)、(17)、(18)、(19)(20)的表达式如下：

$S_r(k+1)=S_r\left(k\right)+\frac{d{S}_r}{\mathrm{dt}}\left(k\right)\×T---\left(15\right)$

$S_s(k+1)=S_s\left(k\right)+\frac{d{S}_s}{\mathrm{dt}}\left(k\right)\×T---\left(16\right)$

${\hat{T}}_i(k+1)=\frac{T_{r1}{e}^{-{\mathrm{\γ}}_r(i-S_r(k+1))}+T_{r2}}{1+e^{-{\mathrm{\γ}}_r(i-S_r(k+1))}},i=\mathrm{1,2},...,f-1---\left(17\right)$

${\hat{T}}_i(k+1)=\frac{T_{s1}{e}^{-{\mathrm{\γ}}_s(i-S_s(k+1))}+T_{s2}}{1+e^{-{\mathrm{\γ}}_s(i-S_s(k+1))}},i=f,f+1,......,n---\left(18\right)$

${\hat{X}}_i(k+1)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_r\left(k\right)e^{\frac{{\hat{T}}_i(k+1)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=\mathrm{1,2},......,f-1---\left(19\right)$

${\hat{X}}_i(k+1)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_s\left(k\right)e^{\frac{{\hat{T}}_i(k+1)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=f,f+1,......,n---\left(20\right)$

式中，T为采样周期，i表示塔板编号(i＝1,2,...,f,f+1,...,n，其中，1为塔顶编号，f为进料板编号，n为塔底编号)，α为相对挥发度，a、b、c为安东尼常数，P_r(k)为k采样时刻的精馏段压强、P_s(k)为k采样时刻的提馏段压强，S_r(k+1)、S_s(k+1)分别为k+1采样时刻内部热耦合精馏塔精馏段、提馏段温度波的拐点，为k+1采样时刻第i块塔板的温度预测值，为k+1采样时刻第i块塔板的液相轻组分浓度的预测值；

2)修正温度与组分浓度预测值，根据所述误差计算模块得到的温度预测误差，由公式(21)、(22)、(23)获得修正后的温度与组分浓度的预测值，公式(21)、(22)、(23)的表达式如下：

${\tilde{T}}_i(k+1)={\hat{T}}_i(k+1)+e_i\left(k\right)---\left(21\right)$

${\tilde{X}}_i(k+1)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_r\left(k\right)e^{\frac{{\tilde{T}}_i(k+1)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=\mathrm{1,2},......,f-1---\left(22\right)$

${\tilde{X}}_i(k+1)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_s\left(k\right)e^{\frac{{\tilde{T}}_i(k+1)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=f,f+1,......,n---\left(23\right)$

式中，i表示塔板编号(i＝1,2,...,f,f+1,...,n，其中，1为塔顶编号，f为进料板编号，n为塔底编号)，α为相对挥发度，a、b、c为安东尼常数，P_r(k)为k采样时刻的精馏段压强、P_s(k)为k采样时刻的提馏段压强，表示修正后的k+1采样时刻第i块塔板的温度预测值， 表示修正后的k+1采样时刻第i块塔板的液相轻组分浓度的预测值；

3)预测未来m个时刻状态预测，即对k+2到k+m(m>2)采样时刻的温度和组分浓度进行预测，通过以下步骤完成：

A、对k+2采样时刻的温度和组分进行预测，首先由公式(24)、(25)得到k+1采样时刻的拐点的移动速度，公式(24)、(25)的表达式如下：

$\frac{d{S}_r}{\mathrm{dt}}(k+1)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}\left[\frac{{\tilde{T}}_i(k+1)-T_i\left(k\right)}{T}\frac{T_{r2}-T_{r1}}{{\mathrm{\γ}}_r({\tilde{T}}_i(k+1)-T_{r2})({\tilde{T}}_i(k+1)-T_{r1})}\right]---\left(24\right)$

$\frac{d{S}_s}{\mathrm{dt}}(k+1)=\underset{i=f}{\overset{n}{\mathrm{\Σ}}}\left[\frac{{\tilde{T}}_i(k+1)-T_i\left(k\right)}{T}\frac{T_{s2}-T_{s1}}{{\mathrm{\γ}}_s({\tilde{T}}_i(k+1)-T_{s2})({\tilde{T}}_i(k+1)-T_{s1})}\right]---\left(25\right)$

进而由公式(26)、(27)得到k+2采样时刻的拐点位置，公式(26)、(27)的表达式如下：

$S_r(k+2)=S_r(k+1)+\frac{d{S}_r}{\mathrm{dt}}(k+1)\×T---\left(26\right)$

$S_s(k+2)=S_s(k+1)+\frac{d{S}_s}{\mathrm{dt}}(k+1)\×T---\left(27\right)$

由于所述内部热耦合精馏塔在一个生产周期内提馏段压强变化很小，则认定为P_s(k+m-1)＝...＝P_s(k+1)＝P_s(k)，精馏段压强是主要的操作变量，由所述控制站的控制策略决定，所述控制站将未来时刻的操作变量值传输至所述数据存储装置，再由所述数据存储装置传输给所述上位机观测器，进而得到未来m-1个时刻的精馏段压强P_r(k+1),P_r(k+2),...,P_r(k+m-1)；根据公式(28)、(29)、(30)、(31)得到k+2采用时刻的温度和组分浓度预测，公式(28)、(29)、(30)、(31)的表达式如下：

${\hat{T}}_i(k+1)=\frac{T_{r1}{e}^{-{\mathrm{\γ}}_r(i-S_r(k+1))}+T_{r2}}{1+e^{-{\mathrm{\γ}}_r(i-S_r(k+1))}},i=\mathrm{1,2},...,f-1---\left(17\right)$

${\hat{T}}_i(k+2)=\frac{T_{s1}{e}^{-{\mathrm{\γ}}_s(i-S_s(k+2))}+T_{s2}}{1+e^{-{\mathrm{\γ}}_s(i-S_s(k+2))}},i=f,f+1,......,n---\left(29\right)$

${\hat{X}}_i(k+2)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_r(k+1)e^{\frac{{\hat{T}}_i(k+2)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=\mathrm{1,2},......,f-1---\left(30\right)$

${\hat{X}}_i(k+2)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_s\left(k\right)e^{\frac{{\hat{T}}_i(k+2)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=f,f+1,......,n---\left(31\right)$

设定误差不发生变化，则由公式(32)、(33)、(34)获得修正后的温度与组分浓度的预测值，公式(32)、(33)、(34)的表达式如下：

${\tilde{T}}_i(k+2)={\hat{T}}_i(k+2)+e_i\left(k\right)---\left(32\right)$

${\tilde{X}}_i(k+2)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_r(k+1)e^{\frac{{\tilde{T}}_i(k+2)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=1,2,......,f-1---\left(33\right)$

${\tilde{X}}_i(k+2)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_s\left(k\right)e^{\frac{{\tilde{T}}_i(k+2)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=f,f+1,......,n---\left(34\right)$

B、通过步骤A中所述的方法，进而可以递推求出k+3到k+m采样时刻的温度和组分浓 度预测值，由公式(35)、(36)、(37)获得第k+m采样时刻的温度与组分浓度的预测值，公式(35)、(36)、(37)的表达式如下：

${\tilde{T}}_i(k+m)={\hat{T}}_i(k+m)+e_i\left(k\right)---\left(35\right)$

${\tilde{X}}_i(k+m)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_r(k+m-1)e^{\frac{{\tilde{T}}_i(k+m)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=1,2,......,f-1---\left(36\right)$

${\tilde{X}}_i(k+m)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_s\left(k\right)e^{\frac{{\tilde{T}}_i(k+m)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=f,f+1,......,n---\left(37\right)$

式中，i表示塔板编号(i＝1,2,...,f,f+1,...,n，其中，1为塔顶编号，f为进料板编号，n为塔底编号)，α为相对挥发度，a、b、c为安东尼常数，表示修正后的k+m采样时刻第i块塔板的温度预测值，表示修正后的k+m采样时刻第i块塔板的液相轻组分浓度的预测值。

作为优选，所述上位机观测器还用于设定采用周期T以及预测时域m，并将上位机观测器输出的当前时刻温度、浓度测量值，以及未来m个时刻的温度、浓度预测值，传输到上位机显示界面，并将以上值传输给数据存储装置，方便操作人员进行生产质量控制，提高生产控制品质。

本发明的有益效果是：本发明对内部热耦合精馏过程中重要的温度、浓度状态变量进行在线观测，克服已有的在先观测器测量精度不高、在线运行效率低下、对噪声敏感度低的不足，通过对温度波特性的准确把握，建立内部热耦合精馏过程的高效率、高精度非线性模型，实现对温度、浓度进行高效率、高精度的在线观测。本发明采用基于内部热耦合精馏过程的温度波特性的非线性预测模型，能够高效率、高精度地跟踪内部热耦合的复杂动态特性，在线运行效率高，预测比传统基于数据的统计模型以及基于线性化的辨识模型精确，抗干扰能力强，实现了高效而准确的动态观测效果，实时误差滚动校正保证在有噪声的影响下，观测器仍能正常工作。

附图说明

附图1为本发明具体实施例的基本结构示意图。

附图2为本发明具体实施例上位机观测器的原理图。

具体实施方式

以下结合附图对本发明作进一步说明。

如图1和图2所示，一种基于温度波特性的内部热耦合精馏在线观测器，包括内部热耦合精馏塔1、智能仪表2、控制站3、数据存储装置4、上位机观测器5及上位机显示界面6，所述智能仪表2与所述内部热耦合精馏塔1连接，用于进行数据采集；所述控制站3与所述内部热耦合精馏塔1连接，用于实现对内部热耦合精馏塔1的控制；所述数据存储装置4与所述智能仪表2和所述控制站3连接，用于实现数据存储；所述上位机观测器5与所述数据存储装置4连接，用于实现在线观测并预测未来时刻内部热耦合精馏塔1的状态，所述上位机显示界面6与所述上位机观测器5连接，用于显示内部热耦合精馏塔1的状态；所述上位机观测器5包括用于观测浓度及浓度梯度的浓度梯度描述模块7、用于观测温度梯度的温度梯度描述模块8、用于观测温度波静态的温度波静态描述模块9、用于观测温度波动态的温度波动态描述模块10、用于观测计算误差的误差计算模块11及用于观测未来时刻状态的未来时刻状态预测模块12。

本实施例所述上位机观测器5在进行观测时，其中：

(一)所述浓度梯度描述模块7观测浓度及浓度梯度的步骤为：通过智能仪表2中的温度检测元件、压力检测元件、流量检测元件采集相应的温度、压力、流量参数，传输至数据存储装置4，再由数据存储装置4传输至所述浓度梯度描述模块7，通过所述浓度梯度描述模块7确定出浓度及浓度梯度，所述浓度梯度描述模块7包括各塔板的浓度观测和各塔板的浓度梯度观测两部分；

1)各塔板的浓度观测，根据公式(1)、(2)获得各塔板当前时刻的浓度值，并将结果传输至数据存储装置及上位机显示界面，公式(1)、(2)的表达式如下：

$X_i\left(k\right)=\frac{P_r\left(k\right)\mathrm{\α}{e}^{\frac{T_i\left(k\right)+c}{b}-a}-1}{\mathrm{\α}-1},i=\mathrm{1,2},......,f-1---\left(1\right)$

$X_i\left(k\right)=\frac{P_s\left(k\right)\mathrm{\α}{e}^{\frac{T_i\left(k\right)+c}{b}-a}-1}{\mathrm{\α}-1},i=f,f+1,......,n---\left(2\right)$

式中，k为当前采样时刻，P_r(k)为k采样时刻的精馏段压强、P_s(k)为k采样时刻的提馏段压强，T_i(k)为k采样时刻第i块塔板的温度，i表示塔板编号(i＝1,2,...,f,f+1,...,n，其中，1为塔顶编号，f为进料板编号，n为塔底编号)，P_r(k)、P_s(k)及T_i(k)由智能仪表测得，α为相对挥发度，a、b、c为安东尼常数，X_i(k)为k采样时刻第i块板塔的液相轻组分的浓度测量值；

2)各塔板的浓度梯度观测，根据公式(3)、(4)、(5)、(6)获得各塔板当前时刻的浓度 梯度值，并将结果传输至温度梯度描述模块、数据存储装置及上位机显示界面，公式(3)、(4)、(5)、(6)的表达式如下：

$\frac{d{X}_1\left(k\right)}{\mathrm{dt}}=\frac{1}{H}[V_2\left(k\right)Y_2\left(k\right)-V_1\left(k\right)Y_1\left(k\right)-L_1\left(k\right)X_1\left(k\right)]---\left(3\right)$

$\frac{d{X}_i\left(k\right)}{\mathrm{dt}}=\frac{1}{H}[V_{i+1}\left(k\right)Y_{i+1}\left(k\right)-V_i\left(k\right)Y_i\left(k\right)+L_{i-1}\left(k\right)X_{i-1}\left(k\right)-L_i\left(k\right)X_i\left(k\right)]---\left(4\right)$

(i＝2,...,n-1并且i≠f)

$\frac{d{X}_f\left(k\right)}{\mathrm{dt}}=\frac{1}{H}[V_{f+1}\left(k\right)Y_{f+1}\left(k\right)-V_f\left(k\right)Y_f\left(k\right)+L_{f-1}\left(k\right)X_{f-1}\left(k\right)-L_f\left(k\right)X_f\left(k\right)+F\left(k\right)Z_f\left(k\right)]---\left(5\right)$

$\frac{d{X}_n\left(k\right)}{\mathrm{dt}}=\frac{1}{H}[-V_n\left(k\right)Y_n\left(k\right)+L_{n-1}\left(k\right)X_{n-1}\left(k\right)-L_n\left(k\right)X_n\left(k\right)]---\left(6\right)$

式中，H为持液量，V_i(k)为k采样时刻第i块板塔的气相流率，L_i(k)为k采样时刻第i块板塔的液相流率，X_i(k)为k采样时刻第i块板塔的液相轻组分的浓度测量值，Y_i(k)为k采样时刻第i块板塔的气相轻组分浓度，为k采样时刻第i块板塔的液相轻组分的浓度梯度值，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号，n为塔底编号)，F(k)为k采样时刻的进料流量，Z_f(k)为k采样时刻的进料组分；所述Y_i(k)由公式(7)得到，公式(7)的表达式如下：

Y_i(k)＝αX_i(k)/[(α-1)X_i(k)+1]  i＝1,2,...,f,f+1,...,n    (7) 

式中，α为相对挥发度，X_i(k)为k采样时刻第i块板塔的液相轻组分的浓度测量值。

(二)所述温度梯度描述模块8观测温度梯度的步骤为：通过数据存储装置4提取智能仪表2中的温度检测元件、压力检测元件收集的温度、压力参数，以及所述浓度梯度描述模块7得到的浓度梯度信息，通过所述温度梯度描述模块8确定温度梯度，所述温度梯度描述模块8根据公式(8)、(9)获得温度梯度，公式(8)、(9)的表达式如下：

$\frac{d{T}_i\left(k\right)}{\mathrm{dt}}=-\frac{d{X}_i\left(k\right)}{\mathrm{dt}}\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_r\left(k\right)\mathrm{\αb}{e}^{\frac{b}{T_i\left(k\right)+c}-a}},i=\mathrm{1,2},......,f-1---\left(8\right)$

$\frac{d{T}_i\left(k\right)}{\mathrm{dt}}=-\frac{d{X}_i\left(k\right)}{\mathrm{dt}}\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_s\left(k\right)\mathrm{\αb}{e}^{\frac{b}{T_i\left(k\right)+c}-a}},i=f,f+1,......,n---\left(9\right)$

式中，为k采样时刻第i块板塔的液相轻组分的浓度梯度值，α为相对挥发度，a、b、c为安东尼常数，T_i(k)为k采样时刻第i块塔板的温度，P_r(k)为k采样时刻的精馏段压强、 P_s(k)为k采样时刻的提馏段压强，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号，n为塔底编号)，为k采样时刻第i块板塔的温度梯度值。

(三)所述温度波静态描述模块9观测温度波静态步骤为：根据内部热耦合精馏塔1的温度波特性，通过温度检测元件采集各塔板当前时刻的温度测量值，并将各温度值在坐标轴中对应的点连成连续的光滑曲线，获得当前时刻精馏段与提馏段的温度波形，进而得到温度波静态描述函数的常系数T_r1、T_r2、T_s1、T_s2、γ_r、γ_s的值，以及精馏段、提馏段温度波的拐点初始值S_r(0)、S_s(0)，通过温度波静态描述函数公式(10)、(11)获得各塔板的温度预测初始值公式(10)、(11)的表达式如下：

${\hat{T}}_i\left(k\right)=\frac{T_{r1}{e}^{-{\mathrm{\γ}}_r(i-S_r\left(k\right))}+T_{r2}}{1+e^{-{\mathrm{\γ}}_r(i-S_r\left(k\right))}},i=\mathrm{1,2},...,f-1---\left(10\right)$

${\hat{T}}_i\left(k\right)=\frac{T_{s1}{e}^{-{\mathrm{\γ}}_s(i-S_s\left(k\right))}+T_{s2}}{1+e^{-{\mathrm{\γ}}_s(i-S_s\left(k\right))}},i=f,f+1,......,n---\left(11\right)$

式中，k为当前采样时刻，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号,n为塔底编号)，为k采样时刻第i块板塔的温度预测值，S_r(k)、S_s(k)分别为k采样时刻内部热耦合精馏塔精馏段、提馏段温度波的拐点，T_r1、T_r2、T_s1、T_s2、γ_r、γ_s为温度波静态描述函数的常系数，T_r1、T_r2、T_s1、T_s2分别表示精馏段与提馏段温度波两端的渐进浓度，γ_r、γ_s分别表征精馏段与提馏段温度波拐点处的斜率大小。

(四)所述温度波动态描述模块10观测温度波动态的步骤为：根据所述温度梯度描述模块8和所述温度波静态描述模块9所获得的信息，确定各塔板温度波在未来时刻的变化趋势，该变化趋势用拐点的移动速度来表示，具体由公式(12)、(13)获得，公式(12)、(13)的表达式如下：

$\frac{d{S}_r}{\mathrm{dt}}\left(k\right)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}\left[\frac{d{T}_i\left(k\right)}{\mathrm{dt}}\frac{T_{r2}-T_{r1}}{{\mathrm{\γ}}_r(T_i\left(k\right)-T_{r2})(T_i\left(k\right)-T_{r1})}\right]---\left(12\right)$

$\frac{d{S}_s}{\mathrm{dt}}\left(k\right)=\underset{i=f}{\overset{n}{\mathrm{\Σ}}}\left[\frac{d{T}_i\left(k\right)}{\mathrm{dt}}\frac{T_{s2}-T_{s1}}{{\mathrm{\γ}}_s(T_i\left(k\right)-T_{s2})(T_i\left(k\right)-T_{s1})}\right]---\left(13\right)$

式中，分别为k采样时刻精馏段和提馏段温度波的拐点移动速度，为k采样时刻第i块板塔的温度梯度值，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号,n为塔底编号)，T_i(k)为k采样时刻第i块塔板的温度，T_r1、T_r2、T_s1、T_s2、 γ_r、γ_s为温度波静态描述函数的常系数。

(五)所述误差计算模块11用观测计算误差的步骤为：根据公式(14)计算k采用时刻的温度预测误差，公式(14)的表达式如下：

$e_i\left(k\right)=T_i\left(k\right)-{\hat{T}}_i\left(k\right)---\left(14\right)$

式中，e_i(k)为k采样时刻的温度预测误差，T_i(k)为k采样时刻第i块塔板的温度测量值，为k采样时刻第i块塔板的温度预测值。

(六)所述未来时刻状态模块12其观测未来时刻状体的步骤为：根据当前时刻的参数信息预测未来时刻的状态，包括以下三部分：

1)根据公式(15)、(16)、(17)、(18)、(19)(20)获得下一时刻拐点位置以及温度和组分浓度，公式(15)、(16)、(17)、(18)、(19)(20)的表达式如下：

$S_r(k+1)=S_r\left(k\right)+\frac{d{S}_r}{\mathrm{dt}}\left(k\right)\×T---\left(15\right)$

$S_s(k+1)=S_s\left(k\right)+\frac{d{S}_s}{\mathrm{dt}}\left(k\right)\×T---\left(16\right)$

${\hat{T}}_i(k+1)=\frac{T_{r1}{e}^{-{\mathrm{\γ}}_r(i-S_r(k+1))}+T_{r2}}{1+e^{-{\mathrm{\γ}}_r(i-S_r(k+1))}},i=\mathrm{1,2},...,f-1---\left(17\right)$

${\hat{T}}_i(k+1)=\frac{T_{s1}{e}^{-{\mathrm{\γ}}_s(i-S_s(k+1))}+T_{s2}}{1+e^{-{\mathrm{\γ}}_s(i-S_s(k+1))}},i=f,f+1,......,n---\left(18\right)$

${\hat{X}}_i(k+1)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_r\left(k\right)e^{\frac{{\hat{T}}_i(k+1)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=\mathrm{1,2},......,f-1---\left(19\right)$

${\hat{X}}_i(k+1)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_s\left(k\right)e^{\frac{{\hat{T}}_i(k+1)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=f,f+1,......,n---\left(20\right)$

式中，T为采样周期，i表示塔板编号(i＝1,2,...,f,f+1,...,n，其中，1为塔顶编号，f为进料板编号，n为塔底编号)，α为相对挥发度，a、b、c为安东尼常数，P_r(k)为k采样时刻的精馏段压强、P_s(k)为k采样时刻的提馏段压强，S_r(k+1)、S_s(k+1)分别为k+1采样时刻内部热耦合精馏塔精馏段、提馏段温度波的拐点，为k+1采样时刻第i块塔板的温度预测值，为k+1采样时刻第i块塔板的液相轻组分浓度的预测值；

2)修正温度与组分浓度预测值，根据所述误差计算模块11得到的温度预测误差，由公式(21)、(22)、(23)获得修正后的温度与组分浓度的预测值，公式(21)、(22)、(23)的表达式如下：

${\tilde{T}}_i(k+1)={\hat{T}}_i(k+1)+e_i\left(k\right)---\left(21\right)$

${\tilde{X}}_i(k+1)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_r\left(k\right)e^{\frac{{\tilde{T}}_i(k+1)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=\mathrm{1,2},......,f-1---\left(22\right)$

${\tilde{X}}_i(k+1)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_s\left(k\right)e^{\frac{{\tilde{T}}_i(k+1)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=f,f+1,......,n---\left(23\right)$

式中，i表示塔板编号(i＝1,2,...,f,f+1,...,n，其中，1为塔顶编号，f为进料板编号，n为塔底编号)，α为相对挥发度，a、b、c为安东尼常数，P_r(k)为k采样时刻的精馏段压强、P_s(k)为k采样时刻的提馏段压强，表示修正后的k+1采样时刻第i块塔板的温度预测值， 表示修正后的k+1采样时刻第i块塔板的液相轻组分浓度的预测值；

3)预测未来m个时刻状态预测，即对k+2到k+m(m>2)采样时刻的温度和组分浓度进行预测，通过以下步骤完成：

B、对k+2采样时刻的温度和组分进行预测，首先由公式(24)、(25)得到k+1采样时刻的拐点的移动速度，公式(24)、(25)的表达式如下：

$\frac{d{S}_r}{\mathrm{dt}}(k+1)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}\left[\frac{{\tilde{T}}_i(k+1)-T_i\left(k\right)}{T}\frac{T_{r2}-T_{r1}}{{\mathrm{\γ}}_r({\tilde{T}}_i(k+1)-T_{r2})({\tilde{T}}_i(k+1)-T_{r1})}\right]---\left(24\right)$

$\frac{d{S}_s}{\mathrm{dt}}(k+1)=\underset{i=f}{\overset{n}{\mathrm{\Σ}}}\left[\frac{{\tilde{T}}_i(k+1)-T_i\left(k\right)}{T}\frac{T_{s2}-T_{s1}}{{\mathrm{\γ}}_s({\tilde{T}}_i(k+1)-T_{s2})({\tilde{T}}_i(k+1)-T_{s1})}\right]---\left(25\right)$

进而由公式(26)、(27)得到k+2采样时刻的拐点位置，公式(26)、(27)的表达式如下：

$S_r(k+2)=S_r(k+1)+\frac{d{S}_r}{\mathrm{dt}}(k+1)\×T---\left(26\right)$

$S_s(k+2)=S_s(k+1)+\frac{d{S}_s}{\mathrm{dt}}(k+1)\×T---\left(27\right)$

由于所述内部热耦合精馏塔在一个生产周期内提馏段压强变化很小，则认定为P_s(k+m-1)＝...＝P_s(k+1)＝P_s(k)，精馏段压强是主要的操作变量，由所述控制站3的控制策略决定，所述控制站3将未来时刻的操作变量值传输至所述数据存储装置4，再由所述数据存储装置4传输给所述上位机观测器5，进而得到未来m-1个时刻的精馏段压强P_r(k+1),P_r(k+2),...,P_r(k+m-1)；根据公式(28)、(29)、(30)、(31)得到k+2采用时刻的温度和组分浓度预测，公式(28)、(29)、(30)、(31)的表达式如下：

${\hat{T}}_i(k+2)=\frac{T_{r1}{e}^{-{\mathrm{\γ}}_r(i-S_r(k+2))}+T_{r2}}{1+e^{-{\mathrm{\γ}}_r(i-S_r(k+2))}},i=\mathrm{1,2},...,f-1---\left(28\right)$

${\hat{T}}_i(k+1)=\frac{T_{s1}{e}^{-{\mathrm{\γ}}_s(i-S_s(k+1))}+T_{s2}}{1+e^{-{\mathrm{\γ}}_s(i-S_s(k+1))}},i=f,f+1,......,n---\left(18\right)$

${\hat{X}}_i(k+2)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_r(k+1)e^{\frac{{\hat{T}}_i(k+2)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=\mathrm{1,2},......,f-1---\left(30\right)$

${\hat{X}}_i(k+2)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_s\left(k\right)e^{\frac{{\hat{T}}_i(k+2)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=f,f+1,......,n---\left(31\right)$

设定误差不发生变化，则由公式(32)、(33)、(34)获得修正后的温度与组分浓度的预测值，公式(32)、(33)、(34)的表达式如下：

${\tilde{T}}_i(k+2)={\hat{T}}_i(k+2)+e_i\left(k\right)---\left(32\right)$

${\tilde{X}}_i(k+2)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_r(k+1)e^{\frac{{\tilde{T}}_i(k+2)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=1,2,......,f-1---\left(33\right)$

${\tilde{X}}_i(k+2)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_s\left(k\right)e^{\frac{{\tilde{T}}_i(k+2)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=f,f+1,......,n---\left(34\right)$

B、通过步骤A中所述的方法，进而可以递推求出k+3到k+m采样时刻的温度和组分浓度预测值，由公式(35)、(36)、(37)获得第k+m采样时刻的温度与组分浓度的预测值，公式(35)、(36)、(37)的表达式如下：

${\tilde{T}}_i(k+m)={\hat{T}}_i(k+m)+e_i\left(k\right)---\left(35\right)$

${\tilde{X}}_i(k+m)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_r(k+m-1)e^{\frac{{\tilde{T}}_i(k+m)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=1,2,......,f-1---\left(36\right)$

${\tilde{X}}_i(k+m)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_s\left(k\right)e^{\frac{{\tilde{T}}_i(k+m)+c}{b}-a}-\frac{1}{\mathrm{\α}}),i=f,f+1,......,n---\left(37\right)$

式中，i表示塔板编号(i＝1,2,...,f,f+1,...,n，其中，1为塔顶编号，f为进料板编号，n为塔底编号)，α为相对挥发度，a、b、c为安东尼常数，表示修正后的k+m采样时刻第i块塔板的温度预测值，表示修正后的k+m采样时刻第i块塔板的液相轻组分浓度的预测值。

本实施例中，所述上位机观测器5还用于设定采用周期T以及预测时域m，并将上位机观测器5输出的当前时刻温度、浓度测量值，以及未来m个时刻的温度、浓度预测值，传输到上位机显示界面6，并将以上值传输给数据存储装置4。

上述实施例用来解释本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明做出的任何修改和改变，都落入本发明的保护范围。

Claims (2)

1.一种基于温度波特性的内部热耦合精馏在线观测器，其特征在于：包括内部热耦合精馏塔、智能仪表、控制站、数据存储装置、上位机观测器及上位机显示界面，所述智能仪表与所述内部热耦合精馏塔连接，用于进行数据采集；所述控制站与所述内部热耦合精馏塔连接，用于实现对内部热耦合精馏塔的控制；所述数据存储装置与所述智能仪表和所述控制站连接，用于实现数据存储；所述上位机观测器与所述数据存储装置连接，用于实现在线观测并预测未来时刻内部热耦合精馏塔的状态，所述上位机显示界面与所述上位机观测器连接，用于显示内部热耦合精馏塔的状态；所述上位机观测器包括用于观测浓度及浓度梯度的浓度梯度描述模块、用于观测温度梯度的温度梯度描述模块、用于观测温度波静态的温度波静态描述模块、用于观测温度波动态的温度波动态描述模块、用于观测计算误差的误差计算模块及用于观测未来时刻状态的未来时刻状态预测模块；其中，

(一)所述浓度梯度描述模块观测浓度及浓度梯度的步骤为：通过智能仪表中的温度检测元件、压力检测元件、流量检测元件采集相应的温度、压力、流量参数，传输至数据存储装置，再由数据存储装置传输至所述浓度梯度描述模块，通过所述浓度梯度描述模块确定出浓度及浓度梯度，所述浓度梯度描述模块包括各塔板的浓度观测和各塔板的浓度梯度观测两部分；

1)各塔板的浓度观测，根据公式(1)、(2)获得各塔板当前时刻的浓度值，并将结果传输至数据存储装置及上位机显示界面，公式(1)、(2)的表达式如下：

$X_i\left(k\right)=\frac{P_r\left(k\right)\mathrm{\α}{e}^{\frac{T_i\left(k\right)+c}{b}-a}-1}{\mathrm{\α}-1},$ i＝1,2,......,f-1               (1)

$X_i\left(k\right)=\frac{P_s\left(k\right)\mathrm{\α}{e}^{\frac{T_i\left(k\right)+c}{b}-a}-1}{\mathrm{\α}-1},$ i＝f,f+1,......,n               (2)

式中，k为当前采样时刻，P_r(k)为k采样时刻的精馏段压强、P_s(k)为k采样时刻的提馏段压强，T_i(k)为k采样时刻第i块塔板的温度，i表示塔板编号(i＝1,2,...,f,f+1,...,n，其中，1为塔顶编号，f为进料板编号，n为塔底编号)，P_r(k)、P_s(k)及T_i(k)由智能仪表测得，α为相对挥发度，a、b、c为安东尼常数，X_i(k)为k采样时刻第i块板塔的液相轻组分的浓度测量值；

2)各塔板的浓度梯度观测，根据公式(3)、(4)、(5)、(6)获得各塔板当前时刻的浓度梯度值，并将结果传输至温度梯度描述模块、数据存储装置及上位机显示界面，公式(3)、(4)、(5)、(6)的表达式如下：

$\frac{d{X}_1\left(k\right)}{\mathrm{dt}}=\frac{1}{H}[V_2\left(k\right)Y_2\left(k\right)-V_1\left(k\right)Y_1\left(k\right)-L_1\left(k\right)X_1\left(k\right)]---\left(3\right)$

$\frac{d{X}_i\left(k\right)}{\mathrm{dt}}=\frac{1}{H}[V_{i+1}\left(k\right)Y_{i+1}\left(k\right)-V_i\left(k\right)Y_i\left(k\right)+L_{i-1}\left(k\right)X_{i-1}\left(k\right)-L_i\left(k\right)X_i\left(k\right)]---\left(4\right)$

                         (i＝2,...,n-1并且i≠f)

$\frac{d{X}_f\left(k\right)}{\mathrm{dt}}=\frac{1}{H}[V_{f+1}\left(k\right)Y_{f+1}\left(k\right)-V_f\left(k\right)Y_f\left(k\right)+L_{f-1}\left(k\right)X_{f-1}\left(k\right)-L_f\left(k\right)X_f\left(k\right)+F\left(k\right)Z_f\left(k\right)]---\left(5\right)$

$\frac{d{X}_n\left(k\right)}{\mathrm{dt}}=\frac{1}{H}[-V_n\left(k\right)Y_n\left(k\right)+L_{n-1}\left(k\right)X_{n-1}\left(k\right)-L_n\left(k\right)X_n\left(k\right)]---\left(6\right)$

式中，H为持液量，V_i(k)为k采样时刻第i块板塔的气相流率，L_i(k)为k采样时刻第i块板塔的液相流率，X_i(k)为k采样时刻第i块板塔的液相轻组分的浓度测量值，Y_i(k)为k采样时刻第i块板塔的气相轻组分浓度，为k采样时刻第i块板塔的液相轻组分的浓度梯度值，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号，n为塔底编号)，F(k)为k采样时刻的进料流量，Z_f(k)为k采样时刻的进料组分；所述Y_i(k)由公式(7)得到，公式(7)的表达式如下：

Y_i(k)＝αX_i(k)/[(α-1)X_i(k)+1]  i＝1,2,...,f,f+1,...,n          (7)

式中，α为相对挥发度，X_i(k)为k采样时刻第i块板塔的液相轻组分的浓度测量值；

(二)所述温度梯度描述模块观测温度梯度的步骤为：通过数据存储装置提取智能仪表中的温度检测元件、压力检测元件收集的温度、压力参数，以及所述浓度梯度描述模块得到的浓度梯度信息，通过所述温度梯度描述模块确定温度梯度，所述温度梯度描述模块根据公式(8)、(9)获得温度梯度，公式(8)、(9)的表达式如下：

$\frac{d{T}_i\left(k\right)}{\mathrm{dt}}=-\frac{d{X}_i\left(k\right)}{\mathrm{dt}}\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_r\left(k\right){\mathrm{\αbe}}^{\frac{b}{T_i\left(k\right)+c}-a}},$ i＝1,2,......,f-1              (8)

$\frac{d{T}_i\left(k\right)}{\mathrm{dt}}=-\frac{d{X}_i\left(k\right)}{\mathrm{dt}}\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_s\left(k\right){\mathrm{\αbe}}^{\frac{b}{T_i\left(k\right)+c}-a}},$ i＝f,f+1,......,n             (9)

式中，为k采样时刻第i块板塔的液相轻组分的浓度梯度值，α为相对挥发度，a、b、c为安东尼常数，T_i(k)为k采样时刻第i块塔板的温度，P_r(k)为k采样时刻的精馏段压强、P_s(k)为k采样时刻的提馏段压强，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号，n为塔底编号)，为k采样时刻第i块板塔的温度梯度值；

(三)所述温度波静态描述模块观测温度波静态步骤为：根据内部热耦合精馏塔的温度波特性，通过温度检测元件采集各塔板当前时刻的温度测量值，并将各温度值在坐标轴中对应的点连成连续的光滑曲线，获得当前时刻精馏段与提馏段的温度波形，进而得到温度波静态描述函数的常系数T_r1、T_r2、T_s1、T_s2、γ_r、γ_s的值，以及精馏段、提馏段温度波的拐点初始值S_r(0)、S_s(0)，通过温度波静态描述函数公式(10)、(11)获得各塔板的温度预测初始值公式(10)、(11)的表达式如下：

${\hat{T}}_i\left(k\right)=\frac{T_{r1}{e}^{-{\mathrm{\γ}}_r(i-S_r\left(k\right))}+T_{r2}}{1+e^{-{\mathrm{\γ}}_r(i-S_r\left(k\right))}},$ i＝1,2,...,f-1              (10)

${\hat{T}}_i\left(k\right)=\frac{T_{s1}{e}^{-{\mathrm{\γ}}_s(i-S_s\left(k\right))}+T_{s2}}{1+e^{-{\mathrm{\γ}}_s(i-S_s\left(k\right))}},$ i＝f,f+1,......,n             (11)

式中，k为当前采样时刻，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号,n为塔底编号)，为k采样时刻第i块板塔的温度预测值，S_r(k)、S_s(k)分别为k采样时刻内部热耦合精馏塔精馏段、提馏段温度波的拐点，T_r1、T_r2、T_s1、T_s2、γ_r、γ_s为温度波静态描述函数的常系数，T_r1、T_r2、T_s1、T_s2分别表示精馏段与提馏段温度波两端的渐进浓度，γ_r、γ_s分别表征精馏段与提馏段温度波拐点处的斜率大小；

(四)所述温度波动态描述模块观测温度波动态的步骤为：根据所述温度梯度描述模块和所述温度波静态描述模块所获得的信息，确定各塔板温度波在未来时刻的变化趋势，该变化趋势用拐点的移动速度来表示，具体由公式(12)、(13)获得，公式(12)、(13)的表达式如下：

$\frac{d{S}_r}{\mathrm{dt}}\left(k\right)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}\left[\frac{d{T}_i\left(k\right)}{\mathrm{dt}}\frac{T_{r2}-T_{r1}}{{\mathrm{\γ}}_r(T_i\left(k\right)-T_{r2})(T_i\left(k\right)-T_{r1})}\right]---\left(12\right)$

$\frac{d{S}_s}{\mathrm{dt}}\left(k\right)=\underset{i=f}{\overset{n}{\mathrm{\Σ}}}\left[\frac{d{T}_i\left(k\right)}{\mathrm{dt}}\frac{T_{s2}-T_{s1}}{{\mathrm{\γ}}_s(T_i\left(k\right)-T_{s2})(T_i\left(k\right)-T_{s1})}\right]---\left(13\right)$

式中，分别为k采样时刻精馏段和提馏段温度波的拐点移动速度，为k采样时刻第i块板塔的温度梯度值，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号,n为塔底编号)，T_i(k)为k采样时刻第i块塔板的温度，T_r1、T_r2、T_s1、T_s2、γ_r、γ_s为温度波静态描述函数的常系数；

(五)所述误差计算模块用观测计算误差的步骤为：根据公式(14)计算k采用时刻的温度预测误差，公式(14)的表达式如下：

$e_i\left(k\right)=T_i\left(k\right)-{\hat{T}}_i\left(k\right)---\left(14\right)$

式中，e_i(k)为k采样时刻的温度预测误差，T_i(k)为k采样时刻第i块塔板的温度测量值，为k采样时刻第i块塔板的温度预测值；

(六)所述未来时刻状态模块其观测未来时刻状态的步骤为：根据当前时刻的参数信息预测未来时刻的状态，包括以下三部分：

1)根据公式(15)、(16)、(17)、(18)、(19)(20)获得下一时刻拐点位置以及温度和组分浓度，公式(15)、(16)、(17)、(18)、(19)(20)的表达式如下：

$S_r(k+1)=S_r\left(k\right)+\frac{d{S}_r}{\mathrm{dt}}\left(k\right)\×T---\left(15\right)$

$S_s(k+1)=S_s\left(k\right)+\frac{d{S}_s}{\mathrm{dt}}\left(k\right)\×T---\left(16\right)$

${\hat{T}}_i(k+1)=\frac{T_{r1}{e}^{-{\mathrm{\γ}}_r(i-S_r(k+1))}+T_{r2}}{1+e^{-{\mathrm{\γ}}_r(i-S_r(k+1))}},$ i＝1,2,...,f-1              (17)

${\hat{T}}_i(k+1)=\frac{T_{s1}{e}^{-{\mathrm{\γ}}_s(i-S_s(k+1))}+T_{s2}}{1+e^{-{\mathrm{\γ}}_s(i-S_s(k+1))}},$ i＝f,f+1,......,n            (18)

${\hat{X}}_i(k+1)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_r\left(k\right)e^{\frac{{\hat{T}}_i(k+1)+c}{b}-a}-\frac{1}{\mathrm{\α}}),$ i＝1,2,......,f-1           (19)

${\hat{X}}_i(k+1)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_s\left(k\right)e^{\frac{{\hat{T}}_i(k+1)+c}{b}-a}-\frac{1}{\mathrm{\α}}),$ i＝f,f+1,......,n           (20)

式中，T为采样周期，i表示塔板编号(i＝1,2,...,f,f+1,...,n，其中，1为塔顶编号，f为进料板编号，n为塔底编号)，α为相对挥发度，a、b、c为安东尼常数，P_r(k)为k采样时刻的精馏段压强、P_s(k)为k采样时刻的提馏段压强，S_r(k+1)、S_s(k+1)分别为k+1采样时刻内部热耦合精馏塔精馏段、提馏段温度波的拐点，为k+1采样时刻第i块塔板的温度预测值，为k+1采样时刻第i块塔板的液相轻组分浓度的预测值；

2)修正温度与组分浓度预测值，根据所述误差计算模块得到的温度预测误差，由公式(21)、(22)、(23)获得修正后的温度与组分浓度的预测值，公式(21)、(22)、(23)的表达式如下：

${\tilde{T}}_i(k+1)={\hat{T}}_i(k+1)+e_i\left(k\right)---\left(21\right)$

${\tilde{X}}_i(k+1)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_r\left(k\right)e^{\frac{{\tilde{T}}_i(k+1)+c}{b}-a}-\frac{1}{\mathrm{\α}}),$ i＝1,2,......,f-1            (22)

${\tilde{X}}_i(k+1)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_s\left(k\right)e^{\frac{{\tilde{T}}_i(k+1)+c}{b}-a}-\frac{1}{\mathrm{\α}}),$ i＝f,f+1,......,n            (23)

式中，i表示塔板编号(i＝1,2,...，f,f+1,...,n，其中，1为塔顶编号，f为进料板编号，n为塔底编号)，α为相对挥发度，a、b、c为安东尼常数，P_r(k)为k采样时刻的精馏段压强、P_s(k)为k采样时刻的提馏段压强，表示修正后的k+1采样时刻第i块塔板的温度预测值，表示修正后的k+1采样时刻第i块塔板的液相轻组分浓度的预测值；

3)预测未来m个时刻状态预测，即对k+2到k+m(m>2)采样时刻的温度和组分浓度进行预测，通过以下步骤完成：

A、对k+2采样时刻的温度和组分进行预测，首先由公式(24)、(25)得到k+1采样时刻的拐点的移动速度，公式(24)、(25)的表达式如下：

$\frac{d{S}_r}{\mathrm{dt}}(k+1)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}\left[\frac{{\tilde{T}}_i(k+1)-T_i\left(k\right)}{T}\frac{T_{r2}-T_{r1}}{{\mathrm{\γ}}_r({\tilde{T}}_i(k+1)-T_{r2})({\tilde{T}}_i(k+1)-T_{r1})}\right]---\left(24\right)$

$\frac{d{S}_s}{\mathrm{dt}}(k+1)=\underset{i=f}{\overset{n}{\mathrm{\Σ}}}\left[\frac{{\tilde{T}}_i(k+1)-T_i\left(k\right)}{T}\frac{T_{s2}-T_{s1}}{{\mathrm{\γ}}_s({\tilde{T}}_i(k+1)-T_{s2})({\tilde{T}}_i(k+1)-T_{s1})}\right]---\left(25\right)$

进而由公式(26)、(27)得到k+2采样时刻的拐点位置，公式(26)、(27)的表达式如下：

$S_r(k+2)=S_r(k+1)+\frac{d{S}_r}{\mathrm{dt}}(k+1)\×T---\left(26\right)$

$S_s(k+2)=S_s(k+1)+\frac{d{S}_s}{\mathrm{dt}}(k+1)\×T---\left(27\right)$

由于所述内部热耦合精馏塔在一个生产周期内提馏段压强变化很小，则认定为P_s(k+m-1)＝...＝P_s(k+1)＝P_s(k)，精馏段压强是主要的操作变量，由所述控制站的控制策略决定，所述控制站将未来时刻的操作变量值传输至所述数据存储装置，再由所述数据存储装置传输给所述上位机观测器，进而得到未来m-1个时刻的精馏段压强P_r(k+1),P_r(k+2),...,P_r(k+m-1)；根据公式(28)、(29)、(30)、(31)得到k+2采用时刻的温度和组分浓度预测，公式(28)、(29)、(30)、(31)的表达式如下：

${\hat{T}}_i(k+2)=\frac{T_{r1}{e}^{-{\mathrm{\γ}}_r(i-S_r(k+2))}+T_{r2}}{1+e^{-{\mathrm{\γ}}_r(i-S_r(k+2))}},$ i＝1,2,...,f-1             (28)

${\hat{T}}_i(k+2)=\frac{T_{s1}{e}^{-{\mathrm{\γ}}_s(i-S_s(k+2))}+T_{s2}}{1+e^{-{\mathrm{\γ}}_s(i-S_s(k+2))}},$ i＝f,f+1,......,n            (29)

${\hat{X}}_i(k+2)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_r(k+1)e^{\frac{{\hat{T}}_i(k+2)+c}{b}-a}-\frac{1}{\mathrm{\α}}),$ i＝1,2,......,f-1           (30)

${\hat{X}}_i(k+2)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_s(k+1)e^{\frac{{\hat{T}}_i(k+2)+c}{b}-a}-\frac{1}{\mathrm{\α}}),$ i＝f,f+1,......,n           (31)

设定误差不发生变化，则由公式(32)、(33)、(34)获得修正后的温度与组分浓度的预测值，公式(32)、(33)、(34)的表达式如下：

${\tilde{T}}_i(k+2)={\hat{T}}_i(k+2)+e_i\left(k\right)---\left(32\right)$

${\tilde{X}}_i(k+2)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_r(k+1)e^{\frac{{\tilde{T}}_i(k+2)+c}{b}-a}-\frac{1}{\mathrm{\α}}),$ i＝1,2,......,f-1        (33)

${\tilde{X}}_i(k+2)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_s(k+1)e^{\frac{{\tilde{T}}_i(k+2)+c}{b}-a}-\frac{1}{\mathrm{\α}}),$ i＝f,f+1,......,n           (34)

B、通过步骤A中所述的方法，进而可以递推求出k+3到k+m采样时刻的温度和组分浓度预测值，由公式(35)、(36)、(37)获得第k+m采样时刻的温度与组分浓度的预测值，公式(35)、(36)、(37)的表达式如下：

${\tilde{T}}_i(k+m)={\hat{T}}_i(k+m)+e_i\left(k\right)---\left(35\right)$

${\tilde{X}}_i(k+m)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_r(k+m-1)e^{\frac{{\tilde{T}}_i(k+m)+c}{b}-a}-\frac{1}{\mathrm{\α}}),$ i＝1,2,......,f-1       (36)

${\tilde{X}}_i(k+m)=\frac{\mathrm{\α}}{\mathrm{\α}-1}(P_s\left(k\right)e^{\frac{{\tilde{T}}_i(k+m)+c}{b}-a}-\frac{1}{\mathrm{\α}}),$ i＝f,f+1,......,n           (37)

式中，i表示塔板编号(i＝1,2,...,f,f+1,...,n，其中，1为塔顶编号，f为进料板编号，n为塔底编号)，α为相对挥发度，a、b、c为安东尼常数，表示修正后的k+m采样时刻第i块塔板的温度预测值，表示修正后的k+m采样时刻第i块塔板的液相轻组分浓度的预测值。

2.根据权利要求1所述的基于温度波特性的内部热耦合精馏在线观测器，其特征在于：所述上位机观测器还用于设定采用周期T以及预测时域m，并将上位机观测器输出的当前时刻温度、浓度测量值，以及未来m个时刻的温度、浓度预测值，传输到上位机显示界面，并将以上值传输给数据存储装置。