CN104587695A - 基于温度波特性的内部热耦合精馏塔控制装置 - Google Patents

Abstract

本发明涉及一种基于温度波特性的内部热耦合精馏塔控制装置，包括内部热耦合精馏塔、智能仪表、控制站、数据存储装置及上位机，智能仪表与内部热耦合精馏塔连接，控制站与内部热耦合精馏塔连接，数据存储装置与智能仪表和控制站连接，上位机与数据存储装置和控制站连接，上位机包括浓度梯度描述模块、温度梯度描述模块、温度波静态描述模块、温度波动态描述模块、设定值转换模块及控制参数求解模块。本发明非线性控制方案建立在高精度非线性模型基础上，能够及时抑制干扰作用；较好地处理了耦合问题，能够快速准确地跟踪设定值变化；以测量精度高、测量延时小的温度为主要控制手段，大大提高了控制速度与控制性能。

Description

基于温度波特性的内部热耦合精馏塔控制装置

技术领域

本发明属于精馏节能过程中的非线性控制领域，具体地说，涉及一种基于温度波特性的内部热耦合精馏塔控制装置。

背景技术

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

内部热耦合精馏技术充分利用精馏段与提馏段之间的热交换，比常规精馏节能30％以上。然而，内部热耦合精馏过程的热耦合导致该过程具有显著的非线性动态特性，使得该塔的控制策略设计显得尤为困难。传统的PID控制方案等已经不能满足要求，在内部热耦合精馏塔的过程控制当中，这些方案已经很难使精馏过程稳定。而基于线性辨识模型的控制方案只能工作在稳态工作点附近，稍微增大干扰幅度，或者设定值阶跃变化，系统控制质量则出现明显下降。因此，基于内部热耦合精馏塔的非线性特性，并在此基础上实现内部热耦合精馏塔高效节能过程有效的的非线性控制方案，是提高内部热耦合精馏过程的产品品质的保障，已经成为一项关键的精馏节能技术，具有十分重要的意义。

发明内容

本发明针对现有内部热耦合精馏的控制装置存在的在线运行效率低下、抑制干扰能力差、控制效果差、对噪声敏感度低等上述不足，提供了一种基于温度波特性的内部热耦合精馏塔控制装置，该控制装置能够实现精确迅速的设定值跟踪，具有在线运行速度快，抗噪声能力强、控制效果好等优点。

