CN104635493A - 基于温度波模型预测控制的内部热耦合精馏控制装置 - Google Patents

Abstract

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

Description

基于温度波模型预测控制的内部热耦合精馏控制装置

技术领域

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

背景技术

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

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

发明内容

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

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

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

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

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

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

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

2)各塔板的浓度梯度观测，根据公式(3)、(4)、(5)、(6)获得各塔板当前时刻的浓度梯度与进料热状况之间的关系，并将结果传输至温度梯度描述模块，公式(3)、(4)、(5)、(6)的表达式如下：

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

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

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

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

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

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

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

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

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

V₁(t)＝F(t)(1-q(t))    (8) 

L_n(t)＝F(t)q(t)    (9)

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

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

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

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

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

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

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

$\frac{d{X}_i\left(t\right)}{\mathrm{dt}}=f_i\left(q\left(t\right)\right)(i=\mathrm{1,2},...,n)---\left(15\right)$

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

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

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

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

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

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

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

$\frac{d{T}_i\left(t\right)}{\mathrm{dt}}=-f_i\left(q\left(t\right)\right)\frac{{(\mathrm{\α}-1)(T_i\left(t\right)+c)}^2}{P_s\left(t\right)\mathrm{\αb}{e}^{\frac{b}{T_i\left(t\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},...2f-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)$

式中，t为当前采样时刻，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号,n为塔底编号)，为t采样时刻第i块板塔的温度预测值，S_r(t)、S_s(t)分别为t采样时刻内部热耦合精馏塔精馏段、提馏段温度波的拐点，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(t\right)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}\left[\frac{{\mathrm{dT}}_i\left(t\right)}{\mathrm{dt}}\frac{T_{r2}-T_{r1}}{{\mathrm{\γ}}_r(T_i\left(t\right)-T_{r2})(T_i\left(t\right)-T_{r1})}\right]---\left(22\right)$

$\frac{{\mathrm{dS}}_s}{\mathrm{dt}}\left(t\right)=\underset{i=f}{\overset{n}{\mathrm{\Σ}}}\left[\frac{{\mathrm{dT}}_i\left(t\right)}{\mathrm{dt}}\frac{T_{s2}-T_{s1}}{{\mathrm{\γ}}_s(T_i\left(t\right)-T_{s2})(T_i\left(t\right)-T_{s1})}\right]---\left(23\right)$

式中，分别为t采样时刻精馏段和提馏段温度波的拐点移动速度，为t采样时刻第i块板塔的温度梯度值，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号，n为塔底编号)，T_i(t)为t采样时刻第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(t\right)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}[-f_i\left(q\left(t\right)\right)\frac{(\mathrm{\α}-1){(T_i\left(t\right)+c)}^2}{P_r\left(t\right)\mathrm{\αb}{e}^{\frac{b}{T_i\left(t\right)+c}-a}}\frac{T_{r2}-T_{r1}}{{\mathrm{\γ}}_r(T_i\left(t\right)-T_{r2})(T_i\left(t\right)-T_{r1})}]---\left(24\right)$

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

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

$\frac{b}{a-\ln [p_r\left(t\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(t\right)}{{X_n}^{*}+(1-{X_n}^{*})/\mathrm{\α}}}-c=\frac{T_{r1}{e}^{-{\mathrm{\γ}}_r(1-{S_r}^{*})}+T_{r2}}{1+e^{-{\mathrm{\γ}}_r(1-{S_r}^{*})}}---\left(27\right)$

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

(六)所述控制参数求解模块求解控制参数的步骤为：

1)计算当前时刻拐点位置的预测值与实际值的误差，来修正预测模型，由公式(28)、(29)得到当前时刻拐点位置的实际值，公式(28)、(29)的表达式如下：

$S_r\left(t\right)=\frac{1}{{\mathrm{\γ}}_r}\ln \frac{T_{r2}-T_1\left(t\right)}{T_1\left(t\right)-T_{r1}}+1---\left(28\right)$

