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

摘要

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

主权项

一种基于温度波模型预测控制的内部热耦合精馏控制装置，其特征在于：包括内部热耦合精馏塔、智能仪表、控制站、数据存储装置及上位机，所述智能仪表与所述内部热耦合精馏塔连接，用于进行数据采集；所述控制站与所述内部热耦合精馏塔连接，用于实现对内部热耦合精馏塔的控制；所述数据存储装置与所述智能仪表和所述控制站连接，用于实现数据存储；所述上位机与所述数据存储装置和所述控制站连接，用于实现控制参数的求解，所述上位机包括用于观测浓度及浓度梯度的浓度梯度描述模块、用于观测温度梯度的温度梯度描述模块、用于观测温度波静态的温度波静态描述模块、用于观测温度波动态的温度波动态描述模块、用于实现设定值转换的设定值转换模块及用于求解控制参数的控制参数求解模块；其中，(一)所述浓度梯度描述模块观测浓度及浓度梯度的步骤为：通过智能仪表中的温度检测元件、压力检测元件、流量检测元件采集相应的温度、压力、流量参数，传输至数据存储装置，再由数据存储装置传输至所述浓度梯度描述模块，通过所述浓度梯度描述模块确定出浓度梯度与进料热状况之间的关系，所述浓度梯度描述模块包括各塔板的浓度观测和各塔板的浓度梯度观测两部分；1)各塔板的浓度观测，根据公式(1)、(2)获得各塔板当前时刻的浓度值，并将结果传输至数据存储装置，公式(1)、(2)的表达式如下：$X_i\left(t\right)=\frac{P_r\left(t\right)\alpha e^{\frac{T_i\left(t\right)+c}{b}-a}-1}{\alpha -1},i=\mathrm{1,2},......,f-1---\left(1\right)$$X_i\left(t\right)=\frac{P_s\left(t\right)\alpha e^{\frac{T_i\left(t\right)+c}{b}-a}-1}{\alpha -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}{\Sigma }}{Q}_j\left(t\right)/\lambda ,(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}{\Sigma }}{Q}_j\left(t\right)/\lambda ,(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}{\Sigma }}{Q}_j\left(t\right)/\lambda ,(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{(\alpha -1){(T_i\left(t\right)+c)}^2}{P_r\left(t\right)\alpha {\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{(\alpha -1){(T_i\left(t\right)+c)}^2}{P_s\left(t\right)\alpha {\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{(\alpha -1){(T_i\left(t\right)+c)}^2}{P_r\left(t\right)\alpha {\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{(\alpha -1){(T_i\left(t\right)+c)}^2}{P_s\left(t\right)\alpha {\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}^{-{\gamma }_r(i-S_r\left(k\right))}+T_{r2}}{1+e^{-{\gamma }_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}^{-{\gamma }_s(i-S_s\left(k\right))}+T_{s2}}{1+e^{-{\gamma }_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}{\Sigma }}\left[\frac{d{T}_i\left(t\right)}{\mathrm{dt}}\frac{T_{r2}-T_{r1}}{{\gamma }_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}{\Sigma }}\left[\frac{d{T}_i\left(t\right)}{\mathrm{dt}}\frac{T_{s2}-T_{s1}}{{\gamma }_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}{\Sigma }}[-f_i\left(q\left(t\right)\right)\frac{(\alpha -1){(T_i\left(t\right)+c)}^2}{P_r\left(t\right)\alpha {\mathrm{be}}^{\frac{b}{T_i\left(t\right)+c}-a}}\frac{T_{r2}-T_{r1}}{{\gamma }_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}{\Sigma }}[-f_i\left(q\left(t\right)\right)\frac{(\alpha -1){(T_i\left(t\right)+c)}^2}{P_s\left(t\right)\alpha {\mathrm{be}}^{\frac{b}{T_i\left(t\right)+c}-a}}\frac{T_{s2}-T_{s1}}{{\gamma }_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)\times (\alpha -(\alpha -1){Y_1}^{*})]}-c=\frac{T_{r1}{e}^{-{\gamma }_r(1-{S_r}^{*})}+T_{r2}}{1+e^{-{\gamma }_r(1-{S_r}^{*})}}---\left(2\right)$$\frac{b}{a-\ln \frac{p_s\left(t\right)}{{X_n}^{*}+(1-{X_n}^{*})/\alpha }}-c=\frac{T_{s1}{e}^{-{\gamma }_s(n-{S_s}^{*})}+T_{s2}}{1+e^{-{\gamma }_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}{{\gamma }_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}{{\gamma }_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)\times T---\left(30\right)$${\hat{S}}_s\left(t\right)=S_s(t-T)+\frac{d{S}_s}{\mathrm{dt}}(t-T)\times 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(\tau \right)=\underset{i=1}{\overset{f-1}{\Sigma }}[-f_i\left(q\left(\tau \right)\right)\frac{(\alpha -1){(T_i\left(\tau \right)+c)}^2}{P_r\left(\tau \right)\alpha {\mathrm{be}}^{\frac{b}{T_i\left(\tau \right)+c}-a}}\frac{T_{r2}-T_{r1}}{{\gamma }_r(T_i\left(\tau \right)-T_{r2})(T_i\left(\tau \right)-T_{r1})}]---\left(34\right)$$\frac{d{S}_s}{\mathrm{dt}}\left(\tau \right)=\underset{i=f}{\overset{n}{\Sigma }}[-f_i\left(q\left(\tau \right)\right)\frac{(\alpha -1){(T_i\left(\tau \right)+c)}^2}{P_s\left(\tau \right)\alpha {\mathrm{be}}^{\frac{b}{T_i\left(\tau \right)+c}-a}}\frac{T_{s2}-T_{s1}}{{\gamma }_s(T_i\left(\tau \right)-T_{s2})(T_i\left(\tau \right)-T_{s1})}]---\left(35\right)$式中，τ表示未来一段时间，即对于提馏段的压强保持不变，即P_s(τ)＝P_s(t)，$T_i\left(\tau \right)=T_i\left(t\right)+{\int }_t^{\tau }\frac{d{T}_i\left(t^{\prime }\right)}{\mathrm{dt}}{\mathrm{dt}}^{\prime },$而由公式(18)、(19)得到；定义：$u\left(\tau \right)=\left[\begin{array}{c}q\left(\tau \right)\\ P_r\left(\tau \right)\end{array}\right]---\left(36\right)$式中，u(τ)表示未来一段时间的控制参数，即需要求解的控制参数；将公式(34)、(35)用以下公式简化表示：$\stackrel{.}{S}\left(\tau \right)=\left[\begin{array}{c}{\stackrel{.}{S}}_r\left(\tau \right)\\ {\stackrel{.}{S}}_s\left(\tau \right)\end{array}\right]=h(u,\tau )---\left(37\right)$根据公式(32)、(33)得到的误差，定义：φ(τ)＝(S(τ)‑S^*+e(t))²+(u(τ))²   (38)这样对于控制参数的求解就可以转化为以下动态优化问题：$\mathrm{Min\; J}={\int }_t^{\tau }\phi \left(t^{\prime }\right){\mathrm{dt}}^{\prime }---\left(39\right)$$s.t.\stackrel{.}{S}\left(\tau \right)=h(u,\tau )$式中，J表示目标函数；3)实时求解优化问题，梯度信息用如下方法得到：构造函数：H(τ)＝φ(τ)+v^Th(u,τ)   (40)其中，$\frac{\mathrm{dv}}{\mathrm{d\tau }}=-\frac{\partial H}{\partial S}=-\frac{\partial \phi }{\partial S}-v^T\frac{\partial h}{\partial S}---\left(41\right)$$\stackrel{.}{S}\left(\tau \right)=\frac{\partial H}{\partial v}=h\left(\tau \right)---\left(42\right)$进而得到梯度公式的表达式：$g\left(u\right)=\frac{\partial H}{\partial u}---\left(43\right)$给定初始的控制变量u⁰、初始步长α⁰、迭代截止条件ε以及初始迭代计数k＝0，通过以下步骤完成对控制参数的实时优化：A、计算$g^k=\frac{\partial H}{\partial 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(\tau \right)=u^k=\left[\begin{array}{c}q\left(\tau \right)\\ P_r\left(\tau \right)\end{array}\right].$