2022, Vol. 45
表1(Table 1)
2a. 重庆大学 化学化工学院, 重庆 400044;
2b. 重庆大学 化工过程强化与反应国家地方联合工程实验室, 重庆 400044;
2c. 重庆大学 低品位能源利用技术及系统教育部重点实验室, 重庆 400044
2a. School of Chemistry and Chemical Engineering, Ministry of Education, Chongqing University, Chongqing 400044, P. R. China;
2b. National-Municipal Joint Engineering Laboratory for Chemical Process Intensification and Reaction, Ministry of Education, Chongqing University, Chongqing 400044, P. R. China;
2c. Key Laboratory of Low-grade Energy Utilization Technologies & Systems, Ministry of Education, Chongqing University, Chongqing 400044, P. R. China
模型预测控制(model predictive control，MPC)由于能够很好地处理多变量约束控制问题，自20世纪70年代以来已被广泛应用于石油化工生产过程中[1-2]。MPC是一种基于模型的控制算法，模型是否能够准确预测过程的动态特征是MPC实现良好控制性能的内在关键因素，MPC控制器相关参数的整定是影响控制性能的外在主要因素。
MPC中多变量之间的耦合，使采样时间、预测时域、控制时域、输入和输出变量的权重参数整定极具挑战性。与过程动态特征相关的采样时间、预测时域和控制时域一般可通过分析过程动态响应特征得到，而输入和输出变量的权重参数之间由于存在相互耦合，很难通过分析的方法获得最优值，因此，如何整定权重参数对MPC取得良好控制性能具有重要意义。目前主流的MPC权重参数整定方法可分为两类[3]：1)将过程模型等价为一阶、二阶模型，应用动态响应性能分析的方法获得权重参数[4-7]；2)将权重参数整定转化为优化问题，以优化控制性能为目标而获得最优权重参数[8-12]。前者较适于单输入单输出MPC的权重参数整定，而后者由于更易于实施，且适用于多变量系统，不需要操作者具备良好的控制理论知识而备受关注。Yamashita等[12-13]以最小化每个输出变量的参考轨迹跟踪积分方差为优化目标，应用多目标优化的方法获得输出变量和输入变量变化速率的最优权重值。Giraldo等[14]提出了基于双层优化的广义预测控制器权重参数整定方法。Gutiérrez-Urquídez等[15]将多目标优化应用于拉盖尔多项式参数化的MPC权重参数的整定。Santamaría等[16]及Vallerio等[17]则将多目标优化应用于非线性MPC权重参数的整定。
多目标优化权重参数整定方法通常为双层结构，其中内层为整定时域内每个采样时间点的MPC子优化问题的求解，而外层为整定时域内以权重参数为决策变量的控制性能目标函数的优化。这种双层结构的多目标优化问题的求解不仅存在内外层优化问题之间的迭代，而且需要大量的随机样本点来评估目标函数以获得最优Pareto解，计算结构复杂，时间成本较高。Tran等[18]将无约束MPC优化问题的解析解引入MPC权重参数整定优化问题中，提出了基于单层非线性规划的MPC权重参数整定方法，可快速获得最优权重参数，减小计算时间成本，然而该方法却不适用于有约束MPC权重参数的整定。为此，笔者提出了一种针对有约束MPC权重参数整定的非线性规划方法，该方法将MPC子优化问题转化为MPC权重参数整定优化问题的一阶最优KKT(Karush-Kuhn-Tucker)条件，可快速获得MPC的最优权重参数，与其他多目标优化方法相比可显著降低计算耗时成本。
1 模型预测控制算法 1.1 被控对象模型线性模型由于其鲁棒性能较强、计算负荷较小，被广泛应用于过程工业中MPC的实施，本研究的对象是方程(1)所示的离散线性时不变状态空间模型。
| $ \left\{\begin{array}{l} \boldsymbol{x}_{k+1}=\boldsymbol{A}_{0} \boldsymbol{x}_{k}+\boldsymbol{B}_{0} u_{k}+\boldsymbol{B}_{\mathrm{d} 0} \boldsymbol{d}_{k}+\boldsymbol{w} ; \\ \boldsymbol{y}_{k+1}=\boldsymbol{C}_{0} \boldsymbol{x}_{k+1}+\boldsymbol{v}。\end{array}\right. $ | (1) |
式中：xk∈Rn、uk∈Rm、yk∈Rq和dk∈Rl分别为第k个时间采样点的状态、输入、输出和可测扰动变量维度；A0∈Rn×n、B0∈Rn×m、C0∈Rq×n和Bd0∈Rn×l分别为状态、输入、输出和可测扰动变量的转化矩阵，其中n、m、q和l均为向量或矩阵的行数或列数；w和v分别为状态和可测输出的不确定性噪音。此处假设方程(1)的矩阵对(A0，B0)和(A0，C0)分别可控和可观测，且w和v服从高斯分布。
1.2 增广模型在方程(1)中引入差分算子后可得：
| $ \left\{\begin{array}{l} \Delta \boldsymbol{x}_{k+1}=\boldsymbol{A}_{0} \Delta \boldsymbol{x}_{k}+\boldsymbol{B}_{0} \Delta \boldsymbol{u}_{k}+\boldsymbol{B}_{\mathrm{d} 0} \Delta \boldsymbol{d}_{k}+\boldsymbol{w} ; \\ \Delta \boldsymbol{y}_{k+1}=\boldsymbol{C}_{0} \Delta \boldsymbol{x}_{k+1}+\boldsymbol{v}。\end{array}\right. $ | (2) |
式中：Δxk+1=xk+1-xk，Δyk+1=yk+1-yk，Δuk=uk-uk-1，Δdk=dk-dk-1。
将输出维度yk集成到状态维度Δxk中，得到新的状态空间模型：
| $ \left\{\begin{array}{l} \boldsymbol{z}_{k+1}=\boldsymbol{A} \boldsymbol{z}_{k}+\boldsymbol{B} \Delta \boldsymbol{u}_{k}+\boldsymbol{B}_{\mathrm{d}} \Delta \boldsymbol{d}_{k}+\boldsymbol{\xi} ; \\ \boldsymbol{y}_{k}=\boldsymbol{C} \boldsymbol{z}_{k} 。\end{array}\right. $ | (3) |
式中：
方程(1)所示的被控对象模型通常是通过系统辨识或非线性机理模型线性化得到，其状态变量通常不可测，故需在每个采样时间点以当前测得的输出值和工艺扰动量估计当前的实际状态。以方程(4)所示的状态观测器估计当前的状态：
| $ \hat{\boldsymbol{z}}_{k}=\boldsymbol{A} \hat{\boldsymbol{z}}_{k-1}+\boldsymbol{B} \Delta \boldsymbol{u}_{k}+\boldsymbol{B}_{\mathrm{d}} \Delta \boldsymbol{d}_{k}+\boldsymbol{K}\left(\boldsymbol{y}_{k}^{\mathrm{m}}-\boldsymbol{C} \hat{\boldsymbol{z}}_{k-1}\right) 。$ | (4) |
式中：
| $ \left\{\begin{array}{l} \boldsymbol{G}=\boldsymbol{A}\left[\boldsymbol{G}-\boldsymbol{G} \boldsymbol{C}^{\mathrm{T}}\left(\boldsymbol{C} \boldsymbol{GC}+\boldsymbol{R}_{\mathrm{v}}\right)^{-1} \boldsymbol{C} \boldsymbol{G}\right] \boldsymbol{A}^{\mathrm{T}}+\boldsymbol{Q}_{\mathrm{w}} ; \\ \boldsymbol{K}=\boldsymbol{A} \boldsymbol{G} \boldsymbol{C}^{\mathrm{T}}\left(\boldsymbol{C} \boldsymbol{G} \boldsymbol{C}+\boldsymbol{R}_{\mathrm{v}}\right)^{-1}。\end{array}\right. $ | (5) |
式中：Qw和Rv分别为方程(3)的状态与输出维度的噪音方差矩阵，本研究中分别取Rv=Iq×q和Qw=CTRvC。
1.4 模型预测控制器基于以上方程，MPC控制器的优化问题可描述为：
| $ \begin{array}{c} \begin{array}{l} \min \limits_{\Delta \boldsymbol{u}_{k}, \cdots \Delta \boldsymbol{u}_{k+N_{\rm c}-1}} J_{k}=\sum\limits_{i=0}^{N_{\mathrm{p}}}\left(\hat{\boldsymbol{y}}_{k+i}-\boldsymbol{y}_{k}^{\mathrm{sp}}\right)^{\mathrm{T}} \boldsymbol{Q}_{\mathrm{t}}\left(\hat{\boldsymbol{y}}_{k+i}-\boldsymbol{y}_{k}^{\mathrm{sp}}\right)+\sum\limits_{i=0}^{N_{\mathrm{c}}-1} \Delta \boldsymbol{u}_{k+i}^{\mathrm{T}} \boldsymbol{R}_{\mathrm{t}} \Delta \boldsymbol{u}_{k+i}\\ \ \ \ \ {\text{subject to:}} \end{array}\\ \begin{gathered} \hat{\boldsymbol{z}}_{k+i+1}=\boldsymbol{A} \hat{\boldsymbol{z}}_{k+i}+\boldsymbol{B} \Delta \boldsymbol{u}_{k+i}+\boldsymbol{B}_{\mathrm{d}} \Delta \boldsymbol{d}_{k+i}, \quad i=0,1, \cdots, N_{\mathrm{p}}-1 ;\\ \hat{\boldsymbol{y}}_{k+i}=\boldsymbol{C} \hat{\boldsymbol{z}}_{k+i}, \quad i=0,1, \cdots, N_{\mathrm{p}} ; \\ \boldsymbol{u}_{k+i}=\boldsymbol{u}_{k+i-1}+\Delta \boldsymbol{u}_{k+i}, \quad i=0,1, \cdots, N_{\mathrm{c}}-1 ; \\ \boldsymbol{u}_{k+i}=\boldsymbol{u}_{k+i-1}, \quad i=N_{\mathrm{c}}, \cdots, N_{\mathrm{p}}-1 ; \\ \boldsymbol{u}_{\min } \leqslant \boldsymbol{u}_{k+i} \leqslant \boldsymbol{u}_{\max }, \quad i=0,1, \cdots, N_{\mathrm{c}}-1 ; \\ \Delta \boldsymbol{u}_{\min } \leqslant \Delta \boldsymbol{u}_{k+i} \leqslant \Delta \boldsymbol{u}_{\max }, \quad i=0,1, \cdots, N_{\mathrm{c}}-1 ; \\ \boldsymbol{Q}_{\mathrm{t}}=\operatorname{diag}(\boldsymbol{Q}), \boldsymbol{R}_{\mathrm{t}}=\operatorname{diag}(\boldsymbol{R}) 。\end{gathered} \end{array} $ | (6) |
式中：Q和R分别为q个输出权重参数和m个输入变化速率权重组成的对角矩阵；Np和Nc分别为预测时域和控制时域；yksp为tk时刻的输出变量(控制变量)设定值，并假设在整个预测时域内为定值。
2 模型预测控制器权重参数整定 2.1 双层多目标优化整定方法(MOO)在MPC控制器的参数整定过程中，本研究中仅考虑了MPC控制器输出变量权重参数Q和输入变量变化速率权重参数R的最优选择，预测时域、控制时域以及采样时间均通过分析过程的动力学特征得到。Q越大则MPC控制器能够更快速地跟踪参考轨迹设定点；R越大则MPC动态响应越趋于平滑，但不能快速地跟踪参考轨迹设定点。为此，将整定时域内的MPC控制器跟踪参考轨迹的平均积分方差(Ψ1)和输入变量变化速率的二次方的平均值(Ψ2)作为优化目标，构成如下所示的权重参数整定多目标优化问题[19-20]：
| $ \begin{gathered} \begin{array}{l} \begin{gathered} \min \limits_{Q, \boldsymbol{R}} \varPsi_{1}=\frac{1}{N_{\mathrm{t}}} \sum\limits_{j=1}^{N_{\mathrm{t}}}\left(\boldsymbol{y}_{j}^{\mathrm{m}}-\boldsymbol{y}_{j}^{\mathrm{sp}}\right)^{\mathrm{T}}\left(\boldsymbol{y}_{j}^{\mathrm{m}}-\boldsymbol{y}_{j}^{\mathrm{sp}}\right), \\ \min \limits_{Q, \boldsymbol{R}} \varPsi_{2}=\frac{1}{N_{\mathrm{t}}} \sum\limits_{j=1}^{N_{\mathrm{t}}} \Delta \boldsymbol{u}_{j, 0}^{\mathrm{T}} \Delta \boldsymbol{u}_{j, 0}, \end{gathered}\\ {\text{subject to:}} \end{array}\\ \boldsymbol{Q}_{\text {min }} \leqslant \boldsymbol{Q} \leqslant \boldsymbol{Q}_{\text {max }}, \\ \boldsymbol{R}_{\min } \leqslant \boldsymbol{R} \leqslant \boldsymbol{R}_{\max }, \\ \left.\begin{array}{l} \boldsymbol{x}_{j}^{\mathrm{m}}=\boldsymbol{A}_{0} \boldsymbol{x}_{j-1}^{\mathrm{m}}+\boldsymbol{B}_{0} \boldsymbol{u}_{j-1,0}+\boldsymbol{B}_{\mathrm{d} 0} \boldsymbol{d}_{j-1}, \\ \boldsymbol{y}_{j}^{\mathrm{m}}=\boldsymbol{C}_{0} \boldsymbol{x}_{j}^{\mathrm{m}}, \\ \boldsymbol{u}_{j, 0}=\boldsymbol{u}_{j-1,0}+\Delta \boldsymbol{u}_{j, 0}, \\ \Delta \boldsymbol{u}_{j, 0} \in \min \limits_{\Delta \boldsymbol{u}_{j, 0} \cdots \boldsymbol{u}_{j}, {N}_{\mathrm{c}-1}} J_{j}, \end{array}\right\} j=1,2, \cdots, N_{\mathrm{t} }。\end{gathered} $ | (7) |
式中：xjm为被控对象模型状态变量；yjm和yjsp分别为tj时刻工艺输出值和参考轨迹设定值；Nt为选取的整定时域采样点个数。在每个采样时间点tj，MPC控制器需要求解方程(6)所示的MPC子优化问题，并将得到的控制动作序列中的第一个值(Δuj, 0)作为被控对象模型的输入以得到输出变量yjm。
方程(7)所示的多目标优化问题应用Shama等[21]开发的基于MS Excel的多目标优化求解器(EMOO)求解，该求解器采用NSGA-II算法[22]，能够很好地处理含有约束和整型变量的多目标优化问题，已被广泛应用于换热器设计[23]及化工过程优化[24-26]中。多目标优化得到的Pareto解为一系列最优解的集合，因此需进一步选择一组最佳的解作为MPC最优权重参数。Wang等[27]比较了典型的10种用于从Pareto解中选取最优解的方法，表明灰色关联法(gray relational analysis，GRA)具有较好的选择性能，且该方法可直接选出最优解，而不需要进一步设定每个目标的权重值。因此，用GRA方法从多目标优化所得的Pareto解中选取最优解。
2.2 单层非线性规划整定方法(NLP)方程(7)所示的多目标优化问题在求解过程中不仅存在内层MPC子优化问题和外层权重参数整定优化问题之间的迭代，而且需要大量的子代数和种群数量作为评估目标方程的随机样本点以求得最优Pareto解，因此MOO权重参数整定方法通常计算复杂，耗时较长，不能快速求解权重参数整定优化问题。然而在方程(7)中，有约束的MPC子优化问题可转化为外层MPC权重参数整定优化问题的一阶KKT最优条件，从而将方程(7)所示的双层结构MPC权重参数整定优化问题转化为单层非线性规划问题，运用基于梯度的非线性规划方法进行求解，可极大地缩短计算时间。
将方程(6)所示的MPC子优化问题转化为标准二次规划形式：
| $ \begin{gathered} \begin{array}{l} \ \ \min \limits_{\Delta \boldsymbol{u}_{k}} J_{k}=\frac{1}{2} \Delta \boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{H} \Delta \boldsymbol{u}_{k}+\boldsymbol{g}^{\mathrm{T}} \Delta \boldsymbol{u}_{k}+\text { Cons, }\\ {\text{subject to:}} \end{array}\\ \boldsymbol{A}_{\mathrm{b}} \Delta \boldsymbol{u}_{k} \leqslant \boldsymbol{b} \text { 。} \end{gathered} $ | (8) |
式中：
| $ L\left(\Delta \boldsymbol{u}_{k}, \boldsymbol{\lambda}, \boldsymbol{s}\right)=\frac{1}{2} \Delta \boldsymbol{u}_{k}^{\mathrm{T}} \boldsymbol{H} \Delta \boldsymbol{u}_{k}+\boldsymbol{g}^{\mathrm{T}} \Delta \boldsymbol{u}_{k}+\sum\limits_{i=1}^{n_{\rm c}} \lambda_{i}\left(A_{b i} \Delta \boldsymbol{u}_{k}-b_{i}+s_{i}\right) 。$ | (9) |
式中：λi和si分别为第i个线性约束的拉格朗日乘子和松弛因子；nc为MPC子优化问题线性约束不等式的个数；Abi和bi分别为第i个线性约束表达式的系数和约束值。故方程(8)的最优KKT条件可表示为：
| $ \left\{\begin{array}{l} \boldsymbol{H} \Delta \boldsymbol{u}_{k}+\boldsymbol{g}^{\mathrm{T}} \Delta \boldsymbol{u}_{k}+\boldsymbol{A}_{\mathrm{b}}^{\mathrm{T}} \boldsymbol{\lambda}=0 ; \\ \boldsymbol{A}_{\mathrm{b}} \Delta \boldsymbol{u}_{k}-\boldsymbol{b}+\boldsymbol{s}=0 ; \\ \lambda_{i} s_{i}=0, \quad \lambda_{i} \geqslant 0, s_{i} \geqslant 0, i=1,2, \cdots, n_{\mathrm{c} } 。\end{array}\right. $ | (10) |
方程(7)所示的优化目标Ψ1和Ψ2可通过权重化集成为单目标，并将方程(10)作为MPC权重参数整定优化问题的约束条件，可得方程(11)所示的MPC权重参数整定单层非线性规划法。
| $ \begin{gathered} \begin{array}{l} \min \limits_{Q, \boldsymbol{R}, \Delta \boldsymbol{U}_{1}, \cdots, \Delta \boldsymbol{u}_{N_{\rm t}}, \lambda_{1}, \cdots, \lambda_{N_{\rm t}}, s_{1}, \cdots, s_{N_{\rm t}}} \phi=\gamma \varPsi_{1}+(1-\gamma) \varPsi_{2}, \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\text{subject to:}} \end{array}\\ \boldsymbol{Q}_{\min } \leqslant \boldsymbol{Q} \leqslant \boldsymbol{Q}_{\max }, \\ \boldsymbol{R}_{\min } \leqslant \boldsymbol{R} \leqslant \boldsymbol{R}_{\max },\\ \left.\begin{array}{l} \boldsymbol{x}_{j}^{\mathrm{m}}=\boldsymbol{A}_{0} \boldsymbol{x}_{j-1}^{\mathrm{m}}+\boldsymbol{B}_{0} \boldsymbol{u}_{j-1}+\boldsymbol{B}_{\mathrm{d} 0} \boldsymbol{d}_{j-1}, \\ \boldsymbol{y}_{j}^{\mathrm{m}}=\boldsymbol{C}_{0} \boldsymbol{x}_{j}^{\mathrm{m}}, \\ \boldsymbol{u}_{0, j}=\boldsymbol{u}_{0, j-1}+\Delta \boldsymbol{u}_{0, j}, \\ \Delta \boldsymbol{u}_{0, j}=\boldsymbol{M} \Delta \boldsymbol{u}_{j}, \\ \boldsymbol{H} \Delta \boldsymbol{u}_{j}+\boldsymbol{g}^{\mathrm{T}} \Delta \boldsymbol{u}_{j}-\boldsymbol{A}_{\mathrm{b}} \lambda_{j}=0, \\ \boldsymbol{A}_{\mathrm{b}} \Delta \boldsymbol{u}_{j}-\boldsymbol{b}+\boldsymbol{s}_{j}=0, \\ \lambda_{i, j} s_{i, j}=0, \quad \lambda_{i, j} \geqslant 0, s_{i, j} \geqslant 0, i=1,2, \cdots, n_{\mathrm{c}}, \end{array}\right\} j=1,2, \cdots, N_{\mathrm{t} }。\end{gathered} $ | (11) |
式中：γ∈[0, 1]为目标函数的权重因子，其值越大，优化所得的权重参数使得MPC能够更快地跟踪参考轨迹设定值；Δu0, j为控制动作序列Δuj中的第一个控制动作，由转换矩阵M得到。
| $ \boldsymbol{M}=\operatorname{diag}\left(\left[\begin{array}{cc} \overbrace{1, \cdots, 1}^{m} & \overbrace{0, \cdots, 0}^{m \times\left(N_{\rm c}-1\right)} \end{array}\right]\right)。$ |
在Matlab 2019a中编译MPC控制算法，并通过Matlab中的Spreadlink工具箱实现Matlab与在MS Excel中建立的EMOO多目标优化求解器之间的数据传输。在Matlab 2019a中应用Casadi算法网执行非线性规划权重参数整定算法，并用ipopt求解器求解方程(11)所示的非线性规划问题。仿真案例的测试在Simulink 2019a中实施。本节分别以一个单输入单输出(SISO)二阶时间延时传递函数和连续搅拌反应器(CSTR)模型阐述所提出的MPC权重参数整定算法的有效性。
3.1 单输入单输出模型方程(12)给出了一个SISO二阶时间延时传递函数模型[12]。
| $ G(s)=\frac{0.164\ 9 s+1}{19\ s^{2}+6.5 s+1} \mathrm{e}^{(-2 s)}。$ | (12) |
式中：G为被控对象模型；s为复变量。
该传递函数为无约束系统，首先以采样时间ts=0.5 min离散化为方程(1)所示的状态空间模型，并取Np和Nc分别为60和3个采样点，故对于MPC控制器还需确定权重参数Q和R。为了简化计算，对于SISO系统，取输出权重Q=1，则优化过程中只需确定相应的R即可。故可通过优化整定时域内MPC控制器对参考轨迹的跟踪性能而获得最优权重参数值，选取整定时域为100采样点，起始点u0=0，y0=0，ysp=0；当t=2.5 min时，设定ysp=2；当t=25.0 min时，设定ysp=-2。
在MOO整定方法中，将遗传算法的代数和种群数分别设定为100和50，而在NLP整定方法中将权重因子γ分别设定为1.00、0.75和0.50以考察γ对MPC整定结果的影响。表 1为优化所得的R值。由于应用GRA方法从Pareto解中选取的最优解更倾向于使积分跟踪均方差Ψ1最小化，因此MOO整定方法所得的最优R值与γ=1.00时NLP整定方法所得的最优解相近。
|
表 1 SISO系统整定参数 Table 1 Tuning parameters for SISO system |
由图 1描述的以不同整定方法所得的最优权重R作为MPC权重参数时SISO系统的动态响应可以看出，随着权重因子γ逐渐增大，MPC趋于更快速地跟踪参考轨迹设定点(图 1(a))，但控制动作的变化幅度也同时增大；反之，MPC输出响应更趋于平缓。因此，在NLP整定方法中，可通过调节γ的大小优化MPC动态响应特性。当γ=1.00时，NLP和MOO两种整定方法所得的MPC动态响应性能相近。
|
图 1 SISO系统设定点跟踪MPC动态响应 Fig. 1 Dynamic responses of MPC controller for SISO system under setpoint changes |
不同整定方法所得的MPC的动态响应性能由方程(13)所示的积分绝对偏差(integral absolute error，IAE)予以定量评估。
| $ \mathrm{IAE}=\int_{0}^{t}\left|\boldsymbol{y}_{\tau}^{\mathrm{m}}-\boldsymbol{y}_{\tau}^{\mathrm{sp}}\right| \mathrm{d} \tau 。$ | (13) |
式中：t为动态模拟时间长度；yτm和yτsp分别表示τ时刻的工艺输出值和控制器设定点。
表 2给出了应用不同整定方法所得的最优权重R时MPC输出动态响应的IAE值。可以看出，应用NLP整定方法时随着γ值减小，MPC输出响应的IAE逐渐变大；而当γ=1.00时，应用NLP和MOO整定方法所得的MPC的输出动态响应的IAE值非常接近。NLP整定方法的计算用时仅需5~10 s，而MOO整定方法则需约1 h，由此可见，NLP整定方法能够快速地整定MPC权重参数，且整定性能与MOO整定方法相近。
|
|
表 2 SISO系统输出动态响应IAE值 Table 2 IAE values of dynamic responses of the output variable for SISO system |
一连续搅拌反应器如图 2所示，反应物A以体积流量F、摩尔浓度CA0、温度T0进入反应器内，通过液相二阶可逆反应转化为产物B，反应产生的热量经夹套冷剂以热流速率QR移走。反应器可描述为方程(14)和(15)所示的数学模型，其参数值见表 3[28]。
| $ \frac{\mathrm{d} C_{\mathrm{A}}}{\mathrm{d} t}=\frac{F}{V}\left(C_{\mathrm{A} 0}-C_{\mathrm{A}}\right)-k_{0} \mathrm{e}^{-E / R_{\mathrm{g}} T} C_{\mathrm{A}}^{2} ; $ | (14) |
| $ \frac{\mathrm{d} T}{\mathrm{~d} t}=\frac{F}{V}\left(T_{0}-T\right)-\frac{\Delta H k_{0}}{\rho_{\mathrm{L}} C_{\mathrm{p}}} \mathrm{e}^{-E / R_{\mathrm{g}} T} C_{\mathrm{A}}^{2}+\frac{Q_{\mathrm{R}}}{\rho_{\mathrm{L}} C_{\mathrm{p}} V}。$ | (15) |
|
图 2 连续搅拌反应器示意图 Fig. 2 Schematic diagram of CSTR |
|
|
表 3 CSTR模型参数 Table 3 Model parameters of CSTR |
式中：V为反应器体积；T为反应器内液体温度；k0为反应速率常数；E为反应活化能；Cp为反应器内液体比热容；ρL为反应器内液体密度；ΔH为反应热；Rg为气体常数。
反应物A的温度T0设为可测扰动量，因此，反应器的操作目标为通过调节反应物摩尔浓度CA0和热流速率QR使反应器内反应物浓度CA和温度T维持在设定点。故反应器操作过程中，其控制变量、操作变量及其可测扰动可表示为方程(16)~(18)，其中下标s表示每个变量的稳态值(见表 3)。
| $ \boldsymbol{y}=\left[y_{1}, y_{2}\right]^{\mathrm{T}}=\left[C_{\mathrm{A}}-C_{\mathrm{As}}, T-T_{\mathrm{s}}\right]^{\mathrm{T}} ; $ | (16) |
| $ \boldsymbol{u}=\left[u_{1}, u_{2}\right]^{\mathrm{T}}=\left[C_{\mathrm{A} 0}-C_{\mathrm{A} 0 \mathrm{~s}}, Q_{\mathrm{R}}-Q_{\mathrm{Rs}}\right]^{\mathrm{T}} ; $ | (17) |
| $ d=\left[T_{0}-T_{0 \mathrm{~s}}\right]。$ | (18) |
操作变量、控制变量和可测扰动的范围见表 4。
|
|
表 4 变量范围 Table 4 Boundary conditions of variables |
方程(14)和(15)在稳态操作点经线性化并以采样时间为0.02 h离散化后可得方程(1)所示的线性动态模型，其参数分别为：
| $ \boldsymbol{A}_{0}=\left[\begin{array}{cccc} 0.498\ 3 & -0.007\ 8 \\ 2.240\ 1 & 1.292\ 4 \end{array}\right], \ \ \ \ \boldsymbol{B}_{0}=\left[\begin{array}{ccc} 0.074\ 1 & -3.497\ 9 \times 10^{-7} \\ 2.240\ 1 & 9.980\ 6 \times 10^{-5} \end{array}\right], $ |
| $ \boldsymbol{B}_{\mathrm{d} 0}=\left[\begin{array}{c} -4.040\ 1 \times 10^{-4} \\ 0.115\ 3 \end{array}\right], \ \ \ \ \boldsymbol{C}_{0}=\left[\begin{array}{cc} -1 & 0 \\ 0 & -1 \end{array}\right] \text { 。} $ |
为了方便权重参数整定优化问题的求解，模型参数需根据变量各自的范围区间转化为无维度值，并取Np和Nc分别为10和3个采样点，整定时域选取100个采样点。在MOO整定方法中，遗传算法的代数和种群数分别设定为100和50，NLP整定方法中γ分别设定为1.00、0.75和0.50，以考察γ对于MPC整定结果的影响。该系统为抗扰动干扰的MPC权重参数整定，因此，在整定过程中设定，当t=0.02 h时，可测扰动原料温度升高5 K，即方程(1)中d=5(转化后的无维度值)。
表 5给出了优化所得的权重参数Q和R的值。图 3描述了原料温度T0升高5 K时，MPC控制器输入与输出变量动态响应。可以看出，NLP和MOO两种整定方法所得的MPC控制器都具有很好的抗可测温度扰动的能力，能够快速平滑地回到设定点。表 6给出了原料温度T0升高5 K后CA和T的动态响应IAE值。可以看出，虽然MOO整定方法所得的MPC控制器CA的IAE值较NLP整定方法所得的IAE值小，但MOO整定方法所得的MPC控制器T的IAE值较γ=1.00时NLP整定方法所得的IAE值大。此外，表 6中γ=1.00时NLP整定方法所得的MPC输出响应IAE值总和远较MOO整定方法所得的IAE值总和小。因此，NLP整定方法所得的MPC控制性能更好，或与MOO整定方法所得的MPC控制器性能相近。
|
|
表 5 CSTR系统整定参数 Table 5 Tuning parameters for CSTR system |
|
图 3 CSTR系统T0升高5 K时MPC动态响应 Fig. 3 Dynamic responses of MPC controller for CSTR system under 5 K increase of T0 |
|
|
表 6 原料温度升高5 K时CSTR系统输出响应IAE值 Table 6 IAE values of dynamic responses of the output variables for CSTR system under 5 K increase of feed temperature |
图 4描述了MPC控制器对反应器浓度CA的设定点改变的跟踪性能。此处假设当t=0.5 h时，CA的设定点由初始稳态值1.22 kmol · m-3变为1.35 kmol · m-3；当t=4.0 h时，CA的设定点由1.35 kmol · m-3变为1.15 kmol · m-3。可以看出随着γ值减小，NLP整定方法所得的MPC趋于缓慢、平滑地跟踪CA的设定点，而T的瞬时偏差则更小。MOO整定方法所得的MPC虽然能快速地跟踪CA的设定点，但T的瞬时偏差较γ=1.00的NLP整定方法所得的MPC大。表 7给出了CA设定点改变后CA和T动态响应的IAE值。与原料温度变化时MPC的动态响应性能相似，随着γ值减小，NLP整定方法所得的MPC控制器CA和T的IAE值都逐渐变大；当γ=1.00时NLP整定方法所得的MPC输出动态响应IAE值与MOO整定方法的MPC的IAE相近，NLP整定方法的总IAE值略小于MOO整定方法的总IAE值。
|
图 4 CSTR系统CA设定点跟踪MPC动态响应 Fig. 4 Dynamic responses of MPC controller for CSTR system under the changes of CA setpoint |
|
|
表 7 CA设定点跟踪CSTR系统输出动态响应IAE值 Table 7 IAE values of dynamic responses of the output variables for CSTR system under the changes of CA setpoint |
综上所述，NLP整定方法所得的MPC控制性能优于或接近于MOO整定方法所得的MPC，且在整定过程中，NLP整定方法求解优化问题所需的计算时间仅为10~90 s，而MOO整定方法所需的计算时间则约为1.5 h。由图 5所示的MOO整定方法所得的Pareto解可以看到GRA方法选取的最优点(蓝色五角星)趋向于最小化参考轨迹跟踪目标Ψ1，这与方程(11)中γ=1.00时的优化意义一致，从另一层面表明了两种整定方法的一致性。
|
图 5 MPC权重参数整定多目标优化Pareto解 Fig. 5 Pareto solution of multi-objective optimization for MPC weight parameter tuning |
权重参数的整定是MPC控制器取得良好控制性能的重要因素。针对基于双层结构多目标优化的MPC权重参数整定方法存在求解过程较慢、耗时较长的问题，建立了一种非线性规划的权重参数整定方法。该方法将MPC权重参数整定中每个时间采样点的MPC子优化问题等价为外层MPC权重参数整定优化问题的最优KKT(Karush-Kuhn-Tucker)条件，将MPC权重参数整定的双层多目标优化问题转化为单层非线性规划问题，并应用一个SISO时间延时传递函数和CSTR反应器仿真案例评估该算法的有效性。结果表明：
1) 基于单层结构的非线性规划整定方法能够快速实现MPC权重参数的整定，极大地降低权重参数整定优化问题的求解时间。对于SISO系统和CSTR系统，非线性规划权重参数整定优化问题的求解时间仅为5~10 s和10~90 s，而多目标优化的整定方法则分别为1.0 h和1.5 h。
2) 基于单层非线性规划整定方法的MPC控制器控制性能好于或接近基于多目标优化整定方法的MPC控制器控制性能。
| [1] |
Qin S J, Badgwell T A. A survey of industrial model predictive control technology[J]. Control Engineering Practice, 2003, 11(7): 733-764. DOI:10.1016/S0967-0661(02)00186-7 |
| [2] |
罗雄麟, 于洋, 许鋆. 化工过程预测控制的在线优化实现机制[J]. 化工学报, 2014, 65(10): 3984-3992. Luo X L, Yu Y, Xu J. Online optimization implementation on model predictive control in chemical process[J]. CIESC Journal, 2014, 65(10): 3984-3992. (in Chinese) DOI:10.3969/j.issn.0438-1157.2014.10.032 |
| [3] |
Garriga J L, Soroush M. Model predictive control tuning methods: a review[J]. Industrial & Engineering Chemistry Research, 2010, 49(8): 3505-3515. |
| [4] |
Trierweiler J O, Farina L A. RPN tuning strategy for model predictive control[J]. Journal of Process Control, 2003, 13(7): 591-598. DOI:10.1016/S0959-1524(02)00093-8 |
| [5] |
Hinde R F Jr, Cooper D J Jr. A pattern-based approach to excitation diagnostics for adaptive process control[J]. Chemical Engineering Science, 1994, 49(9): 1403-1415. DOI:10.1016/0009-2509(94)85069-0 |
| [6] |
Shridhar R, Cooper D J. A tuning strategy for unconstrained SISO model predictive control[J]. Industrial & Engineering Chemistry Research, 1997, 36(3): 729-746. |
| [7] |
Shridhar R, Cooper D J. A tuning strategy for unconstrained multivariable model predictive control[J]. Industrial & Engineering Chemistry Research, 1998, 37(10): 4003-4016. |
| [8] |
Gous G Z, de Vaal P L. Using MV overshoot as a tuning metric in choosing DMC move suppression values[J]. ISA Transactions, 2012, 51(5): 657-664. DOI:10.1016/j.isatra.2012.05.006 |
| [9] |
Davtyan A, Hoffmann S, Scheuring R. Optimization of model predictive control by means of sequential parameter optimization[C]//Computational Intelligence in Control and Automation (CICA), April 11-15, 2011, Paris, France. IEEE, 2011: 11-16.
|
| [10] |
Exadaktylos V, Taylor C J. Multi-objective performance optimisation for model predictive control by goal attainment[J]. International Journal of Control, 2010, 83(7): 1374-1386. DOI:10.1080/00207171003736295 |
| [11] |
Nery G A Jr, Martins M A F, Kalid R. A PSO-based optimal tuning strategy for constrained multivariable predictive controllers with model uncertainty[J]. ISA Transactions, 2014, 53(2): 560-567. DOI:10.1016/j.isatra.2013.12.019 |
| [12] |
Yamashita A S, Alexandre P M, Zanin A C, et al. Reference trajectory tuning of model predictive control[J]. Control Engineering Practice, 2016, 50: 1-11. DOI:10.1016/j.conengprac.2016.02.003 |
| [13] |
Yamashita A S, Zanin A C, Odloak D. Tuning the model predictive control of a crude distillation unit[J]. ISA Transactions, 2016, 60: 178-190. DOI:10.1016/j.isatra.2015.10.017 |
| [14] |
Giraldo S A C, Melo P A, Secchi A R. Tuning of model predictive control based on hybrid optimization[J]. IFAC-PapersOnLine, 2019, 52(1): 136-141. DOI:10.1016/j.ifacol.2019.06.050 |
| [15] |
Gutiérrez-Urquídez R C, Valencia-Palomo G, Rodríguez-Elias O M, et al. Systematic selection of tuning parameters for efficient predictive controllers using a multiobjective evolutionary algorithm[J]. Applied Soft Computing, 2015, 31: 326-338. DOI:10.1016/j.asoc.2015.02.033 |
| [16] |
Lozano Santamaría F, Gómez J M. An algorithm for tuning NMPC controllers with application to chemical processes[J]. Industrial & Engineering Chemistry Research, 2016, 55(34): 9215-9228. |
| [17] |
Vallerio M, van Impe J, Logist F. Tuning of NMPC controllers via multi-objective optimisation[J]. Computers & Chemical Engineering, 2014, 61: 38-50. |
| [18] |
Tran Q N, Özkan L, Backx A C P M. Generalized predictive control tuning by controller matching[J]. Journal of Process Control, 2015, 25: 1-18. DOI:10.1016/j.jprocont.2014.10.002 |
| [19] |
Feng Z M, Shen W F, Rangaiah G P, et al. Proportional-integral control and model predictive control of extractive dividing-wall column based on temperature differences[J]. Industrial & Engineering Chemistry Research, 2018, 57(31): 10572-10590. |
| [20] |
Feng Z M, Shen W F, Rangaiah G P, et al. Closed-loop identification and model predictive control of extractive dividing-wall column[J]. Chemical Engineering and Processing -Process Intensification, 2019, 142: 107552. DOI:10.1016/j.cep.2019.107552 |
| [21] |
Sharma S, Rangaiah G P, Cheah K S. Multi-objective optimization using MS Excel with an application to design of a falling-film evaporator system[J]. Food and Bioproducts Processing, 2012, 90(2): 123-134. DOI:10.1016/j.fbp.2011.02.005 |
| [22] |
Deb K, Pratap A, Agarwal S, et al. A fast and elitist multiobjective genetic algorithm: NSGA-Ⅱ[J]. IEEE Transactions on Evolutionary Computation, 2002, 6(2): 182-197. DOI:10.1109/4235.996017 |
| [23] |
Wong J Y Q, Sharma S, Rangaiah G P. Design of shell-and-tube heat exchangers for multiple objectives using elitist non-dominated sorting genetic algorithm with termination criteria[J]. Applied Thermal Engineering, 2016, 93: 888-899. DOI:10.1016/j.applthermaleng.2015.10.055 |
| [24] |
Oni A O, Fadare D A, Sharma S, et al. Multi-objective optimisation of a double contact double absorption sulphuric acid plant for cleaner operation[J]. Journal of Cleaner Production, 2018, 181: 652-662. DOI:10.1016/j.jclepro.2018.01.239 |
| [25] |
da Cunha S, Rangaiah G P, Hidajat K. Design, optimization, and retrofit of the formic acid process I: base case design and dividing-wall column retrofit[J]. Industrial & Engineering Chemistry Research, 2018, 57(29): 9554-9570. |
| [26] |
da Cunha S, Rangaiah G P, Hidajat K. Design, optimization, and retrofit of the formic acid process Ⅱ: reactive distillation and reactive dividing-wall column retrofits[J]. Industrial & Engineering Chemistry Research, 2018, 57(43): 14665-14679. |
| [27] |
Wang Z Y, Rangaiah G P. Application and analysis of methods for selecting an optimal solution from the Pareto-optimal front obtained by multiobjective optimization[J]. Industrial & Engineering Chemistry Research, 2017, 56(2): 560-574. |
| [28] |
Alanqar A, Durand H, Christofides P D. On identification of well-conditioned nonlinear systems: application to economic model predictive control of nonlinear processes[J]. AIChE Journal, 2015, 61(10): 3353-3373. DOI:10.1002/aic.14942 |