本发明的技术方案是：一种基于温度波特性的内部热耦合精馏塔控制装置，包括内部热耦合精馏塔、智能仪表、控制站、数据存储装置及上位机，所述智能仪表与所述内部热耦合精馏塔连接，用于进行数据采集；所述控制站与所述内部热耦合精馏塔连接，用于实现对内部热耦合精馏塔的控制；所述数据存储装置与所述智能仪表和所述控制站连接，用于实现数据存储；所述上位机与所述数据存储装置和所述控制站连接，用于实现控制参数的求解，所述上位机包括用于观测浓度及浓度梯度的浓度梯度描述模块、用于观测温度梯度的温度梯度 描述模块、用于观测温度波静态的温度波静态描述模块、用于观测温度波动态的温度波动态描述模块、用于实现设定值转换的设定值转换模块及用于求解控制参数的控制参数求解模块；其中，

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

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{{\mathrm{dX}}_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{{\mathrm{dX}}_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{{\mathrm{dX}}_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{{\mathrm{dX}}_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)

式中，α为相对挥发度，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号，n为塔底编号)，X_i(k)为k采样时刻第i块板塔的液相轻组分的浓度测量值；

所述气液相流率由公式(8)、(9)、(10)、(11)、(12)、(13)、(14)得到，公式(8)、(9)、(10)、(11)、(12)、(13)、(14)的表达式如下：

V₁(k)＝F(k)(1-q(k+1))      (8) 

L_n(k)＝F(k)q(k+1)        (9) 

$L_i\left(k\right)=\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{Q}_j\left(k\right)/\mathrm{\λ},(i=1,...,f-1)---\left(10\right)$

V_i+1(k)＝V₁(k)+L_i(k)(i＝1,...,f-1)      (11)

$L_{f+i-1}\left(k\right)=L_{f-1}\left(k\right)+F\left(k\right)q(k+1)-\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{Q}_j\left(k\right)/\mathrm{\λ},(i=1,...,f-2)---\left(12\right)$

$L_{f+i}\left(k\right)=V_f\left(k\right)-F\left(k\right)(1-q(k+1))-\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{Q}_j\left(k\right)/\mathrm{\λ},(i=1,...,f-2)---\left(13\right)$

Q_i(k)＝UA×(T_i(k)-T_i+f-1(k)),i＝1,...,f-1       (14) 

式中，i表示塔板编号，f为进料板编号，Q_i(k)为第i块塔板上的热耦合量，UA为传热速率，λ为汽化潜热，q(k+1)为第k+1采样时刻的进料热状况，其中，q(k+1)为所述控制装置下一时刻的两个控制参数之一；

通过公式(3)-(14)，建立浓度梯度与进料热状况之间的关系，浓度梯度与进料热状况之间的关系由公式(15)简化表示，公式(15)的表达式如下：

$\frac{{\mathrm{dX}}_i\left(k\right)}{\mathrm{dt}}=f_i\left(q(k+1)\right),(i=\mathrm{1,2},...,n)---\left(15\right)$

式中，i表示塔板编号(1为塔顶编号，n为塔底编号)，f_i表示由公式(3)-(14)得到的结构已 知的非线性函数关系；

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

$\frac{{\mathrm{dT}}_i\left(k\right)}{\mathrm{dt}}=-\frac{{\mathrm{dX}}_i\left(k\right)}{\mathrm{dt}}\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_r(k+1)\mathrm{\α}{\mathrm{be}}^{\frac{b}{T_i\left(k\right)+c}-a}},i=\mathrm{1,2},......,f-1---\left(16\right)$

$\frac{{\mathrm{dT}}_i\left(k\right)}{\mathrm{dt}}=-\frac{{\mathrm{dX}}_i\left(k\right)}{\mathrm{dt}}\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_s\left(k\right)\mathrm{\α}{\mathrm{be}}^{\frac{b}{T_i\left(k\right)+c}-a}},i=f,f+1,......,n---\left(17\right)$

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

将公式(15)代入公式(16)、(17)，进一步得到公式(18)、(19)，公式(18)、(19)的表达式如下：

$\frac{{\mathrm{dT}}_i\left(k\right)}{\mathrm{dt}}=-f_i\left(q(k+1)\right)\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_r(k+1)\mathrm{\α}{\mathrm{be}}^{\frac{b}{T_i\left(k\right)+c}-a}},i=\mathrm{1,2},......,f-1---\left(18\right)$

$\frac{{\mathrm{dT}}_i\left(k\right)}{\mathrm{dt}}=-f_i\left(q(k+1)\right)\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_s\left(k\right)\mathrm{\α}{\mathrm{be}}^{\frac{b}{T_i\left(k\right)+c}-a}},i=f,f+1,......,n---\left(19\right)$

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