$S_s\left(t\right)=\frac{1}{{\mathrm{\γ}}_s}\ln \frac{T_{s2}-T_n\left(t\right)}{T_n\left(t\right)-T_{s1}}+n---\left(29\right)$

并由公式(30)、(31)得到当前时刻拐点位置的预测值，公式(30)、(31)的表达式如下：

${\hat{S}}_r\left(t\right)=S_r(t-T)+\frac{{\mathrm{dS}}_r}{\mathrm{dt}}(t-T)\×T---\left(30\right)$

${\hat{S}}_s\left(t\right)=S_s(t-T)+\frac{{\mathrm{dS}}_s}{\mathrm{dt}}(t-T)\×T---\left(31\right)$

式中，T为采样周期，分别为t采样时刻内部热耦合精馏塔精馏段、提馏段温度波的拐点位置的预测值，S_r(t-T)、S_s(t-T)由t-T采样时刻的公式(28)、(29)得到， 由t-T采样时刻的公式(24)、(25)得到；

误差由公式(32)、(33)得到，公式(32)、(33)的表达式如下：

$e_r\left(t\right)=S_r\left(t\right)-{\hat{S}}_r\left(t\right)---\left(32\right)$

$e_s\left(t\right)=S_s\left(t\right)-{\hat{S}}_s\left(t\right)---\left(33\right)$

将公式(32)、(33)两式记作 $e\left(t\right)=\left[\begin{array}{c}{e}_r\left(t\right)\\ e_s\left(t\right)\end{array}\right];$

2)建立求解控制参数的预测模型，由公式(24)、(25)可得：

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

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

式中，τ表示未来一段时间，即对于提馏段的压强保持不变，即P_s(τ)＝P_s(t)，  $T_i\left(\mathrm{\τ}\right)=T_i\left(t\right)+{\∫}_t^{\mathrm{\τ}}\frac{d{T}_i\left(t^{\′}\right)}{\mathrm{dt}}{\mathrm{dt}}^{\′},$ 而由公式(18)、(19)得到；

定义：

$\left(\mathrm{\τ}\right)=\left[\begin{array}{c}q\left(\mathrm{\τ}\right)\\ P_r\left(\mathrm{\τ}\right)\end{array}\right]---\left(36\right)$

式中，u(τ)表示未来一段时间的控制参数，即需要求解的控制参数；

将公式(34)、(35)用以下公式简化表示：

$\stackrel{.}{s}\left(\mathrm{\τ}\right)=\left[\begin{array}{c}{\stackrel{.}{S}}_r\left(\mathrm{\τ}\right)\\ {\stackrel{.}{S}}_s\left(\mathrm{\τ}\right)\end{array}\right]=h(u,\mathrm{\τ})---\left(37\right)$

根据公式(32)、(33)得到的误差，定义：

φ(τ)＝(S(τ)-S^*+e(t))²+(u(τ))²    (38) 

这样对于控制参数的求解就可以转化为以下动态优化问题：

$\mathrm{MinJ}={\∫}_t^{\mathrm{\τ}}\mathrm{\φ}\left(t^{\′}\right)d{t}^{\′}---\left(39\right)$

$s.t.\stackrel{.}{S}\left(\mathrm{\τ}\right)=h(u,\mathrm{\τ})$

式中，J表示目标函数；

3)实时求解优化问题，梯度信息用如下方法得到：

构造函数： 

H(τ)＝φ(τ)+v^Th(u,τ)    (40) 

其中，

$\frac{\mathrm{dv}}{\mathrm{dt}}=-\frac{\∂H}{\∂S}=-\frac{\∂\mathrm{\φ}}{\∂S}-v^T\frac{\∂h}{\∂S}---\left(41\right)$

$\stackrel{.}{S}\left(\mathrm{\τ}\right)=\frac{\∂H}{\∂v}=h\left(\mathrm{\τ}\right)---\left(42\right)$

进而得到梯度公式的表达式：

$g\left(u\right)=\frac{\∂H}{\∂u}---\left(43\right)$

给定初始的控制变量u⁰、初始步长α⁰、迭代截止条件ε以及初始迭代计数k＝0，通过以下步骤完成对控制参数的实时优化：

A、计算

B、如果k＝0，跳到第3步；否则，将u^k代入目标函数J，如果|J^k-J^k+1|≤ε，停止迭代并输出u^k，如果|J^k-J^k+1|＞ε，则计算其中s^k-1＝u^k-u^k-1，y^k-1＝g^k-g^k-1；

C、计算u^k+1＝u^k+α^k·(-g^k)；

D、增加迭代计数k＝k+1，并返回第1步进行下一次迭代；

其中，上标k表示迭代计数，通过以上步骤即可得到控制参数

作为优选，所述上位机还用于设定塔顶的气相轻组分浓度和塔底的液相轻组分浓度的设定值Y₁ ^*、X_n ^*，以及设定初始的控制变量u⁰、初始步长α⁰、迭代截止条件ε，显示当前时刻温度波拐点位置和求解出的下一段时间的控制参数，并将控制参数传递给所述控制站，所述控制站根据得到的控制参数对控制器进行调整，进而实现对内部热耦合精馏塔的控制调整；所述上位机还将以上信息传递给数据存储装置，方便操作人员查阅历史记录，提高生产控制品质。

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

附图说明

附图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(t\right)=\frac{P_r\left(t\right)\mathrm{\α}{e}^{\frac{T_i\left(t\right)+c}{b}-a}-1}{\mathrm{\α}-1},i=\mathrm{1,2},......,f-1---\left(1\right)$

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

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

2)各塔板的浓度梯度观测，根据公式(3)、(4)、(5)、(6)获得各塔板当前时刻的浓度梯度与进料热状况之间的关系，并将结果传输至温度梯度描述模块，公式(3)、(4)、(5)、(6)的表达式如下：

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

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

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

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

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

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

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

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

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

V₁(t)＝F(t)(1-q(t))    (8) 

L_n(t)＝F(t)q(t)    (9)

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

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

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

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

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

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

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

$\frac{d{X}_i\left(t\right)}{\mathrm{dt}}=f_i\left(q\left(t\right)\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{d{T}_i\left(t\right)}{\mathrm{dt}}=-\frac{d{X}_i\left(t\right)}{\mathrm{dt}}\frac{{(\mathrm{\α}-1)(T_i\left(t\right)+c)}^2}{P_r\left(t\right)\mathrm{\αb}{e}^{\frac{b}{T_i\left(t\right)+c}-a}},i=\mathrm{1,2},......,f-1---\left(16\right)$

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

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

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

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

$\frac{d{T}_i\left(t\right)}{\mathrm{dt}}=-f_i\left(q\left(t\right)\right)\frac{{(\mathrm{\α}-1)(T_i\left(t\right)+c)}^2}{P_s\left(t\right)\mathrm{\αb}{e}^{\frac{b}{T_i\left(t\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},...2f-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)$

式中，t为当前采样时刻，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号,n为塔底编号)，为t采样时刻第i块板塔的温度预测值，S_r(t)、S_s(t)分别为t采样时刻内部热耦合精馏塔精馏段、提馏段温度波的拐点，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}{\mathrm{dt}}\left(t\right)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}\left[\frac{{\mathrm{dT}}_i\left(t\right)}{\mathrm{dt}}\frac{T_{r2}-T_{r1}}{{\mathrm{\γ}}_r(T_i\left(t\right)-T_{r2})(T_i\left(t\right)-T_{r1})}\right]---\left(22\right)$

$\frac{{\mathrm{dS}}_s}{\mathrm{dt}}\left(t\right)=\underset{i=f}{\overset{n}{\mathrm{\Σ}}}\left[\frac{{\mathrm{dT}}_i\left(t\right)}{\mathrm{dt}}\frac{T_{s2}-T_{s1}}{{\mathrm{\γ}}_s(T_i\left(t\right)-T_{s2})(T_i\left(t\right)-T_{s1})}\right]---\left(23\right)$