${\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(20\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(21\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分别表征精馏段与提馏段温度波拐点处的斜率大小；

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

$\frac{{\mathrm{dS}}_r}{\mathrm{dt}}\left(k\right)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}\left[\frac{{\mathrm{dT}}_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(22\right)$

$\frac{{\mathrm{dS}}_s}{\mathrm{dt}}\left(k\right)=\underset{i=f}{\overset{n}{\mathrm{\Σ}}}\left[\frac{{\mathrm{dT}}_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(23\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为温度波静态描述函数的常系数；

将公式(18)、(19)代入公式(22)、(23)，进一步得到公式(24)、(25)，公式(24)、(25)的表达式如下：

$\frac{{\mathrm{dS}}_r}{\mathrm{dt}}\left(k\right)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}[-f_i\left(q(k+1)\right)\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_r(k+1)\mathrm{\α}{\mathrm{be}}^{\frac{b}{T_i\left(k\right)+c}-a}}\frac{T_{r2}-T_{r1}}{{\mathrm{\γ}}_r(T_i\left(k\right)-T_{r2})(T_i\left(k\right)-T_{r1})}]---\left(24\right)$

$\frac{{\mathrm{dS}}_s}{\mathrm{dt}}\left(k\right)=\underset{i=f}{\overset{n}{\mathrm{\Σ}}}[-f_i\left(q(k+1)\right)\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_s\left(k\right)\mathrm{\α}{\mathrm{be}}^{\frac{b}{T_i\left(k\right)+c}-a}}\frac{T_{s2}-T_{s1}}{{\mathrm{\γ}}_s(T_i\left(k\right)-T_{s2})(T_i\left(k\right)-T_{s1})}]---\left(25\right)$

(五)所述设定值转换模块实现设定值转换的步骤为：根据温度波特静态描述函数，通过转换公式(26)、(27)将浓度设定值转换为温度波的拐点设定值，转换公式(26)、(27)的表达式如下：

$\frac{b}{a-\ln [P_r\left(k\right)\×(\mathrm{\α}-(\mathrm{\α}-1){Y_1}^{*})]}-c=\frac{T_{r1}{e}^{-{\mathrm{\γ}}_r(1-{S_r}^{*})}+T_{r2}}{1+e^{-{\mathrm{\γ}}_r(1-{S_r}^{*})}}---\left(26\right)$

$\frac{b}{a-\ln \frac{p_s\left(k\right)}{{X_n}^{*}+(1-{X_n}^{*})/\mathrm{\α}}}-c=\frac{T_{s1}{e}^{-{\mathrm{\γ}}_s(n-{S_s}^{*})}+T_{s2}}{1+e^{-{\mathrm{\γ}}_s(n-{S_s}^{*})}}---\left(27\right)$

式中，Y₁ ^*、X_n ^*分别为塔顶的气相轻组分浓度的设定值和塔底的液相轻组分浓度的设定值，S_r ^*、S_s ^*分别为精馏段和提馏段温度波拐点的设定值，P_r(k)为k采样时刻精馏段压强、P_s(k)为k采样时刻提馏段压强，α为相对挥发度，a、b、c为安东尼常数；

(六)控制参数求解模块求解控制参数的步骤为：根据所述各个描述模块和所述设定值转换模块得到的信息求取下一时刻的控制参数，其中，控制率由公式(28)、(29)描述，公式(28)、(29)的表达式如下：

$\frac{{\mathrm{dS}}_r}{\mathrm{dt}}\left(k\right)=K_1({S_r}^{*}-S_r\left(k\right))+K_{2}T\underset{i=1}{\overset{k}{\mathrm{\Σ}}}({S_r}^{*}-S_r\left(i\right))---\left(28\right)$

$\frac{{\mathrm{dS}}_s}{\mathrm{dt}}\left(k\right)=K_3({S_s}^{*}-S_s\left(k\right))+K_{4}T\underset{i=1}{\overset{k}{\mathrm{\Σ}}}({S_s}^{*}-S_s\left(i\right))---\left(29\right)$

式中，T为采样周期，S_r(k)、S_r(i)分别为k和i采样时刻内部热耦合精馏塔精馏段温度波的拐点位置，S_s(k)、S_s(i)分别为k和i采样时刻内部热耦合精馏塔提馏段温度波的拐点位置，S_r ^*、S_s ^*分别为精馏段和提馏段温度波拐点的设定值，K₁、K₂、K₃、K₄为控制器参数参数，可根据实际的控制需求调整，其中，K₁和K₃在20-80之间取值，K₂和K₄在50-800之间取值；将公式(24)、(25)代入到公式(28)、(29)，即可求解出所述控制装置下一时刻的控制参数，即第k+1采样时刻的进料热状况q(k+1)和精馏段压强P_r(k+1)。

作为优选，所述上位机还用于设定采用周期T，设定K₁、K₂、K₃、K₄四个系统参数的值，以及设定塔顶的气相轻组分浓度和塔底的液相轻组分浓度的设定值Y₁ ^*、X_n ^*，显示精馏段和提馏段温度波拐点的设定值S_r ^*、S_s ^*及当前时刻温度波拐点位置和控制参数求解模块求解出的下一时刻的控制参数，并将控制参数传递给所述控制站，所述控制站根据得到的控制参数对控制器进行调整，进而实现对内部热耦合精馏塔的控制调整；所述上位机还将以上信息传递给数据存储装置。

本发明的有益效果是：本发明对内部热耦合精馏过程中的温度波特性进行精确描述，成功准确地把握内部热耦合精馏塔的非线性动态特性，克服已有的控制装置在线运行效率底下、抑制干扰能力差、控制效果差、对噪声敏感度低的不足，具有运行速度快、抗噪声能力强、控制效果好的优点，可以实现精确迅速的设定值跟踪。本发明非线性控制方案建立在高精度 非线性模型基础上，能够及时抑制干扰作用；较好地处理了耦合问题，能够快速准确地跟踪设定值变化；以测量精度高、测量延时小的温度为主要控制手段，大大提高了控制速度与控制性能。

附图说明

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

附图2为本发明具体实施例上位机实现观测及求解控制参数的原理图。

具体实施方式

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

如图1和图2所示，一种基于温度波特性的内部热耦合精馏塔控制装置，包括内部热耦合精馏塔1、智能仪表2、控制站3、数据存储装置4及上位机5，所述智能仪表2与所述内部热耦合精馏塔1连接，用于进行数据采集；所述控制站3与所述内部热耦合精馏塔1连接，用于实现对内部热耦合精馏塔1的控制；所述数据存储装置4与所述智能仪表2和所述控制站3连接，用于实现数据存储；所述上位机5与所述数据存储装置4和所述控制站3连接，用于实现控制参数的求解，所述上位机5包括用于观测浓度及浓度梯度的浓度梯度描述模块6、用于观测温度梯度的温度梯度描述模块7、用于观测温度波静态的温度波静态描述模块8、用于观测温度波动态的温度波动态描述模块9、用于实现设定值转换的设定值转换模块10及用于求解控制参数的控制参数求解模块11；其中，

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

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{{\mathrm{dX}}_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{{\mathrm{dX}}_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)-K_i\left(k\right)X_i\left(k\right)]---\left(4\right)$

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

$\frac{{\mathrm{dX}}_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{{\mathrm{dX}}_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) 

式中，α为相对挥发度，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号，n为塔底编号)，X_i(k)为k采样时刻第i块板塔的液相轻组分的浓度测量值；

所述气液相流率由公式(8)、(9)、(10)、(11)、(12)、(13)、(14)得到，公式(8)、(9)、(10)、(11)、(12)、(13)、(14)的表达式如下：

V₁(k)＝F(k)(1-q(k+1))       (8) 

L_n(k)＝F(k)q(k+1)          (9) 

$L_i\left(k\right)=\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{Q}_j\left(k\right)/\mathrm{\λ},(i=1,...,f-1)---\left(10\right)$

V_i+1(k)＝V₁(k)+L_i(k)  (i＝1,...,f-1)       (11)

$L_{f+i-1}\left(k\right)=L_{f-1}\left(k\right)+F\left(k\right)q(k+1)-\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{Q}_j\left(k\right)/\mathrm{\λ},(i=1,...,f-2)---\left(12\right)$

$L_{f+i}\left(k\right)=V_f\left(k\right)-F\left(k\right)(1-q(k+1))-\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{Q}_j\left(k\right)/\mathrm{\λ},(i=1,...,f-2)---\left(13\right)$

Q_i(k)＝UA×(T_i(k)-T_i+f-1(k)),i＝1,...,f-1      (14) 

式中，i表示塔板编号，f为进料板编号，Q_i(k)为第i块塔板上的热耦合量，UA为传热速率，λ为汽化潜热，q(k+1)为第k+1采样时刻的进料热状况，其中，q(k+1)为所述控制装置下一时刻的两个控制参数之一；

通过公式(3)-(14)，建立浓度梯度与进料热状况之间的关系，浓度梯度与进料热状况之间的关系由公式(15)简化表示，公式(15)的表达式如下：

$\frac{{\mathrm{dX}}_i\left(k\right)}{\mathrm{dt}}=f_i\left(q(k+1)\right),(i=\mathrm{1,2},...,n)---\left(15\right)$

式中，i表示塔板编号(1为塔顶编号，n为塔底编号)，f_i表示由公式(3)-(14)得到的结构已知的非线性函数关系；

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

$\frac{{\mathrm{dT}}_i\left(k\right)}{\mathrm{dt}}=-\frac{{\mathrm{dX}}_i\left(k\right)}{\mathrm{dt}}\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_r(k+1)\mathrm{\α}{\mathrm{be}}^{\frac{b}{T_i\left(k\right)+c}-a}},i=\mathrm{1,2},......,f-1---\left(16\right)$

$\frac{{\mathrm{dT}}_i\left(k\right)}{\mathrm{dt}}=-\frac{{\mathrm{dX}}_i\left(k\right)}{\mathrm{dt}}\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_s\left(k\right)\mathrm{\α}{\mathrm{be}}^{\frac{b}{T_i\left(k\right)+c}-a}},i=f,f+1,......,n---\left(17\right)$

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

将公式(15)代入公式(16)、(17)，进一步得到公式(18)、(19)，公式(18)、(19)的表达式如下：

$\frac{{\mathrm{dT}}_i\left(k\right)}{\mathrm{dt}}=-f_i\left(q(k+1)\right)\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_r(k+1)\mathrm{\α}{\mathrm{be}}^{\frac{b}{T_i\left(k\right)+c}-a}},i=\mathrm{1,2},......,f-1---\left(18\right)$

$\frac{{\mathrm{dT}}_i\left(k\right)}{\mathrm{dt}}=-f_i\left(q(k+1)\right)\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_s\left(k\right)\mathrm{\α}{\mathrm{be}}^{\frac{b}{T_i\left(k\right)+c}-a}},i=f,f+1,......,n---\left(19\right)$

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