式中，分别为t采样时刻精馏段和提馏段温度波的拐点移动速度，为t采样时刻第i块板塔的温度梯度值，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号，n为塔底编号)，T_i(t)为t采样时刻第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(t\right)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}[-f_i\left(q\left(t\right)\right)\frac{(\mathrm{\α}-1){(T_i\left(t\right)+c)}^2}{P_r\left(t\right)\mathrm{\αb}{e}^{\frac{b}{T_i\left(t\right)+c}-a}}\frac{T_{r2}-T_{r1}}{{\mathrm{\γ}}_r(T_i\left(t\right)-T_{r2})(T_i\left(t\right)-T_{r1})}]---\left(24\right)$

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

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

$\frac{b}{a-\ln [p_r\left(t\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(t\right)}{{X_n}^{*}+(1-{X_n}^{*})/\mathrm{\α}}}-c=\frac{T_{r1}{e}^{-{\mathrm{\γ}}_r(1-{S_r}^{*})}+T_{r2}}{1+e^{-{\mathrm{\γ}}_r(1-{S_r}^{*})}}---\left(27\right)$

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

(六)所述控制参数求解模块11求解控制参数的步骤为：

1)计算当前时刻拐点位置的预测值与实际值的误差，来修正预测模型，由公式(28)、(29)得到当前时刻拐点位置的实际值，公式(28)、(29)的表达式如下：

$S_r\left(t\right)=\frac{1}{{\mathrm{\γ}}_r}\ln \frac{T_{r2}-T_1\left(t\right)}{T_1\left(t\right)-T_{r1}}+1---\left(28\right)$

$S_s\left(t\right)=\frac{1}{{\mathrm{\γ}}_s}\ln \frac{T_{s2}-T_n\left(t\right)}{T_n\left(t\right)-T_{s1}}+n---\left(29\right)$

并由公式(30)、(31)得到当前时刻拐点位置的预测值，公式(30)、(31)的表达式如下：

${\hat{S}}_r\left(t\right)=S_r(t-T)+\frac{{\mathrm{dS}}_r}{\mathrm{dt}}(t-T)\×T---\left(30\right)$

${\hat{S}}_s\left(t\right)=S_s(t-T)+\frac{{\mathrm{dS}}_s}{\mathrm{dt}}(t-T)\×T---\left(31\right)$

式中，T为采样周期，分别为t采样时刻内部热耦合精馏塔精馏段、提馏段温度波的拐点位置的预测值，S_r(t-T)、S_s(t-T)由t-T采样时刻的公式(28)、(29)得到， 由t-T采样时刻的公式(24)、(25)得到；

误差由公式(32)、(33)得到，公式(32)、(33)的表达式如下：

$e_r\left(t\right)=S_r\left(t\right)-{\hat{S}}_r\left(t\right)---\left(32\right)$

$e_s\left(t\right)=S_s\left(t\right)-{\hat{S}}_s\left(t\right)---\left(33\right)$

将公式(32)、(33)两式记作 $e\left(t\right)=\left[\begin{array}{c}{e}_r\left(t\right)\\ e_s\left(t\right)\end{array}\right];$

2)建立求解控制参数的预测模型，由公式(24)、(25)可得：

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

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

式中，τ表示未来一段时间，即对于提馏段的压强保持不变，即P_s(τ)＝P_s(t)，  $T_i\left(\mathrm{\τ}\right)=T_i\left(t\right)+{\∫}_t^{\mathrm{\τ}}\frac{d{T}_i\left(t^{\′}\right)}{\mathrm{dt}}{\mathrm{dt}}^{\′},$ 而由公式(18)、(19)得到；

定义：

$\left(\mathrm{\τ}\right)=\left[\begin{array}{c}q\left(\mathrm{\τ}\right)\\ P_r\left(\mathrm{\τ}\right)\end{array}\right]---\left(36\right)$

式中，u(τ)表示未来一段时间的控制参数，即需要求解的控制参数；

将公式(34)、(35)用以下公式简化表示：

$\stackrel{.}{s}\left(\mathrm{\τ}\right)=\left[\begin{array}{c}{\stackrel{.}{S}}_r\left(\mathrm{\τ}\right)\\ {\stackrel{.}{S}}_s\left(\mathrm{\τ}\right)\end{array}\right]=h(u,\mathrm{\τ})---\left(37\right)$

根据公式(32)、(33)得到的误差，定义：

φ(τ)＝(S(τ)-S^*+e(t))²+(u(τ))²    (38) 

这样对于控制参数的求解就可以转化为以下动态优化问题：

$\mathrm{MinJ}={\∫}_t^{\mathrm{\τ}}\mathrm{\φ}\left(t^{\′}\right)d{t}^{\′}---\left(39\right)$

$s.t.\stackrel{.}{S}\left(\mathrm{\τ}\right)=h(u,\mathrm{\τ})$

式中，J表示目标函数；

3)实时求解优化问题，梯度信息用如下方法得到：

构造函数： 

H(τ)＝φ(τ)+v^Th(u,τ)    (40) 

其中，

$\frac{\mathrm{dv}}{\mathrm{dt}}=-\frac{\∂H}{\∂S}=-\frac{\∂\mathrm{\φ}}{\∂S}-v^T\frac{\∂h}{\∂S}---\left(41\right)$

$\stackrel{.}{S}\left(\mathrm{\τ}\right)=\frac{\∂H}{\∂v}=h\left(\mathrm{\τ}\right)---\left(42\right)$

进而得到梯度公式的表达式：

$g\left(u\right)=\frac{\∂H}{\∂u}---\left(43\right)$

给定初始的控制变量u⁰、初始步长α⁰、迭代截止条件ε以及初始迭代计数k＝0，通过以下步骤完成对控制参数的实时优化：

A、计算

B、如果k＝0，跳到第3步；否则，将u^k代入目标函数J，如果|J^k-J^k+1|≤ε，停止迭代并输出u^k，如果|J^k-J^k+1|＞ε，则计算其中s^k-1＝u^k-u^k-1，y^k-1＝g^k-g^k-1；

C、计算u^k+1＝u^k+α^k·(-g^k)；

D、增加迭代计数k＝k+1，并返回第1步进行下一次迭代；

其中，上标k表示迭代计数，通过以上步骤即可得到控制参数

本实施例中，所述上位机5还用于设定塔顶的气相轻组分浓度和塔底的液相轻组分浓度的设定值Y₁ ^*、X_n ^*，以及设定初始的控制变量u⁰、初始步长α⁰、迭代截止条件ε，显示当前时刻温度波拐点位置和求解出的下一段时间的控制参数，并将控制参数传递给所述控制站3，所述控制站3根据得到的控制参数对控制器进行调整，进而实现对内部热耦合精馏塔1的控制调整；所述上位机5还将以上信息传递给数据存储装置4，方便操作人员查阅历史记录，提高生产控制品质。

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

Claims (2)

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

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

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

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

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

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

2)各塔板的浓度梯度观测，根据公式(3)、(4)、(5)、(6)获得各塔板当前时刻的浓度梯度与进料热状况之间的关系，并将结果传输至温度梯度描述模块，公式(3)、(4)、(5)、(6)的表达式如下：

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

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

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

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

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

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

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

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

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

V₁(t)＝F(t)(1-q(t))   (8)

L_n(t)＝F(t)q(t)   (9)

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

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

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

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

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

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

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

$\frac{{\mathrm{dX}}_i\left(t\right)}{\mathrm{dt}}=f_i\left(q\left(t\right)\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(t\right)}{\mathrm{dt}}=-\frac{{\mathrm{dX}}_i\left(t\right)}{\mathrm{dt}}\frac{(\mathrm{\α}-1){(T_i\left(t\right)+c)}^2}{P_r\left(t\right)\mathrm{\α}{\mathrm{be}}^{\frac{b}{T_i\left(t\right)+c}-a}},i=\mathrm{1,2},......,f-1---\left(16\right)$

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

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

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

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