${\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(20\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(21\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分别表征精馏段与提馏段温度波拐点处的斜率大小；

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

$\frac{{\mathrm{dS}}_r\left(k\right)}{\mathrm{dt}}=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}\left[\frac{{\mathrm{dT}}_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(22\right)$

$\frac{{\mathrm{dS}}_s\left(k\right)}{\mathrm{dt}}=\underset{i=f}{\overset{n}{\mathrm{\Σ}}}\left[\frac{{\mathrm{dT}}_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(23\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为温度波静态描述函数的常系数；

将公式(18)、(19)代入公式(22)、(23)，进一步得到公式(24)、(25)，公式(24)、(25)的表达式如下：

$\frac{{\mathrm{dS}}_r}{\mathrm{dt}}\left(k\right)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}[-f_i\left(q(k+1)\right)\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_r(k+1)\mathrm{\α}{\mathrm{be}}^{\frac{b}{T_i\left(k\right)+c}-a}}\frac{T_{r2}-T_{r1}}{{\mathrm{\γ}}_r(T_i\left(k\right)-T_{r2})(T_i\left(k\right)-T_{r1})}]---\left(24\right)$

$\frac{{\mathrm{dS}}_s}{\mathrm{dt}}\left(k\right)=\underset{i=f}{\overset{n}{\mathrm{\Σ}}}[-f_i\left(q(k+1)\right)\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_s\left(k\right)\mathrm{\α}{\mathrm{be}}^{\frac{b}{T_i\left(k\right)+c}-a}}\frac{T_{s2}-T_{s1}}{{\mathrm{\γ}}_s(T_i\left(k\right)-T_{s2})(T_i\left(k\right)-T_{s1})}]---\left(25\right)$

(五)所述设定值转换模块10实现设定值转换的步骤为：根据温度波特静态描述函数，通过转换公式(26)、(27)将浓度设定值转换为温度波的拐点设定值，转换公式(26)、(27)的表达式如下：

$\frac{b}{a-\ln [P_r\left(k\right)\×(\mathrm{\α}-(\mathrm{\α}-1){Y_1}^{*})]}-c=\frac{T_{r1}{e}^{-{\mathrm{\γ}}_r(1-{S_r}^{*})}+T_{r2}}{1+e^{-{\mathrm{\γ}}_r(1-{S_r}^{*})}}---\left(26\right)$

$\frac{b}{a-\ln \frac{p_s\left(k\right)}{{X_n}^{*}+(1-{X_n}^{*})/\mathrm{\α}}}-c=\frac{T_{s1}{e}^{-{\mathrm{\γ}}_s(n-{S_s}^{*})}+T_{s2}}{1+e^{-{\mathrm{\γ}}_s(n-{S_s}^{*})}}---\left(27\right)$

式中，Y₁ ^*、X_n ^*分别为塔顶的气相轻组分浓度的设定值和塔底的液相轻组分浓度的设定值，S_r ^*、S_s ^*分别为精馏段和提馏段温度波拐点的设定值，P_r(k)为k采样时刻精馏段压强、P_s(k)为k采样时刻提馏段压强，α为相对挥发度，a、b、c为安东尼常数；

(六)控制参数求解模块11求解控制参数的步骤为：根据所述各个描述模块和所述设定值转换模块10得到的信息求取下一时刻的控制参数，其中，控制率由公式(28)、(29)描述，公式(28)、(29)的表达式如下：

$\frac{{\mathrm{dS}}_r}{\mathrm{dt}}\left(k\right)=K_1({S_r}^{*}-S_r\left(k\right))+K_{2}T\underset{i=1}{\overset{k}{\mathrm{\Σ}}}({S_r}^{*}-S_r\left(i\right))---\left(28\right)$

$\frac{{\mathrm{dS}}_s}{\mathrm{dt}}\left(k\right)=K_3({S_s}^{*}-S_s\left(k\right))+K_{4}T\underset{i=1}{\overset{k}{\mathrm{\Σ}}}({S_s}^{*}-S_s\left(i\right))---\left(29\right)$

式中，T为采样周期，S_r(k)、S_r(i)分别为k和i采样时刻内部热耦合精馏塔精馏段温度波的拐点位置，S_s(k)、S_s(i)分别为k和i采样时刻内部热耦合精馏塔提馏段温度波的拐点位置， S_r ^*、S_s ^*分别为精馏段和提馏段温度波拐点的设定值，K₁、K₂、K₃、K₄为控制器参数参数，可根据实际的控制需求调整，其中，K₁和K₃在20-80之间取值，K₂和K₄在50-800之间取值；

将公式(24)、(25)代入到公式(28)、(29)，即可求解出所述控制装置下一时刻的控制参数，即第k+1采样时刻的进料热状况q(k+1)和精馏段压强P_r(k+1)。

本实施例中，所述上位机5还用于设定采用周期T，设定K₁、K₂、K₃、K₄四个系统参数的值，以及设定塔顶的气相轻组分浓度和塔底的液相轻组分浓度的设定值Y₁ ^*、X_n ^*，显示精馏段和提馏段温度波拐点的设定值S_r ^*、S_s ^*及当前时刻温度波拐点位置和控制参数求解模块求解出的下一时刻的控制参数，并将控制参数传递给所述控制站3，所述控制站3根据得到的控制参数对控制器进行调整，进而实现对内部热耦合精馏塔1的控制调整；所述上位机5还将以上信息传递给数据存储装置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=\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{{\mathrm{dX}}_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{{\mathrm{dX}}_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{{\mathrm{dX}}_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}(k{)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{{\mathrm{dX}}_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)

式中，α为相对挥发度，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号，n为塔底编号)，X_i(k)为k采样时刻第i块板塔的液相轻组分的浓度测量值；

所述气液相流率由公式(8)、(9)、(10)、(11)、(12)、(13)、(14)得到，公式(8)、(9)、(10)、(11)、(12)、(13)、(14)的表达式如下：

V₁(k)＝F(k)(1-q(k+1))   (8)

L_n(k)＝F(k)q(k+1)   (9)

$L_i\left(k\right)=\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{Q}_j\left(k\right)/\mathrm{\λ},(i=1,...,f-1)---\left(10\right)$

V_i+1(k)＝V₁(k)+L_i(k) (i＝1,...,f-1)   (11)

$L_{f+i-1}\left(k\right)=L_{f-1}\left(k\right)+F\left(k\right)q(k+1)-\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{Q}_j\left(k\right)/\mathrm{\λ},(i=1,...,f-2)---\left(12\right)$

$L_{f+i}\left(k\right)=V_f\left(k\right)-F\left(k\right)(1-q(k+1))-\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{Q}_j\left(k\right)/\mathrm{\λ},(i=1,...,f-2)---\left(13\right)$

Q_i(k)＝UA×(T_i(k)-T_i+f-1(k)),i＝1,...,f-1   (14)

式中，i表示塔板编号，f为进料板编号，Q_i(k)为第i块塔板上的热耦合量，UA为传热速率，λ为汽化潜热，q(k+1)为第k+1采样时刻的进料热状况，其中，q(k+1)为所述控制装置下一时刻的两个控制参数之一；

通过公式(3)-(14)，建立浓度梯度与进料热状况之间的关系，浓度梯度与进料热状况之间的关系由公式(15)简化表示，公式(15)的表达式如下：

$\frac{{\mathrm{dX}}_i\left(k\right)}{\mathrm{dt}}=f_i\left(q(k+1)\right),(i=\mathrm{1,2},...,n)---\left(15\right)$

式中，i表示塔板编号(1为塔顶编号，n为塔底编号)，f_i表示由公式(3)-(14)得到的结构已知的非线性函数关系；

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

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

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

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

将公式(15)代入公式(16)、(17)，进一步得到公式(18)、(19)，公式(18)、(19)的表达式如下：

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

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

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

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

${\hat{T}}_i\left(k\right)=\frac{T_{s1}{e}^{-{\mathrm{\γ}}_r(i-S_r\left(k\right))}+T_{s2}}{1+e^{-{\mathrm{\γ}}_r(i-S_r\left(k\right))}},i=f,f+1,...,n---\left(21\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分别表征精馏段与提馏段温度波拐点处的斜率大小；

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

$\frac{{\mathrm{dS}}_r}{\mathrm{dt}}\left(k\right)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}\left[\frac{{\mathrm{dT}}_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(22\right)$

$\frac{{\mathrm{dS}}_s}{\mathrm{dt}}\left(k\right)=\underset{i=f}{\overset{n}{\mathrm{\Σ}}}\left[\frac{{\mathrm{dT}}_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(23\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为温度波静态描述函数的常系数；

将公式(18)、(19)代入公式(22)、(23)，进一步得到公式(24)、(25)，公式(24)、(25)的表达式如下：

$\frac{{\mathrm{dS}}_r}{\mathrm{dt}}\left(k\right)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}[-f_i\left(q(k+1)\right)\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_r(k+1){\mathrm{\αbe}}^{\frac{b}{T_i\left(k\right)+c}-a}}\frac{T_{r2}-T_{r1}}{{\mathrm{\γ}}_r(T_i\left(k\right)-T_{r2})(T_i\left(k\right)-T_{r1})}]---\left(24\right)$

$\frac{{\mathrm{dS}}_s}{\mathrm{dt}}\left(k\right)=\underset{i=f}{\overset{n}{\mathrm{\Σ}}}[-f_i\left(q(k+1)\right)\frac{(\mathrm{\α}-1){(T_i\left(k\right)+c)}^2}{P_s(k+1){\mathrm{\αbe}}^{\frac{b}{T_i\left(k\right)+c}-a}}\frac{T_{s2}-T_{s1}}{{\mathrm{\γ}}_s(T_i\left(k\right)-T_{s2})(T_i\left(k\right)-T_{s1})}]---\left(25\right)$

(五)所述设定值转换模块实现设定值转换的步骤为：根据温度波特静态描述函数，通过转换公式(26)、(27)将浓度设定值转换为温度波的拐点设定值，转换公式(26)、(27)的表达式如下：

$\frac{b}{a-\ln [P_r\left(k\right)\×(\mathrm{\α}-(\mathrm{\α}-1){Y_1}^{*})]}-c=\frac{T_{r1}{e}^{-{\mathrm{\γ}}_r(1-{S_r}^{*})}+T_{r2}}{1+e^{-{\mathrm{\γ}}_r(1-{S_r}^{*})}}---\left(26\right)$

$\frac{b}{a-\ln \frac{p_s\left(k\right)}{{x_n}^{*}+(1-{X_n}^{*})/\mathrm{\α}}}-c=\frac{T_{s1}{e}^{-{\mathrm{\γ}}_s(n-{S_s}^{*})}+T_{s2}}{1+e^{-{\mathrm{\γ}}_s(n-{S_s}^{*})}}---\left(27\right)$

式中，Y₁ ^*、X_n ^*分别为塔顶的气相轻组分浓度的设定值和塔底的液相轻组分浓度的设定值，S_r ^*、S_s ^*分别为精馏段和提馏段温度波拐点的设定值，P_r(k)为k采样时刻精馏段压强、P_s(k)为k采样时刻提馏段压强，α为相对挥发度，a、b、c为安东尼常数；

(六)控制参数求解模块求解控制参数的步骤为：根据所述各个描述模块和所述设定值转换模块得到的信息求取下一时刻的控制参数，其中，控制率由公式(28)、(29)描述，公式(28)、(29)的表达式如下：

$\frac{{\mathrm{dS}}_r}{\mathrm{dt}}\left(k\right)=K_1({S_r}^{*}-S_r\left(k\right))+K_{2}T\underset{i=1}{\overset{k}{\mathrm{\Σ}}}({S_r}^{*}-S_r\left(i\right))---\left(28\right)$

$\frac{{\mathrm{dS}}_s}{\mathrm{dt}}\left(k\right)=K_3({S_s}^{*}-S_s\left(k\right))+K_{4}T\underset{i=1}{\overset{k}{\mathrm{\Σ}}}({S_s}^{*}-S_s\left(i\right))---\left(29\right)$

式中，T为采样周期，S_r(k)、S_r(i)分别为k和i采样时刻内部热耦合精馏塔精馏段温度波的拐点位置，S_s(k)、S_s(i)分别为k和i采样时刻内部热耦合精馏塔提馏段温度波的拐点位置，S_r ^*、S_s ^*分别为精馏段和提馏段温度波拐点的设定值，K₁、K₂、K₃、K₄为控制器参数参数，可根据实际的控制需求调整，其中，K₁和K₃在20-80之间取值，K₂和K₄在50-800之间取值；

将公式(24)、(25)代入到公式(28)、(29)，即可求解出所述控制装置下一时刻的控制参数，即第k+1采样时刻的进料热状况q(k+1)和精馏段压强P_r(k+1)。

2.根据权利要求1所述的基于温度波特性的内部热耦合精馏塔控制装置，其特征在于：所述上位机还用于设定采用周期T，设定K₁、K₂、K₃、K₄四个系统参数的值，以及设定塔顶的气相轻组分浓度和塔底的液相轻组分浓度的设定值Y₁ ^*、X_n ^*，显示精馏段和提馏段温度波拐点的设定值S_r ^*、S_s ^*及当前时刻温度波拐点位置和控制参数求解模块求解出的下一时刻的控制参数，并将控制参数传递给所述控制站，所述控制站根据得到的控制参数对控制器进行调整，进而实现对内部热耦合精馏塔的控制调整；所述上位机还将以上信息传递给数据存储装置。