$\frac{{\mathrm{dT}}_i\left(t\right)}{\mathrm{dt}}=-f_i\left(q\left(t\right)\right)\frac{(\mathrm{\α}-1){(T_i\left(t\right)+c)}^2}{P_s\left(t\right)\mathrm{\α}{\mathrm{be}}^{\frac{b}{T_i\left(t\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)$

式中，t为当前采样时刻，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号,n为塔底编号)，为t采样时刻第i块板塔的温度预测值，S_r(t)、S_s(t)分别为t采样时刻内部热耦合精馏塔精馏段、提馏段温度波的拐点，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{d{S}_r}{\mathrm{dt}}\left(t\right)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}\left[\frac{d{T}_i\left(t\right)}{\mathrm{dt}}\frac{T_{r2}-T_{r1}}{{\mathrm{\γ}}_r(T_i\left(t\right)-T_{r2})(T_i\left(t\right)-T_{r1})}\right]---\left(22\right)$

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

式中，分别为t采样时刻精馏段和提馏段温度波的拐点移动速度，为t采样时刻第i块板塔的温度梯度值，i表示塔板编号(i＝1,2,...,f,f+1,...,n，1为塔顶编号，f为进料板编号，n为塔底编号)，T_i(t)为t采样时刻第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(t\right)=\underset{i=1}{\overset{f-1}{\mathrm{\Σ}}}[-f_i\left(q\left(t\right)\right)\frac{(\mathrm{\α}-1){(T_i\left(t\right)+c)}^2}{P_r\left(t\right)\mathrm{\α}{\mathrm{be}}^{\frac{b}{T_i\left(t\right)+c}-a}}\frac{T_{r2}-T_{r1}}{{\mathrm{\γ}}_r(T_i\left(t\right)-T_{r2})(T_i\left(t\right)-T_{r1})}]---\left(24\right)$

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

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

$\frac{b}{a-\ln [P_r\left(t\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(2\right)$

$\frac{b}{a-\ln \frac{p_s\left(t\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(t)为t采样时刻精馏段压强，P_s(t)为t采样时刻提馏段压强，，α为相对挥发度，a、b、c为安东尼常数；

(六)所述控制参数求解模块求解控制参数的步骤为：

1)计算当前时刻拐点位置的预测值与实际值的误差，来修正预测模型，由公式(28)、(29)得到当前时刻拐点位置的实际值，公式(28)、(29)的表达式如下：

$S_r\left(t\right)=\frac{1}{{\mathrm{\γ}}_r}\ln \frac{T_{r2}-T_1\left(t\right)}{T_1\left(t\right)-T_{r1}}+1---\left(28\right)$

$S_s\left(t\right)=\frac{1}{{\mathrm{\γ}}_s}\ln \frac{T_{s2}-T_n\left(t\right)}{T_n\left(t\right)-T_{s1}}+n---\left(29\right)$

并由公式(30)、(31)得到当前时刻拐点位置的预测值，公式(30)、(31)的表达式如下：

${\hat{S}}_r\left(t\right)=S_r(t-T)+\frac{d{S}_r}{\mathrm{dt}}(t-T)\×T---\left(30\right)$

${\hat{S}}_s\left(t\right)=S_s(t-T)+\frac{d{S}_s}{\mathrm{dt}}(t-T)\×T---\left(31\right)$

式中，T为采样周期，分别为t采样时刻内部热耦合精馏塔精馏段、提馏段温度波的拐点位置的预测值，S_r(t-T)、S_s(t-T)由t-T采样时刻的公式(28)、(29)得到，由t-T采样时刻的公式(24)、(25)得到；

误差由公式(32)、(33)得到，公式(32)、(33)的表达式如下：

$e_r\left(t\right)=S_r\left(t\right)-{\hat{S}}_r\left(t\right)---\left(32\right)$

$e_s\left(t\right)=S_s\left(t\right)-{\hat{S}}_s\left(t\right)---\left(33\right)$ 将公式(32)、(33)两式记作 $e\left(t\right)=\left[\begin{array}{c}{e}_r\left(t\right)\\ e_s\left(t\right)\end{array}\right];$

2)建立求解控制参数的预测模型，由公式(24)、(25)可得：

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

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

式中，τ表示未来一段时间，即对于提馏段的压强保持不变，即P_s(τ)＝P_s(t)， $T_i\left(\mathrm{\τ}\right)=T_i\left(t\right)+{\∫}_t^{\mathrm{\τ}}\frac{d{T}_i\left(t^{\′}\right)}{\mathrm{dt}}{\mathrm{dt}}^{\′},$ 而由公式(18)、(19)得到；

定义：

$u\left(\mathrm{\τ}\right)=\left[\begin{array}{c}q\left(\mathrm{\τ}\right)\\ P_r\left(\mathrm{\τ}\right)\end{array}\right]---\left(36\right)$

式中，u(τ)表示未来一段时间的控制参数，即需要求解的控制参数；

将公式(34)、(35)用以下公式简化表示：

$\stackrel{.}{S}\left(\mathrm{\τ}\right)=\left[\begin{array}{c}{\stackrel{.}{S}}_r\left(\mathrm{\τ}\right)\\ {\stackrel{.}{S}}_s\left(\mathrm{\τ}\right)\end{array}\right]=h(u,\mathrm{\τ})---\left(37\right)$

根据公式(32)、(33)得到的误差，定义：

φ(τ)＝(S(τ)-S^*+e(t))²+(u(τ))²   (38)

这样对于控制参数的求解就可以转化为以下动态优化问题：

$\mathrm{Min\; J}={\∫}_t^{\mathrm{\τ}}\mathrm{\φ}\left(t^{\′}\right){\mathrm{dt}}^{\′}---\left(39\right)$

$s.t.\stackrel{.}{S}\left(\mathrm{\τ}\right)=h(u,\mathrm{\τ})$

式中，J表示目标函数；

3)实时求解优化问题，梯度信息用如下方法得到：

构造函数：

H(τ)＝φ(τ)+v^Th(u,τ)   (40)

其中，

$\frac{\mathrm{dv}}{\mathrm{d\τ}}=-\frac{\∂H}{\∂S}=-\frac{\∂\mathrm{\φ}}{\∂S}-v^T\frac{\∂h}{\∂S}---\left(41\right)$

$\stackrel{.}{S}\left(\mathrm{\τ}\right)=\frac{\∂H}{\∂v}=h\left(\mathrm{\τ}\right)---\left(42\right)$

进而得到梯度公式的表达式：

$g\left(u\right)=\frac{\∂H}{\∂u}---\left(43\right)$

给定初始的控制变量u⁰、初始步长α⁰、迭代截止条件ε以及初始迭代计数k＝0，通过以下步骤完成对控制参数的实时优化：

A、计算 $g^k=\frac{\∂H}{\∂u}{|}_{u=u^k};$

B、如果k＝0，跳到第3步；否则，将u^k代入目标函数J，如果|J^k-J^k+1|≤ε，停止迭代并输出u^k，如果|J^k-J^k+1|＞ε，则计算其中s^k-1＝u^k-u^k-1，y^k-1＝g^k-g^k-1；

C、计算u^k+1＝u^k+α^k·(-g^k)；

D、增加迭代计数k＝k+1，并返回第1步进行下一次迭代；

其中，上标k表示迭代计数，通过以上步骤即可得到控制参数 $u\left(\mathrm{\τ}\right)=u^k=\left[\begin{array}{c}q\left(\mathrm{\τ}\right)\\ P_r\left(\mathrm{\τ}\right)\end{array}\right].$

2.根据权利要求1所述的基于温度波模型预测控制的内部热耦合精馏控制装置，其特征在于：所述上位机还用于设定塔顶的气相轻组分浓度和塔底的液相轻组分浓度的设定值Y₁ ^*、X_n ^*，以及设定初始的控制变量u⁰、初始步长α⁰、迭代截止条件ε，显示当前时刻温度波拐点位置和求解出的下一段时间的控制参数，并将控制参数传递给所述控制站，所述控制站根据得到的控制参数对控制器进行调整，进而实现对内部热耦合精馏塔的控制调整；所述上位机还将以上信息传递给数据存储装置。