-
摘要: 主成分分析(Principal component analysis, PCA) 是一种广泛应用的特征提取与数据降维方法, 其目标函数采用L2范数距离度量方式, 对离群数据及噪声敏感. 而L1范数虽然能抑制离群数据的影响, 但其重构误差并不能得到有效控制. 针对上述问题, 综合考虑投影距离最大及重构误差较小的目标优化问题, 提出一种广义余弦模型的目标函数. 通过极大化矩阵行向量的投影距离与其可调幂的2范数之间的比值, 使得其在数据降维的同时提高了鲁棒性. 在此基础上提出广义余弦二维主成分分析(Generalized cosine two dimensional PCA, GC2DPCA), 给出了其迭代贪婪的求解算法, 并对其收敛性及正交性进行理论证明. 通过选择不同的可调幂参数, GC2DPCA可应用于广泛的含离群数据的鲁棒降维. 人工数据集及多个人脸数据集的实验结果表明, 本文算法在重构误差、相关性及分类率等性能方面均得到了提升, 具有较强的抗噪能力.Abstract: Principal component analysis (PCA) is a widely used feature extraction and dimensionality reduction method, which employs L2-norm distance metric and is sensitive to outliers and noise. Although L1-norm distance metric can suppress outliers, the reconstruction error in the objective function cannot be effectively controlled. To solve these problems, a new objective function based on the generalized cosine model is first proposed by considering the maximum projection distance and the smaller reconstruction error. By maximizing the ratio between the projection distance of row vectors and the 2-norm with the adjustable power, the robustness is improved while reducing dimensions. Then, the generalized cosine two-dimensional PCA (GC2DPCA) is presented by using the model and can be solved by an iterative greedy algorithm. The convergence and orthogonality of the algorithm can be proved theoretically. By choosing a different adjustable power parameter, GC2DPCA can be widely used in robust dimensionality reduction with outliers. The experimental results on artificial synthetic datasets and several face datasets show that the proposed algorithm not only improves the performances of reconstruction error, orthogonal correlation and classification rate, but also has strong anti-noise ability.
-
化工过程普遍存在不确定性, 如何采用有效的优化方法找到不确定条件下系统的真实最优点, 对提高化工企业的经济效益发挥着关键作用[1]. 大规模化工过程的控制系统通常为分层递阶结构[2-3], 控制层(下层)的主要任务是抑制底层扰动, 跟踪优化层传递来的被控变量设定值, 优化层(上层)根据调度层(顶层)传达的生产任务指标等, 对当前工况进行识别, 执行相应的优化算法计算出最优设定值, 传递给控制层执行.
优化层执行的优化算法通常以化工过程的非线性模型为基础, 以传统的“二步法” 实时优化[4-5] 为例, 首先确定模型的不确定参数, 运行过程中采集系统的输出量数据, 对未知参数进行估计, 再基于更新的系统模型进行重优化, 计算出被控变量的最优设定值后传递给控制层. 这一过程通常还需要结合数据调和、稳态检测等技术手段加强优化结果的可靠性, 工业过程的优化周期一般为4 ~ 8小时. 针对传统的“二步法” 的缺点, 近年来涌现出了新的实时优化方法, 如Bonvin课题组提出的修正项自适应方法(Modifier adaptation)[6-7], 通过对标称模型的约束及梯度进行修正, 即使不估计扰动参数也能收敛到真实最优点. 文献[8-10]考虑运行层之间的不同时间尺度, 提出了数据驱动的多速率分层运行优化控制方法, 基于Q学习对基础控制回路的设定值进行在线优化, 使运行层能更好地优化控制性能指标. 自优化控制(Self-optimizing control, SOC)[11-13] 提出通过离线选择控制层的被控变量, 设定值则在线保持不变, 提供了实时优化的另一种研究思路. 在自优化控制中, 被控变量可以是常规物理量的函数, 即构造虚拟量进行控制, 可使系统的操作变量可以在不确定性下进行自寻优. 当底层控制的优化作用较强时(经济损失可接受), 甚至可以省略单独的优化层, 从而简化控制系统. 相比传统的优化方式, 自优化控制的优化在工作频率为秒/分的反馈控制中完成, 因此优化速度得到大幅度提升, 在一系列研究中表现出良好的效果[14-17].
间歇过程是一类批次加工的化工过程, 具有规模小、灵活性高的特点, 在需求多元化的现代市场中应用越来越广泛. 相比连续化工过程, 间歇过程具有“多重时变” 的操作特征[18-19]. 一方面, 间歇过程具有重复特性, 可以引入学习机制从历史批次的数据中提炼出有用的信息, 改进后续批次的跟踪控制和经济指标优化, 典型的如迭代学习控制[18, 20-23]、批间实时优化[23-24] 等控制和优化技术. 另一方面, 由于其时变特性, 间歇过程在批次内无稳定操作点, 相比连续过程的控制和稳态优化更具挑战[25-27]. 自优化控制经过近20年的发展, 针对连续过程已报道了一系列被控变量求解方法[12, 28-31], 但是针对需动态优化的间歇过程仍缺乏足够的研究. 值得注意的是, 由于从批间角度看间歇过程是一个静态过程[32], 近年来文献[33-34]提出了间歇过程的批间自优化控制方法. 此类方法仅利用了间歇过程的重复性, 基于已有的静态自优化控制方法求解被控变量, 然后设计批间控制器调节输入轨迹, 逐批次将被控变量控制于恒定设定点, 实现实时优化. 但批间优化本质上还是静态方法, 由于需要若干个批次才能实现被控变量的跟踪控制, 优化作用慢, 因此未充分发挥自优化控制的优势. 此外, 批间优化只对具有重复特性的扰动具有效果, 当系统受到高频扰动作用时, 批间控制器难以实现有效的实时优化.
最近, Ye等[35]提出了一种针对间歇过程的动态自优化控制方法, 通过考虑批内变量的因果性, 最终得到了具有优化作用的控制律. 设计控制系统时, 选择被控变量和设计控制器通常是两个独立任务[36], 前者主要考虑经济指标的优化, 后者关注于如何更好地跟踪控制被控变量, 保证控制系统的稳定性和鲁棒性. 如何在此前提下求解批内被控变量, 仍是一个开放的课题.
本文研究了间歇过程的批内自优化控制问题, 贡献如下: 1)基于自优化控制策略提出以输出变量的线性组合为被控变量(虚拟变量), 在批次运行过程中对其进行跟踪控制, 以控制手段实现实时优化; 2)根据是否在过程不同阶段切换被控变量, 给出了两种自优化控制策略, 对每种策略又分别给出了两种设定轨线选取方案; 3)引入扩张组合矩阵, 将这些情形统一描述为具有不同结构约束的最优组合矩阵求解问题, 并推导得到了其中一种方案的解析解计算方法. 目前为止, 本文所提方法在国内外文献中未见报道.
1. 连续过程的自优化控制
对连续化工过程, 考虑如下静态优化问题
$$ \begin{split} &\min\limits_{{\boldsymbol{u}}} J({\boldsymbol{u}}, {\boldsymbol{d}})\\ &\;{\rm{s}}.{\rm{t}}. \quad {\boldsymbol{y}} = {\boldsymbol{g}}({\boldsymbol{u}}, {\boldsymbol{d}})\\& \quad\quad\;\; {\boldsymbol{g}}_{\rm{in}}({\boldsymbol{u}}, {\boldsymbol{d}}) \leq0 \end{split} $$ (1) 其中,
$ J $ 为经济指标,$ {\boldsymbol{u}}\in {\bf R}^{n_u} $ ,$ {\boldsymbol{d}}\in {\bf R}^{n_d} $ 和$ {\boldsymbol{y}}\in {\bf R}^{n_y} $ 分别是操纵变量、扰动变量和测量变量,$ {\boldsymbol{g}} $ 和$ {\boldsymbol{g}}_{\rm{in}} $ 为输出变量的模型函数和约束条件.扰动变量
$ {\boldsymbol{d}} $ 变化且在线不可测是化工过程偏离最优点的主要原因. 当扰动变量$ {\boldsymbol{d}} $ 变化时, 式(1)的解是$ {\boldsymbol{d}} $ 的函数, 不妨记为$ {\boldsymbol{u}}^{\rm{opt}}({\boldsymbol{d}}) $ . 实时优化的任务是在$ {\boldsymbol{d}} $ 未知的前提下, 寻找到新的最优值$ {\boldsymbol{u}}^{\rm{opt}} $ , 实现过程的最优操作. 自优化控制(SOC)通过构造虚拟的被控变量$ {\boldsymbol{c}} = H{\boldsymbol{y}} $ , 当反馈控制器将$ {\boldsymbol{c}} $ 控制在恒定设定值$ {\boldsymbol{c}}_s $ 上时, 控制器输出能自动逼近当前的实际最优值$ {\boldsymbol{u}}^{\rm{opt}}({\boldsymbol{d}}) $ . 当组合矩阵$ H $ 每行有且只有一个1, 其余为0时,$ {\boldsymbol{c}} $ 为输出变量$ {\boldsymbol{y}} $ 的子集, 此时退化为传统的以单变量为被控变量的情形. 更一般的情况下,$ H $ 中的非零元素提供了更多优化自由度, 可提高系统的闭环经济性能. 例如, 假设系统自由度$ n_u = 2 $ ,$ {\boldsymbol{y}} = [T\;P\;c_A]^{\rm{T}} $ , 包括温度$ T $ , 压力$ P $ 和物质A的浓度$ c_A $ , 考虑两种情况:$$ {H_1} = \left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&0&1 \end{array}} \right];\;{H_2} = \left[ {\begin{array}{*{20}{c}} {{h_{11}}}&{{h_{12}}}&{{h_{12}}}\\ {{h_{21}}}&{{h_{22}}}&{{h_{23}}} \end{array}} \right] $$ $ H_1 $ 对应的被控变量$ {\boldsymbol{c}} = H_1{\boldsymbol{y}} $ 为$ T $ 和$ c_A $ (单个物理量),$ H_2 $ 的被控变量为3个物理量的线性组合. 显然, 前者为后者的一种特殊形式.为求解一般情形的最优组合矩阵
$ H $ , 研究人员针对不同过程特性和衡量标准提出了求解方法[12, 28-31]. 以一种针对线性系统的局部法(Exact local method)为例[28], 首先定义损失函数$ L $ $$ L = J({\boldsymbol{u}}, {\boldsymbol{d}})-J^{\rm{opt}}({\boldsymbol{d}}) $$ (2) 对给定的
$ {\boldsymbol{d}} $ , 将$ J({\boldsymbol{u}}, {\boldsymbol{d}}) $ 在最优点$ {\boldsymbol{u}}^{\rm{opt}} $ 处进行二阶泰勒展开$$ \begin{split} J({\boldsymbol{u}}, {\boldsymbol{d}}) \approx\;& J^{\rm{opt}}({\boldsymbol{d}})+J_{u}^{{\rm{T}}}\left({\boldsymbol{u}}-{\boldsymbol{u}}^{\rm{opt}}\right)+\\&\frac{1}{2}\left({\boldsymbol{u}}-{\boldsymbol{u}}^{\rm{opt}}\right)^{{\rm{T}}} J_{u u}\left({\boldsymbol{u}}-{\boldsymbol{u}}^{\rm{opt}}\right) \end{split} $$ (3) 式中,
$J_u = \frac{\partial J}{\partial {\boldsymbol{u}}}$ 和$J_{uu} =\frac{ \partial^2J}{\partial {\boldsymbol{u}}^2 }$ 分别为一阶和二阶敏感矩阵. 根据最优性条件, 最优点处$ J_u = 0 $ , 将其代入式(1)可得到二次型损失函数$$ L = \frac{1}{2}\left({\boldsymbol{u}}-{\boldsymbol{u}}^{\rm{opt}}\right)^{{\rm{T}}} J_{u u}\left({\boldsymbol{u}}-{\boldsymbol{u}}^{\rm{opt}}\right) $$ (4) 此外, 输出函数在标称点处的线性化方程为
$$ {\boldsymbol{y}} = G_{y} {\boldsymbol{u}}+G_{yd}{\boldsymbol{d}} $$ (5) 若考虑测量变量含噪声:
$ {\boldsymbol{y}}_m = {\boldsymbol{y}}+{\boldsymbol{n}} $ , 当控制器将被控变量$ {\boldsymbol{c}} = H{\boldsymbol{y}}_m $ 控制在原设定值上时$ (\Delta {\boldsymbol{c}} = 0) $ , 扰动变化$ \Delta{\boldsymbol{d}} $ 将引起的$ {\boldsymbol{u}} $ 变化量为$$ \Delta{\boldsymbol{u}} = -\left(HG_{y}\right)^{-1} HG_{yd} \Delta{\boldsymbol{d}}+\left(HG_{y}\right)^{-1} H{\boldsymbol{n}} $$ (6) 同时, 扰动变化
$ \Delta{\boldsymbol{d}} $ 将引起最优点变化$$ \Delta {\boldsymbol{u}}^{\rm{opt}} = -J_{uu}^{-1}J_{ud} \Delta{\boldsymbol{d}} $$ (7) 式中,
$J_{ud} = \frac{\partial^2J}{\partial {\boldsymbol{u}}\partial{\boldsymbol{d}} }$ . 将式(6)和式(7)代入式(4) 可计算损失函数[28]$$ \begin{split} &L = \frac{1}{2}\|{\boldsymbol{z}}\|^2_2\\ &{\boldsymbol{z}}= V(HG_y)^{-1}H\tilde F\begin{bmatrix}{{{\boldsymbol{d}}'}}\\ {\boldsymbol{n}}' \end{bmatrix} = M\begin{bmatrix}{{{\boldsymbol{d}}'}}\\ {\boldsymbol{n}}' \end{bmatrix} \end{split} $$ (8) 式中,
$ V $ 满足$ V^{\rm{T}}V = J_{uu} $ ,$ \tilde F = [FW_d\;\;\;W_n] $ ,$F = \frac{\partial {\boldsymbol{y}}^{\rm{opt}}}{ \partial {\boldsymbol{d}}} = -(G_yJ_{uu})^{-1}J_{ud}+G_{yd}$ ,$ W_d $ 和$ W_n $ 为对角矩阵, 其对角元素为$ {\boldsymbol{d}} $ 和$ {\boldsymbol{n}} $ 的幅值,$ {\boldsymbol{d}}' $ 和$ {\boldsymbol{n}}' $ 为归一化后的扰动和噪声向量, 其最大范数为1.式(8)中的损失
$ L $ 为单个工况$ (\Delta{\boldsymbol{d}},{\boldsymbol{n}}) $ 引起的损失. 当$ {\boldsymbol{d}} $ 和$ {\boldsymbol{n}} $ 在其变化范围内变化时, 文献[23-24]中分别以$ L $ 的最大值$ L_{\max} $ 和平均值$ L_{\rm{av}} $ 为衡量标准, 提出了$ H $ 的求解方法. 以平均损失$ L_{\rm{av}} $ 为例, 当$ {\boldsymbol{d}}' $ 和$ {\boldsymbol{n}}' $ 为正态分布时, 可得到如下最优化问题求解$ H $ $$ \min\limits_{H} L_{\rm{opt}} = 0.5\|M\|_{\rm{F}}^{2} $$ (9) 引理1[35].
$ L_{\rm{av}}(H) = L_{\rm{av}}(QH) $ , 其中$ Q $ 为任意$ n_u $ 维非奇异方阵.引理1表明, 式(9)的解非唯一(因为控制
$ {\boldsymbol{c}} = H{\boldsymbol{y}} $ 和$ {\boldsymbol{c}} = QH{\boldsymbol{y}} $ 等效). 利用该特性, 可以先求解出式(9)的一个特解, 再推广至通解形式. 文献[24]给出了最优$ H $ 的一个特解, 即$$ H^{\rm{T}} = (\tilde{F}\tilde{F}^{\rm{T}})^{-1}{G}_y({G}_y^{\rm{T}}(\tilde{F}\tilde{F}^{\rm{T}})^{-1}G_y)^{-1}J_{uu}^{\frac{1}{2}} $$ (10) 2. 间歇过程的自优化控制
2.1 间歇过程优化
考虑一类带不确定参数的间歇过程优化问题
$$ \begin{split} \min\limits_{{\boldsymbol{u}}(t)} J &= \phi\left({\boldsymbol{y}}\left(t_{f}\right)\right)\\ {\rm{s.t.}} \quad\; &\dot{{\boldsymbol{x}}} = {\boldsymbol{f}}({\boldsymbol{x}}, {\boldsymbol{u}}, {\boldsymbol{d}}), \quad {\boldsymbol{x}}(0) = {\boldsymbol{x}}_{0}\\ &{\boldsymbol{y}} = {\boldsymbol{f}}_{y}({\boldsymbol{x}},{\boldsymbol{u}}) \\ &{\boldsymbol{u}}_{L} \leq {\boldsymbol{u}}(t) \leq {\boldsymbol{u}}_{U} \\ &{\boldsymbol{T}}({\boldsymbol{x}}, {\boldsymbol{u}}) \leq 0 \end{split} $$ (11) 式中,
$ J $ 为最小化目标,$ {\boldsymbol{u}}(t)\in {\bf R}^{n_u} $ 为操纵变量轨迹($ {\boldsymbol{u}}_L $ 和$ {\boldsymbol{u}}_U $ 分别为输入下上限),$ {\boldsymbol{x}}\in {\bf R}^{n_x} $ ,$ {\boldsymbol{y}}\in {\bf R}^{n_y} $ 和$ {\boldsymbol{d}}\in {\bf R}^{n_d} $ 分别为状态向量(初态$ {\boldsymbol{x}}_0 )$ 、测量变量和不确定扰动,$ t_f $ 为批次运行时间,$ \phi $ 为目标函数,$ {\boldsymbol{f}} $ ,$ {\boldsymbol{f}}_y $ 和$ {\boldsymbol{T}} $ 分别为模型方程、输出方程和过程约束.对式(11)所示的动态优化问题, 通常可以基于数值法将其近似为离散化的非线性规划(Non-linear programming, NLP)问题[37]
$$ \begin{split} \min_{{\boldsymbol{u}}_i,\cdots,{\boldsymbol{u}}_N}& J = \phi\left({\boldsymbol{y}}(N)\right)\\ {\rm{s.t. }} \;\;\quad &{\boldsymbol{x}}(i+1) = {\boldsymbol{\hat f}}({\boldsymbol{x}}(i), {\boldsymbol{u}}(i), {\boldsymbol{d}}(i)), \quad {\boldsymbol{x}}(0) = {\boldsymbol{x}}_{0}\\ &{\boldsymbol{y}}(i) = {\boldsymbol{\hat f}}_{y}({\boldsymbol{x}}(i),{\boldsymbol{u}}(i)) \\ &{\boldsymbol{u}}_{L} \leq {\boldsymbol{u}}(i) \leq {\boldsymbol{u}}_{U} \\ &{\boldsymbol{T}}({\boldsymbol{x}}(i), {\boldsymbol{u}}(i)) \leq 0\\ &\forall\;i = 1,\cdots,N \\[-10pt]\end{split} $$ (12) 式中,
$ N $ 为间歇过程在操作区间$ [0, t_f] $ 内的离散段数,$ {\boldsymbol{\hat f}} $ 和$ {\boldsymbol{\hat f}}_y $ 代表离散后的非线性状态方程和输出方程.对上述间歇过程的优化问题, 文献[27-28]提出了批间自优化控制方法, 即构造被控变量
$ {\boldsymbol{c}} = H{\boldsymbol{y}} $ , 利用间歇过程的重复特性逐批次将$ {\boldsymbol{c}} $ 控制在恒设定值上. 从批间角度看, 间歇过程是一个静态过程, 因此第1节中针对连续过程的被控变量求解方法可以较为直接地拓展至间歇过程. 但批间优化需要经历若干批次实现被控变量的控制, 优化速度较慢. 并且, 若扰动的变化频率较高(如非重复性扰动), 则难以实现被控变量的跟踪控制, 优化效果有限.2.2 批内自优化控制策略
本文研究间歇过程的批内自优化控制方法, 即在单批次中控制被控变量实现实时优化. 与批间优化相比, 批内优化的响应速度更快, 能提高优化效果. 由于跟踪控制在单批次内完成, 批内优化能应对非重复性扰动. 对被控变量
$ {\boldsymbol{c}} = H{\boldsymbol{y}} $ 及其设定值$ {\boldsymbol{c}}_s $ , 考虑如下几种策略:策略1.
$ H $ 和$ {\boldsymbol{c}}_s $ 保持恒定;策略2.
$ H $ 恒定,$ {\boldsymbol{c}}_s $ 时变;策略3.
$ H $ 和$ {\boldsymbol{c}}_s $ 均时变.策略1为连续过程中采用的自优化控制方法, 对具有时变特性的间歇过程, 一般难以取得理想效果. 策略2采用恒定被控变量, 其设定值为动态轨线, 较策略1更适合间歇过程. 策略3进一步考虑具有切换结构的控制系统, 对离散化的间歇过程, 在
$ [t_i, t_{i+1}) $ 时间段内控制一组新的被控变量, 如图1所示.
图 1 间歇过程的离散化变量及自优化控制策略Fig. 1 Discretization of batch processes and self-optimizing control strategy结合间歇过程的时变特性, 本文主要研究策略2和策略3的被控变量求解问题. 对此, 引入如下假设条件:
假设1. 输出变量
$ {\boldsymbol{y}} $ 在时间轴$ [0,t_f] $ 上连续可测.假设2. 对一组选定的被控变量
$ {\boldsymbol{c}}(i) = H(i){\boldsymbol{y}} $ 及设定值$ {\boldsymbol{c}}_s(i) $ , 在对应的时间间隔$ [t_i,t_{i+1}) $ 内, 控制器可以将被控变量$ {\boldsymbol{c}}(i) $ 控制在其设定值$ {\boldsymbol{c}}_s(i) $ 上, 即$ \lim_{t\rightarrow t_{i+1}}H(i){\boldsymbol{y}}(t) = {\boldsymbol{c}}_s(i) $ .采用策略2时, 第1个时间段
$ [t_0,t_1) $ 内的被控变量$ {\boldsymbol{c}}(t) = H{\boldsymbol{y}}(t) $ , 设定值为$ c_s(1) $ , 至终点$ t_1 $ 实现$ H{\boldsymbol{y}}(1) = {\boldsymbol{c}}_s(1) $ ; 第2个时间段$ {\boldsymbol{c}}(t) = H{\boldsymbol{y}}(t) $ 的设定值变为$ {\boldsymbol{c}}_s(2) $ , 至终点$ t_2 $ 实现$ H{\boldsymbol{y}}(2) = {\boldsymbol{c}}_s(2) $ ; 以此类推.对策略2, 进一步考虑两种设定值选取方案:
方案1. 设定值轨线
$ [{\boldsymbol{c}}_s(1),\cdots,{\boldsymbol{c}}_s(N)] $ 为一组既定常数, 对给定的组合矩阵$ H $ , 设定值$ c_s(i) $ 为被控变量$ {\boldsymbol{c}} = H{\boldsymbol{y}}(t) $ 在各时间节点处的标称值$ {\boldsymbol{c}}_s^*(i) $ , 使标称工况的损失为0.方案2. 设定值
$ {\boldsymbol{c}}_s(i) $ 为当前批次运行至$ t_{i-1} $ 时刻的变量${\boldsymbol{y}}(0:i-1) = [{\boldsymbol{y}}^{\rm{T}}(0)\;\cdots\;{\boldsymbol{y}}^{\rm{T}}(i-1)]^{\rm{T}}$ 的函数, 不妨假定为线性关系, 记为$ {\boldsymbol{c}}_s(i) = {\boldsymbol{c}}_s^*(i)- H'(i){\boldsymbol{y}}(0:i-1) $ , 其中$ H'(i) $ 为待确定的系数矩阵.可以看到, 方案1中被控变量的设定轨线固定不变. 而方案2的被控变量设定轨线在当前批次运行过程中不断利用测量值进行修正. 相比方案1, 方案2更加充分地利用了过程信息, 理论上能进一步提高优化效果, 但需求解额外的决策变量
$ H'(i), i = 1,\cdots,N $ .为推导这两种方案中损失函数与组合矩阵
$ H $ 之间的关系, 定义如下超向量$$ \begin{split} &{\boldsymbol{\bar{u}}}^{\rm{T}} = [{\boldsymbol{u}}^{\rm{T}}(0)\;\;\;{\boldsymbol{u}}^{\rm{T}}(1)\; \cdots\; {\boldsymbol{u}}(N-1)^{\rm{T}}] \in {\bf R}^{n_{\bar{u}}= N_{u}}\\ &{\boldsymbol{\bar{d}}}^{\rm{T}} = \left[{\boldsymbol{d}}^{\rm{T}}(0)\;\;\; {\boldsymbol{d}}^{\rm{T}}(1) \;\cdots \;{\boldsymbol{d}}(N-1)^{\rm{T}}\right] \in {\bf R}^{n_{\bar{d}}= N n_{d}} \\ &{\boldsymbol{\bar{y}}}^{\rm{T}} = \left[{\boldsymbol{y}}^{\rm{T}}(0)\;\;\; {\boldsymbol{y}}^{\rm{T}}(1)\; \cdots\; {\boldsymbol{y}}(N)^{\rm{T}}\right] \in {\bf R}^{n_{\bar{y}}= (N+1) n_{y}} \end{split} $$ (13) 式中, 超向量
$ {\boldsymbol{\bar{u}}},{\boldsymbol{\bar{y}}},{\boldsymbol{\bar{d}}} $ 由时间轴$ [0, t_f] $ 上各离散点处的变量堆叠组成. 为便于描述, 将上述超向量的非线性映射关系记为$$ {\boldsymbol{\bar{y}}} = G({\boldsymbol{\bar{u}}},{\boldsymbol{\bar{d}}}) $$ (14) 式中, 映射函数
$ G $ 由式(12)中的状态方程$ {\boldsymbol{\hat f}} $ 和$ {\boldsymbol{\hat f}}_y $ 定义. 使用超向量, 间歇过程中变量的动态关系记为式(14)所示的静态函数.对方案1, 定义如下扩张组合矩阵
$ \bar H\in {\bf R}^{n_{\bar u}\times n_{\bar y}} $ 及总被控变量$ {\boldsymbol{\bar c}}\in {\bf R}^{Nn_{\bar u}} $ $$ \begin{split} &\bar{H} = \left[\begin{array}{ccccc} [0\;H] & 0 & \cdots & 0 \\ 0 & H & \cdots & 0 \\ \vdots & \vdots & \ddots &\vdots \\ 0 & 0 & \cdots & H \end{array}\right]\\ &{\boldsymbol{\bar{c}}} = \bar{H} \bar{y} = \left[\begin{array}{c} H{\boldsymbol{y}}(1) \\ H{\boldsymbol{y}}(2) \\ \vdots \\ H{\boldsymbol{y}}(N) \end{array}\right] \end{split} $$ (15) 方案2也可以定义相同维度的扩张组合矩阵
$ \bar H $ 及总被控变量$ {\boldsymbol{\bar c}} $ :$$ \begin{split}&\bar{H}=\left[\begin{array}{ccccc} H_{1}' & H & 0 & \cdots & 0 \\ \left[H_{2}'\right. & \rightarrow] & H & \cdots & 0\\ \vdots & \vdots& \vdots &\ddots & \vdots \\ \left[\leftarrow\right.& H_{N}'& \rightarrow&\rightarrow]& H \end{array}\right]\\ &{\boldsymbol{\bar{c}}} = \bar{H}{\boldsymbol{\bar{y}}} = \left[\begin{array}{c} H{\boldsymbol{y}}(1)+H_{1}' {\boldsymbol{y}}(0) \\ H{\boldsymbol{y}}(2)+H_{2}' {\boldsymbol{y}}(0: 1) \\ \vdots \\ H{\boldsymbol{y}}(N)+H_{N}' {\boldsymbol{y}}(0: N-1) \end{array}\right] \end{split} $$ (16) 式中,
$ [t_{i-1}, t_i) $ 的被控变量为$ {\boldsymbol{\bar{c}}}(i) = H {\boldsymbol{y}}(i)+ H_{i}'{\boldsymbol{y}}(0: i-1) $ , 等同于被控变量$ {\boldsymbol{\bar{c}}}(i) = H {\boldsymbol{y}}(i) $ , 且设定值修正量为$ -H_{i}'{\boldsymbol{y}}(0: i-1) $ , 因为在$ t-1 $ 时刻$ {\boldsymbol{y}}(0: i-1) $ 为已知量.引入扩张组合矩阵
$ \bar{H} $ 后, 总被控变量$ {\boldsymbol{\bar{c}}} $ 由每个时间节点的被控变量组成, 方案1 和方案2统一地描述为静态自优化控制问题的规范形. 结合第1节已有的结论, 可以得到损失函数$ L_{\rm{av}} $ 与$ \bar{H} $ 的关系表达式, 即求解如下最优化问题$$ \begin{split} & \min\limits_{\bar H} L_{\rm{av}} = 0.5\left\|V\left(\bar{H} G_{y}\right)^{-1} \bar{H} \tilde{F}\right\|_{\rm{F}}^{2}\\ &\; {\rm{s.t.}} \quad 式\;(15)\;或式\;(16) \end{split} $$ (17) 式中,
$ V,G_y $ 和$ \tilde F $ 等矩阵均定义为超向量之间的关系.上述优化问题和静态问题(9)之间的不同之处在于, 此处扩张组合矩阵
$ \bar H $ 应满足式(15)或式(16)所示的结构, 即矩阵$ \bar H $ 具有结构性约束. 式(10)所示的解析解不能直接推广至式(17)求解, 否则不满足约束条件(15)或(16).采用策略3时, 在第1个时间段
$ [t_0,t_1) $ 内, 控制器控制被控变量$ {\boldsymbol{c}}(t) = H(1){\boldsymbol{y}}(t) $ , 设定值为$ {\boldsymbol{c}}_s(1) $ , 至终点$ t_1 $ 实现$ H(1){\boldsymbol{y}}(1) = {\boldsymbol{c}}_s(1) $ ; 在第2个时间段$ [t_1,t_2) $ 内, 被控变量切换为$ {\boldsymbol{c}}(t) = H(2){\boldsymbol{y}}(t) $ , 设定值为$ {\boldsymbol{c}}_s(2) $ , 至终点$ t_2 $ 实现$ H(2){\boldsymbol{y}}(2) = {\boldsymbol{c}}_s(2) $ ; 以此类推. 类似地, 对策略3也考虑两种方案:方案3. 每组被控变量
$ {\boldsymbol{c}}(i) $ 对应的设定值$ [{\boldsymbol{c}}_s(1),\cdots,{\boldsymbol{c}}_s(N)] $ 为既定常数, 为$ H(i){\boldsymbol{y}}(i) $ 在时间节点$ t_i $ 处的标称值$ {\boldsymbol{c}}_s^*(i) $ .方案4. 设定值
$ {\boldsymbol{c}}_s(i) $ 为当前批次运行至$ t_{i-1} $ 时刻的变量$ {\boldsymbol{y}}(0:i-1) = [{\boldsymbol{y}}^{\rm{T}}(0)\;\cdots\;{\boldsymbol{y}}^{\rm{T}}(i-1)]^{\rm{T}} $ 的函数, 记为$ {\boldsymbol{c}}_s(i) = {\boldsymbol{c}}_s^*(i)-H'(i) {\boldsymbol{y}}(0:i-1) $ , 其中$ H'(i) $ 为待确定的系数矩阵.同理, 对策略3的两种方案也分别定义扩张组合矩阵
$ \bar H $ 及总被控变量$ {\boldsymbol{\bar c}} $ .方案3中,
$$\begin{split} &\bar{H} = \left[\begin{array}{ccccc} [0\;H(1)] & 0 & \cdots & 0 \\ 0 & H(2) & \cdots & 0 \\ \vdots & \vdots & \ddots &\vdots \\ 0 & 0 & \cdots & H(N) \end{array}\right] \\ &{\boldsymbol{\bar{c}}} = \bar{H} \bar{y} = \left[\begin{array}{c} H(1){\boldsymbol{y}}(1) \\ H(2){\boldsymbol{y}}(2) \\ \vdots \\ H(N){\boldsymbol{y}}(N) \end{array}\right] \end{split} $$ (18) 方案4中,
$$\begin{split} &\bar{H}=\left[\begin{array}{ccccc} H_{1}' & H(1) & 0 & \cdots & 0 \\ \left[H_{2}'\right. & \rightarrow] & H(2) & \cdots & 0\\ \vdots & \vdots& \vdots &\ddots & \vdots \\ \left[\leftarrow\right.& H_{N}'& \rightarrow&\rightarrow]& H(N) \end{array}\right]\\ &{\boldsymbol{\bar{c}}} = \bar{H}{\boldsymbol{\bar{y}}} = \left[\begin{array}{c} H(1){\boldsymbol{y}}(1)+H_{1}' {\boldsymbol{y}}(0) \\ H(2){\boldsymbol{y}}(2)+H_{2}' {\boldsymbol{y}}(0: 1) \\ \vdots \\ H(N){\boldsymbol{y}}(N)+H_{N}' {\boldsymbol{y}}(0: N-1) \end{array}\right]\end{split} $$ (19) 与策略2相比, 策略3中两种方案的组合矩阵
$ H $ 是时变的, 即需求取$ N $ 个组合矩阵$ H(i), i = 1,\cdots,N $ . 同理, 对策略3求解如下最优化问题$$ \begin{split} &\min\limits_{\bar H} L_{\rm{av}} = 0.5\left\|V\left(\bar{H} G_{y}\right)^{-1} \bar{H} \tilde{F}\right\|_{\rm{F}}^{2}\\& \;{\rm s.t.} \quad 式\;(18)\;或式\;(19) \end{split} $$ (20) 从以上分析看到, 对不同的控制策略和设定值选取方案, 可以统一归结为具有不同结构的扩张组合矩阵
$ \bar H $ 的求解问题, 可以在优化问题中对$ \bar H $ 施加等式约束实现. 一般来说, 具有特定结构的组合矩阵难以求得闭合解, 需使用数值优化算法.注1. 以上提出的4种被控变量选择方案, 从控制角度看, 执行策略2 (方案1)最简单, 但优化效果可能较差; 策略3 (方案4)理论上的优化效果最好, 但被控变量需要不断切换, 并且设定轨线也要在线修正. 针对具体过程, 需结合过程特性和优化性能结果综合考虑这两个因素, 选择最合理的自优化控制方案.
2.3 策略3 (方案4)的
${\boldsymbol{ \bar H}} $ 解析解下面提出一种针对策略3 (方案4)的闭合解求解方法. 如式(19)所示, 此时
$ \bar H $ 为块下三角矩阵. 为表述方便, 将式(19)所示的$ \bar H $ 表达式记为$$ \bar{H}=\left[\begin{array}{cccc} \bar H_1 & 0 & \cdots & 0 \\ \left[\bar H_2\right. & \rightarrow]& \cdots & 0\\ \vdots & \vdots &\ddots & \vdots \\ \left[\leftarrow\right.& \bar H_N &\rightarrow & \rightarrow] \end{array}\right] $$ (21) 式中, 子矩阵
$\bar{H}_i =\left[H_{i}'\;H(i)\right] \in {\bf R}^{n_{u} \times i n_{y}}$ , 同时包含了$ i $ 时刻的被控变量组合矩阵$ H(i) $ 及修正设定值轨线的系数矩阵$ H_{i}' $ .引理 2. 对满足式(21)结构的
$ \bar{H} $ 及非奇异块下三角矩阵$ Q $ , 转化矩阵$ \bar H' = Q\bar{H} $ 同样满足式(21)结构, 并且$ L_{\rm{av}}(\bar H) = L_{\rm{av}}(\bar H') $ .证明. 由于
$ \bar H $ 和$ Q $ 均为块下三角, 显然$ \bar H' $ 也为块下三角矩阵. 将$ \bar H' = Q\bar{H} $ 代入到损失函数表达式$L_{\rm{av}}(\bar{H}') = 0.5\|V(\bar{H}'G_{y})^{-1} \bar{H}'\tilde{F}\|_{\rm{F}}^{2}$ 中,$ Q $ 矩阵前后互消, 结论成立. □与引理1类似, 引理2也可用于先求解
$ \bar H $ 的特解. 注意到敏感矩阵$ G_y $ 为块下三角矩阵, 因此$ \bar HG_y $ 的逆也为块下三角. 将$ V $ 取为满足$V^{\rm{T}}V = J_{uu}$ 的块下三角矩阵, 可以对$ J_{uu} $ 进行Cholesky分解得到.定理 1. 对策略3 (方案4)的
$ \bar H $ 矩阵, 式(20)等同于求解如下问题$$ \begin{split} &\min\limits_{\bar{H}(1), \cdots, \bar H(N)} L_{\rm{av}} = 0.5 \sum\limits_{i = 1}^{N}\left\|\bar{H}(i) \tilde{F}_{i}\right\|_{\rm{F}}^{2}\\ &{\rm{s.t.}} \quad \bar{H}(i) G_{y i} = V_{i}, \;\forall i = 1, \cdots, N \end{split} $$ (22) 式中,
$ \tilde{F}_{i} $ 为$ \tilde{F} $ 的子矩阵 (前$ n_yi $ 行, 前$ n_d+n_yi $ 列),$ G_{yi} $ 为$ G_y $ 的子矩阵 (前$ n_yi $ 行, 前$ n_ui $ 列),$ V_i $ 为$ V $ 的第$ i $ 个分块矩阵.证明. 根据引理2, 可选择任意非奇异的块下三角矩阵
$ Q $ 对矩阵$ \bar H $ 进行转化求取特解, 可选择$ Q = V(\bar HG_y)^{-1} $ , 使$ \bar H' = Q\bar{H} $ 满足$$ \bar{H}'G_{y} = V\left(\bar{H} G_{y}\right)^{-1}\bar{H}G_{y} = V $$ (23) 即对
$ \forall i = 1,\cdots,N $ , 均满足$$ \bar{H}'({{i}}) G_{y i} = V_{i} $$ (24) 不失一般性, 式(24)可作为对决策变量
$ \bar H $ 的约束加入到优化问题中. 此时$$ \begin{split} L_{\rm{av}} =\;& 0.5\left\|V(\bar{H} G_{y})^{-1} \bar{H}\tilde{F}\right\|_{\rm{F}}^{2} = 0.5\left\| \bar{H} \tilde{F}\right\|_{\rm{F}}^{2} = \\ & 0.5 \sum\limits_{i = 1}^{N}\left\| \bar{H}(i) \tilde{F}\right\|_{\rm{F}}^{2} \end{split} $$ (25) □
通过合理利用转化矩阵
$ Q $ , 定理1将目标函数及约束条件分解到每个离散时间节点, 能够沿时间轴依次求解出子矩阵$ \bar H(i) $ . 对$\forall i = 1,\cdots,N$ , 求解如下优化问题$$ \begin{split} &\min\limits_{\bar{H}(i)}\;\;0.5\left\|H_{i} \tilde{F}_{i}\right\|_{\rm{F}}^{2}\\ &{\rm{s.t.}}\;\;\bar{H}(i) G_{y i} = V_{i} \end{split} $$ (26) 式(26)为带等式约束的二次型凸优化问题, 可进一步求得解析解.
定理 2. 对式(26)所示的带等式约束的二次型凸优化问题, 其闭合解为
$$ \bar{H}(i)^{\rm{T}} = (\tilde{F}_{i} \tilde{F}_{i}^{\rm{T}})^{-1} G_{yi}\left(G_{y i}^{\rm{T}}(\tilde{F}_{i} \tilde{F}_{i}^{\rm{T}})^{-1} G_{y i}\right)^{-1} V_{i}^{\rm{T}} $$ (27) 证明. 式(26)在形式上与第1节静态自优化控制问题一致, 闭合解(27)的推导过程可参见文献[28]. □
综上, 本文求取最优扩张组合矩阵
$ \bar H $ 的计算步骤如图2所示, 其中策略3 (方案4)可直接应用定理2 求得闭合解, 其他3种情况则需使用数值优化法求取. 由于目标函数$ L_{\rm{av}} $ 是$ \bar H $ 的非线性函数, 优化问题(17)和(20)不能保证得到全局最优解. 对此, 策略3 (方案4)得到的解析解可作为数值优化的初始解进行寻优.
图 2 最优扩张组合矩阵$\bar H$ 的求解步骤Fig. 2 Procedure for solving the optimal extended combination matrix$\bar H$ 3. 仿真研究
3.1 间歇反应器描述
本节研究一个带副反应的间歇反应器, 主副反应分别为
$ A+B\rightarrow C $ 和$ 2B\rightarrow D $ , 其中反应物$ A $ 在初始时刻投放完毕,$ B $ 在反应过程中实时投放, 实时流量为操纵变量$ u(t) $ . 体系的模型方程为$$ \frac{{\rm{d}}c_A}{{\rm{d}}t} = -k_1c_Ac_B-\frac{c_Au}{V},\quad c_A(0) = c_{A0} $$ (28) $$\begin{split}& \frac{{\rm{d}}c_B}{{\rm{d}}t} = -k_1c_Ac_B-2k_2c_B^2-(c_B-c_{Bin})\frac{u}{V}, \\ & c_B(0) = c_{B0} \end{split} $$ (29) $$ \frac{{\rm{d}}V}{{\rm{d}}t} = u, \quad V(0) = V_0 $$ (30) $$c_C = \frac{c_{A0}V_0-c_AV}{V} $$ (31) $$ c_D = \frac{c_A+c_{Bin}-c_B}{2}-\frac{c_{A0}+c_{Bin}-c_{B0}}{2V} $$ (32) 式中,
$ c_X $ 表示物料$ X $ 的浓度,$ V $ 为持液量, 其他符号含义及标称值列于表1.表 1 间歇反应器参数及标称值Table 1 Parameters for the reactor model and nominal values符号 物理含义 标称值 $ k_1 $ 主反应的反应常数 0.053 L·mol/min $ k_2 $ 副反应的反应常数 0.128 L·mol/min $ u_L $ $ u $下限 0 L/min $ u_U $ $ u $上限 0.001 L/min $ c_{Bin} $ B 进料浓度 5 mol/L $ c_{Ao} $ A 初始浓度 0.72 mol/L $ c_{Bo} $ B 初始浓度 0.0614 mol/L $ V_o $ V 初始值 1 L $ t_f $ 批次运行时间 250 min 操作目标为在
$ [0, t_f] $ 操作时段内最大化产物产量$ C $ 的同时减少副产物$ D $ , 即表示为如下优化问题$$ \begin{split} \max\limits_{u(t)} J& = [c_C(t_f)-c_D(t_f)]V(t_f) \\ {\rm{s.t}}. \quad & 0\leq u(t)\leq 0.001 \;{\rm{L/min}} \end{split} $$ (33) 在表1所示的标称工况下, 使用数值优化方法求解式(33)可得到
$ u(t) $ 的最优输入轨迹(图3). 可以看到, 此时$ u(t) $ 整个轨线处于可行域内, 最优值$ J^{\rm{opt}} $ = 0.271687 mol. 反应常数$ k_1 $ 和$ k_2 $ 为不确定扰动, 变化范围为其标称值的$ \pm $ 40%. 当$ k_1 $ 和$ k_2 $ 变化时,$ u(t) $ 的最优输入轨迹随之改变.4. 被控变量计算示例
为更清晰地阐述本文方法, 以
$ N = 2 $ 为例(即$ [0, t_f] $ 被均匀离散为两段), 介绍如何使用第2节中的方法求解不同被控变量. 离散后的优化变量个数$ n_{\bar u} = 2 $ , 对式(33)进行重优化后得到Hessian矩阵和$ V $ 矩阵$$ \begin{split} & J_{\bar u\bar u} = \left[\begin{array}{cc} 3.70 & 1.74 \\ 1.74 & 3.47 \end{array}\right] \times 10^{5}\\ & V = \left[\begin{array}{cc} 532.2 & 0 \\ 294.9 & 589.0 \end{array}\right] \end{split} $$ (34) 考虑使用
$ c_A $ 和$ c_B $ 构造被控变量, 对离散系统进行线性化, 得到$$ {\boldsymbol{\bar{y}}} = G_{y}{\boldsymbol{\bar{u}}}+G_{yd}{\boldsymbol{d}} $$ (35) 式中
$$ \begin{split} & G_{y} = \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \\ -264.45 & 0 \\ 88.27 & 0 \\ -210.43 & -181.20 \\ 10.48 & 98.16 \end{array}\right] \\ & G_{y d} = \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \\ -2.07 & 0.32 \\ -0.41 & -0.15 \\ -2.90 & 0.52 \\ -0.21 & -0.19 \end{array}\right]\\ & F = \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \\ -4.03 & 0.81 \\ 0.24 & -0.31 \\ -4.99 & 1.02 \\ 0.16 & -0.27 \end{array}\right] \end{split} $$ (36) 其中, 测量变量
${\boldsymbol{\bar{y}}}^{\rm{T}} = \left[y^{\rm{T}}(0)\;\; y^{\rm{T}}(1)\;\; y^{\rm{T}}(2)\right]$ 为$ c_A $ 、$ c_B $ 分别在0, 125 min及250 min时刻的量组成. 得到上述矩阵后, 可以构造式(17)和式(20)所示的优化问题来求解被控变量, 结果如下.1) 策略2 (方案1):
$ H = [-0.0026\;\;0.0035] $ , 即整个时间段内都控制被控变量$ c(t) = -0.0026 c_A+ 0.0035 c_B $ . 经计算, 前125 min的设定值为$ c_s(1) = -0.000303 $ , 后125 min的设定值为$ c_s(2) = -0.000059 $ .2)策略2 (方案2): 求解得到的扩张组合矩阵
$ \bar H $ 为$$ \bar{H} = \left[\begin{array}{cccccc} 0 & 0 & -1.11 & 2.70 & 0 & 0 \\ 0 & 0 & \,\;\;0.51 & 1.93 & -1.11 & 2.70 \end{array}\right] $$ 即整个时间段内, 被控变量为
$ c(t) = -1.11 c_A+ 2.70 c_B $ . 前125 min的设定值为$ c_s(1) = -0.33 $ , 后125 min设定值为$ c_s(2) = 0.20-0.51 c_A(1)-1.93 c_B(1) $ .3)策略3 (方案3):
$ H(1) = [-0.0019\;\;\;0.0057] $ ,$H(2) = [-0.0015\;\;\;0.0074]$ , 即前125 min 被控变量$ c(1) = -0.0019 c_A+0.0057 c_B $ , 后125 min被控变量$ c(2) = -0.0015 c_A+0.0074 c_B $ , 其设定值分别为$c_s(1) =-0.00048$ 和$c_s(2) = -0.000034$ .4)策略3 (方案4): 根据定理2, 求得扩张矩阵
$ \bar H $ 为$$ \bar{H} = \left[\begin{array}{cccccc} 0 & 0 & -1.06 & 2.85 & 0 & 0 \\ 0 & 0 &\;\;\, 0.88 & 2.07 & -1.48 & 3.27 \end{array}\right] $$ 即前125 min被控变量
$ c(1) = -1.06 c_A+ 2.85 c_B $ , 设定值$ c_s(1) = -0.29 $ ; 后125 min被控变量$c(2) = -1.48 c_A+3.27 c_B$ , 设定值$c_s(2) = 0.31 - 0.88 c_A(1)- 2.07 c_B(1)$ .4.1 批内自优化控制效果
由于
$ N = 2 $ 难以逼近整个间歇操作过程, 后文设置$ N = 20 $ 并以相同的方法重新求解被控变量, 同时, 在测量变量中加入体积变量$ V $ 提高优化效果. 从表2可观察到:表 2 损失函数$ L_{\rm{av}} $ Table 2 Loss function$ L_{\rm{av}} $ 策略及方案 $ N = 2 $ $ N = 20 $ 策略 2 (方案 1) 0.0371 0.0083 策略 2 (方案 2) 0.03423 0.0024 策略 3 (方案 3) 0.0368 0.0069 策略 3 (方案 4) 0.03420 0.0022 1) 4种方案的损失
$ L_{\rm{av}} $ 在$N=20 $ 时, 相比$ N = 2 $ 都大幅度降低;2) 策略2 (方案1)的损失函数为0.0083, 策略2 (方案2)通过在线设定值修正, 进一步将损失减少到0.0024;
3) 策略3 (方案3)的损失为0.0069, 略低于策略2 (方案1);
4) 策略2 (方案3)的损失为 0.0024, 与策略3 (方案4)的损失0.0022很接近, 表明不切换被控变量也能得到较好的优化控制效果.
基于表2的结果, 策略2 (方案2)与策略3 (方案4)效果接近, 但前者无需在线切换被控变量, 更易于在线控制, 因此考虑使用策略2 (方案2)对该反应器进行批内自优化控制. 此外, 动态仿真中将与策略2 (方案1)的结果进行对比, 有助于进一步理解本文方法.
策略2 (方案1)的被控变量为
$ c_1(t) = 0.0062 c_A+ 0.002 c_B+0.0831 V $ , 设定值轨线如图4所示. 为进一步获取平滑的设定值轨线, 使操作更为平稳, 对这些离散点进行回归分析, 得到平滑的设定值轨迹方程$ c_s(t) = 0.0877+3.705\times 10^{-5}t-1.97\times 10^{-8}t^2 $ , 为一条随时间$ t $ 变化的连续曲线, 如图4所示. 对该系统可以采用普通的PI控制器对被控变量$ c_1(t) $ 进行跟踪控制.策略2 (方案2)的被控变量为
$ c_2(t) = 0.0026 c_A+ 0.00032 c_B+0.0830 V $ , 设定值轨线在每批次运行过程中采集测量值进行在线修正. 为增强操作平稳性, 在$ t_k $ 时刻计算得到$ t_{k+1} $ 时刻的设定点后, 在$ [t_k,t_{k+1}] $ 时间段内设置斜坡形设定值轨线, 使设定轨线维持连续性. 同样使用PI控制器跟踪控制得到的被控变量$ c_2(t) $ .不确定参数
$ k_1 $ 和$ k_2 $ 分别改变 +20%和 −20%时的优化控制效果如图5所示, 从图5(a)中可以看到, 两种方法分别对$ c_1(t) $ 和$ c_2(t) $ 都实现了较好的闭环跟踪控制, 其中,$ c_2(t) $ 的设定轨线根据批内采集到的测量值进行了调整, 相比自身的标称轨线有一定程度的上移; 图5(b)显示不同方法的控制输入$ u(t) $ 轨迹, 其中, 控制$ c_1(t) $ 时的$ u(t) $ 轨迹相比标称操作更靠近当前工况真实的最优轨线, 性能指标$ J $ 有所提高$( J = 0.34374 \rightarrow 0.34505 )$ , 显示出一定的优化控制效果. 控制$ c_2(t) $ 时的$ u(t) $ 轨迹更靠近最优轨线, 其性能指标$ J = 0.34701 $ 和最优值$ J^{\rm{opt}} = 0.34755 $ 差别不大. 同时注意到控制$ c_2(t) $ 时的$ u(t) $ 轨迹振荡更加剧烈, 这是因为$ c_2(t) $ 的设定轨线不断在线修正, 为了得到满意的控制效果, 使用了高增益PI控制器$( K_p = 20) $ . 这并不影响最终得到满意的优化效果$( L = 0.00054) $ , 从另一个角度说明了间歇过程中控制关键变量的重要性.
图 5 批内自优化控制效果$( k_1 $ : +20%,$ k_2 $ : −20%)Fig. 5 Within-batch self-optimizing performance$ (k_1 $ : +20%,$ k_2 $ : −20%)不确定参数
$ k_1 $ 和$ k_2 $ 分别改变 −40%和 +40%时的优化控制效果如图6所示, 此时系统的不确定性向另一个方向变化, 并且幅度更大. 从图6(a)中可以看到, 两种方法同样对$ c_1(t) $ 和$ c_2(t) $ 都实现了较好的闭环跟踪控制, 其中,$ c_1(t) $ 的设定轨线不变, 而$ c_2(t) $ 的设定轨线相比自身的标称轨线有一定程度的下移. 从图6(b)来看, 虽然控制$ c_1(t) $ 能将$ u(t) $ 轨迹向着真实的最优轨线的方向调节, 其性能指标$ J $ 从标称操作的0.09646提高到0.10312, 但作用有限, 距离最优值$ J^{{\rm{opt}}} = 0.12252 $ 仍有较大差距. 控制$ c_2(t) $ 进一步提高了优化控制效果, 其性能指标为$ J = 0.11602 $ , 相比最优性能只有0.006的损失(此时$ k_1 $ ,$ k_2 $ 的变化较大, 该损失在一定程度上由系统的非线性导致). 此外, 控制$ c_2(t) $ 时的$ u(t) $ 轨迹同样振荡较为剧烈, 但随反应进行,$ u(t) $ 大致围绕着最优轨线上下波动.
图 6 批内自优化控制效果$( k_1 $ : +40%,$ k_2 $ : −40%)Fig. 6 Within-batch self-optimizing performance$( k_1 $ : +40%,$ k_2 $ : −40%)表3进一步统计了100组随机扰动下各方法的非线性损失, 其中随机扰动
$ [k_1\; k_2] $ 均匀分布在各自的变化范围. 可以看到, 相比标称操作(平均损失0.0036)和以单变量$ c_B $ (平均损失0.0042)为被控变量的情形, 两种批内自优化控制方法有效提高了经济性能, 其中, 策略2 (方案1)中控制$ c_1(t) $ 将平均损失减少到0.0026, 策略2 (方案2)中控制$ c_2(t) $ 进一步将平均损失减少到0.0007, 几乎可以忽略不计. 此外, 最大损失和标准差等统计量也呈现出相同的变化趋势, 如表3所示.表 3 100组随机扰动下的非线性损失统计量Table 3 Statistics of nonlinear losses for 100 groups of random disturbances方案 平均损失 最大损失 标准差 标称操作 0.0036 0.0227 0.0068 控制$ c_B $ 0.0042 0.0165 0.0054 策略 2 (方案 1) 0.0026 0.0167 0.0050 策略 2 (方案 2) 0.0007 0.0053 0.0016 5. 结束语
本文研究了间歇过程的批内自优化控制问题, 在单批次运行过程中控制一组虚拟的被控变量(输出变量的线性组合), 实现间歇过程的实时优化. 对此, 给出了两种自优化控制策略(被控变量恒定但设定值时变; 被控变量和设定值均时变). 对它们的设定值选取问题又分别提出两种方案(设定值轨线固定不变; 设定值轨线在线修正), 共计4种方法. 通过引入扩张组合矩阵
$ \bar H $ , 将这4种方法统一描述为具有不同结构约束的最优$ \bar H $ 求解问题, 并推导得到了策略3 (方案4)的$ \bar H $ 解析解计算方法(定理2).本文提出的4种被控变量选择方法, 其对应的闭环控制系统具有不同的复杂度和优化性能. 针对一般的实际间歇过程, 应综合考虑这两个因素并取得合理权衡. 间歇反应器的仿真研究中, 采用策略2 (方案2) (恒定被控变量:
$ c_2(t)) $ 得到的控制结构较为简单, 并且能通过在线修正$ c_2(t) $ 的设定值增强优化效果, 是较为合理的方案. -
表 1 Yale数据集下平均分类率 (40%遮盖噪声) (%)
Table 1 Average classification rate under Yale dataset (40% occluded noise) (%)
NPC PCA PCA-$L_1$ 2DPCA 2DPCA-$L_1$ 2DPCA-
$L_1$(non-greedy)Angle-
2DPCAGC2DPCA
S = 0.5GC2DPCA
S = 1GC2DPCA
S = 1.5GC2DPCA
S = 210 78.9 77.8 84.4 84.4 72.9 73.3 84.2 84.7 84.7 84.9 20 75.3 76.4 82.2 83.1 74.2 70.9 83.3 83.3 83.3 83.3 30 72.0 73.1 79.3 80.0 74.2 73.1 80.7 80.2 80.2 80.2 40 71.6 71.6 76.2 77.1 73.8 74.0 77.8 78.0 77.8 78.0 50 72.0 71.8 74.4 75.1 74.2 75.8 75.3 75.8 76.2 76.0 
下载: 导出CSV
表 2 ORL数据集下平均分类率 (40%遮盖噪声) (%)
Table 2 Average classification rate under ORL dataset (40% occluded noise) (%)
NPC PCA PCA-$L_1$ 2DPCA 2DPCA-$L_1$ 2DPCA-
$L_1$(non-greedy)Angle-
2DPCAGC2DPCA
S = 0.5GC2DPCA
S = 1GC2DPCA
S = 1.5GC2DPCA
S = 210 79.9 79.3 88.1 88.1 77.2 80.3 88.3 88.5 88.6 88.5 20 81.2 80.4 85.9 85.9 82 80.6 86.3 86.4 86.4 86.4 30 81.6 82.3 85.1 85.5 84.8 83.3 85.5 85.6 85.7 85.7 40 81.9 82.8 85.4 85.5 85.7 85.4 85.5 85.5 85.5 85.5 50 81.7 83.2 85.6 85.5 85.6 85.5 85.4 85.6 85.6 85.6
下载: 导出CSV
表 3 FERET数据集下平均分类率 (40%遮盖噪声) (%)
Table 3 Average classification rate under FERET dataset (40% occluded noise) (%)
NPC PCA PCA-$L_1$ 2DPCA 2DPCA-$L_1$ 2DPCA-
$L_1$(non-greedy)Angle-
2DPCAGC2DPCA
S = 0.5GC2DPCA
S = 1GC2DPCA
S = 1.5GC2DPCA
S = 210 50.5 52.5 58.0 58.5 47.8 51.0 59.0 59.8 60.0 60.0 20 50.0 51.0 52.5 53.5 50.8 51.3 54.0 54.3 54.5 54.3 30 51.3 51.3 51.5 50.8 51.8 52.8 50.8 51.8 51.8 51.8 40 48.8 49.8 51.3 50.8 50.5 50.5 50.8 50.8 50.8 50.8 50 48.0 49.5 51.3 50.8 50.8 50.8 50.0 50.0 50.0 50.0
下载: 导出CSV
表 4 Yale数据集下平均分类率 (60%遮盖噪声) (%)
Table 4 Average classification rate under Yale dataset (60% occluded noise) (%)
NPC PCA PCA-$L_1$ 2DPCA 2DPCA-$L_1$ 2DPCA-
$L_1$(non-greedy)Angle-
2DPCAGC2DPCA
S = 0.5GC2DPCA
S = 1GC2DPCA
S = 1.5GC2DPCA
S = 210 61.8 58.7 72.4 73.6 62.4 59.3 73.6 74.4 74.4 74.4 20 60.4 60.0 63.3 65.3 62.9 59.1 66.4 67.6 67.8 67.8 30 60.7 61.3 62.2 62.7 62.4 61.1 63.6 63.6 63.6 63.6 40 61.1 61.6 62.2 62.7 62.2 62.9 62.2 62.4 62.2 62.4 50 60.7 62.0 62.2 62.7 62.4 62.7 62.4 62.4 62.4 62.4
下载: 导出CSV
表 5 ORL数据集平均分类率 (60%遮盖噪声) (%)
Table 5 Average classification rate under ORL dataset (60% occluded noise) (%)
NPC PCA PCA-$L_1$ 2DPCA 2DPCA-$L_1$ 2DPCA-
$L_1$(non-greedy)Angle-
2DPCAGC2DPCA
S = 0.5GC2DPCA
S = 1GC2DPCA
S = 1.5GC2DPCA
S = 210 66.0 65.0 73.6 74.4 72.1 67.0 74.6 74.7 74.6 74.7 20 68.1 67.7 73.4 74.2 71.3 70.8 74.1 74.3 74.3 74.2 30 68.7 69.4 73.7 74.2 72.4 72.9 74.4 74.4 74.3 74.4 40 69.0 70.0 74.4 74.4 73.5 75.0 74.5 74.4 74.3 74.4 50 69.2 70.2 74.1 74.4 74.2 74.4 74.3 74.5 74.4 74.4
下载: 导出CSV
表 6 FERET数据集平均分类率 (60%遮盖噪声) (%)
Table 6 Average classification rate under FERET dataset (60% occluded noise) (%)
NPC PCA PCA-$L_1$ 2DPCA 2DPCA-$L_1$ 2DPCA-
$L_1$(non-greedy)Angle-
2DPCAGC2DPCA
S = 0.5GC2DPCA
S = 1GC2DPCA
S = 1.5GC2DPCA
S = 210 43.8 43.5 50.0 47.8 40.0 40.3 48.8 50.5 50.3 50.0 20 43.8 44.3 48.5 45.0 44.8 45.0 44.8 45.3 45.8 46.0 30 44.0 43.8 46.8 45.0 44.8 45.3 45.0 45.0 45.5 45.5 40 41.5 44.3 45.3 45.0 44.0 44.3 45.0 45.3 45.3 45.3 50 42.0 42.0 44.8 45.0 45.0 45.3 45.3 45.0 45.0 45.0
下载: 导出CSV
表 7 时间复杂度分析
Table 7 Time complexity analysis
算法 时间复杂度 PCA ${{\rm{O}}}({Nh}^{2}{w}^{2}+{h}^{3}{w}^{2}) $ PCA-$L_1$ ${ {\rm{O} } }({Nhwk}{t}_{{\rm{PCA}}\text{-}L_{1} })$ 2DPCA ${{\rm{O}}}({Nh}^{2}{w}+{h}^{3}) $ 2DPCA-$L_1$ ${ {\rm{O} } }({Nhwk} {t}_{{\rm{2DPCA}}\text{-}L_{1} })$ 2DPCA-$L_1$(non-greedy) ${{\rm{O}}}[({Nhw}^{2}+{wk}^{2}+{w}^{2}{k})$
${t}_{{\rm{2DPCA}}\text{-}L_{1}({\rm{non}}\text{-}{\rm{greedy}})}]$Angle-2DPCA ${ {\rm{O} } }[({Nhw}^{2}+{wk}^{2}+{w}^{2}{k}){t}_{{\rm{Angle2DPCA}}}]$ GC2DPCA ${ {\rm{O} } }({Nhwkt}_{{\rm{GC2DPCA}}})$
下载: 导出CSV
表 8 各个算法所需时间对比(s)
Table 8 Comparison of the required time by each algorithm (s)
样本数 PCA PCA-$L_1$ 2DPCA 2DPCA-$L_1$ 2DPCA-$L_1$ (non-greedy) Angle-2DPCA GC2DPCA 20 0.855 0.576 0.310 2.861 9.672 5.585 2.637 40 1.643 0.882 0.363 5.834 17.550 7.753 5.413 60 3.072 1.326 0.452 8.611 28.096 16.567 9.049 80 4.733 1.940 0.507 17.800 35.896 14.040 14.883 100 6.732 2.226 0.620 22.449 88.671 40.092 18.517
下载: 导出CSV
-
[1] Jolliffe I T, Cadima J. Principal component analysis: a review and recent developments. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2016, 374(2065): 20150202. DOI: 10.1098/rsta.2015.0202 [2] Yang J, Zhang D D, Frangi A F, Yang J Y. Two-Dimensional PCA: a new approach to appearance-based face representation and recognition. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2004, 26(1):131-137 doi: 10.1109/TPAMI.2004.1261097 [3] Ji Y, Xie H B. Myoelectric signal classification based on S transform and two-directional two-dimensional principal component analysis. Transactions of the Institute of Measurement and Control, 2018, 40(7): 2387-2395 doi: 10.1177/0142331217704035 [4] Yang W, Sun C, Ricanek K. Sequential row-column 2DPCA for face recognition. Neural Computing and Applications, 2012, 21(7): 1729-1735 doi: 10.1007/s00521-011-0676-5 [5] Zhang D, Zhou Z H, Chen S. Diagonal principal component analysis for face recognition. Pattern Recognition, 2006, 39(1): 140-142 doi: 10.1016/j.patcog.2005.08.002 [6] Wang H, Hu J, Deng W. Face feature extraction: a complete review. IEEE Access, 2018, 6: 6001-6039 doi: 10.1109/ACCESS.2017.2784842 [7] Tirumala S S, Shahamiri S R, Garhwal A S, Wang R. Speaker identification features extraction methods: a systematic review. Expert Systems with Applications, 2017, 90: 250-271 doi: 10.1016/j.eswa.2017.08.015 [8] Soora N R, Deshpande P S. Review of feature extraction techniques for character recognition. IETE Journal of Research, 2018, 64(2): 280-295 doi: 10.1080/03772063.2017.1351323 [9] Neiva D H, Zanchettin C. Gesture recognition: a review focusing on sign language in a mobile context. Expert Systems with Applications, 2018, 103: 159-183 doi: 10.1016/j.eswa.2018.01.051 [10] 史加荣, 周水生, 郑秀云. 多线性鲁棒主成分分析. 电子学报, 2014, 42(08): 1480-1486 doi: 10.3969/j.issn.0372-2112.2014.08.004Shi Jia-Rong, Zhou Shui-Sheng, Zheng Xiu-Yun. Multilinear robust principal component analysis. Acta Electronica Sinica, 2014, 42(08): 1480-1486 doi: 10.3969/j.issn.0372-2112.2014.08.004 [11] Namrata V, Thierry B, Sajid J, Praneeth N. Robust subspace learning: robust PCA, robust subspace tracking, and robust subspace recovery. IEEE Signal Processing Magazine, 2018, 35(4): 32-55 doi: 10.1109/MSP.2018.2826566 [12] Menon V, Kalyant S. Structured and unstructured outlier identification for robust PCA: a fast parameter free algorithm. IEEE Transactions on Signal Processing, 2019, 67(9): 2439-2452 doi: 10.1109/TSP.2019.2905826 [13] Wang W Q, Aggarwal V, Aeron S. Principal component analysis with tensor train subspace. Pattern Recognition Letters, 2019, 122: 86-91 doi: 10.1016/j.patrec.2019.02.024 [14] Markopoulos P P, Chachlakis D G, Papalexakis E E. The exact solution to rank-1 L1-norm TUCKER2 decomposition. IEEE Signal Processing Letters, 2018, 25(4): 511-515 doi: 10.1109/LSP.2018.2790901 [15] Kwak N. Principal component analysis based on L1-norm maximization. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2008, 30(9): 1672-1680 doi: 10.1109/TPAMI.2008.114 [16] Nie F P, Huang H, Ding C, Luo D, Wang H. Robust principal component analysis with non-greedy L1-norm maximization. In: Proceedings of the 22nd International Joint Conference on Artificial Intelligence. Barcelona, Spain: AAAI, 2011. 1433−1438 [17] Lu G F, Zou J, Wang Y, Wang Z. L1-norm-based principal component analysis with adaptive regularization. Pattern Recognition, 2016, 60: 901-907 doi: 10.1016/j.patcog.2016.07.014 [18] Wang H. Block principal component analysis with L1-norm for image analysis. Pattern Recognition Letters, 2012, 33(5): 537-542 doi: 10.1016/j.patrec.2011.11.029 [19] Bing N L, Qiang Y, Rong W, Xiang K, Meng W, Li X. Block principal component analysis with nongreedy L1-norm maximization. IEEE Transactions on Cybernetics, 2016, 46(11): 2543-2547 doi: 10.1109/TCYB.2015.2479645 [20] Li X, Pang Y, Yuan Y. L1-norm-based 2DPCA. IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics), 2010, 40(4): 1170-1175 doi: 10.1109/TSMCB.2009.2035629 [21] Wang R, Nie F, Yang X, Gao F, Yao M. Robust 2DPCA With non-greedy L1-norm maximization for image analysis[J]. IEEE Transactions on Cybernetics, 2015, 45(5): 1108-1112 [22] Wang H, Wang J. 2DPCA with L1-norm for simultaneously robust and sparse modelling. Neural Networks, 2013, 46: 190-198 doi: 10.1016/j.neunet.2013.06.002 [23] Pang Y, Li X, Yuan Y. Robust tensor analysis with L1-norm. IEEE Transactions on Circuits and Systems for Video Technology, 2010, 20(2): 172-178 doi: 10.1109/TCSVT.2009.2020337 [24] Zhao L, Jia W, Wang R, Yu Q. Robust tensor analysis with non-greedy L1-norm maximization. Radioengineering, 2016, 25(1): 200-207 doi: 10.13164/re.2016.0200 [25] Kwak N. Principal component analysis by Lp-norm maximization. IEEE Transactions on Cybernetics, 2014, 44(5): 594-609 doi: 10.1109/TCYB.2013.2262936 [26] 李春娜, 陈伟杰, 邵元海. 鲁棒的稀疏Lp-模主成分分析. 自动化学报, 2017, 43(01): 142-151Li Chun-Na, Chen Wei-Jie, Shao Yuan-Hai. Robust sparse Lp-norm principal component analysis. Acta Automatica Sinica, 2017, 43(01): 142-151 [27] Wang J. Generalized 2-D principal component analysis by Lp-norm for image analysis. IEEE Transactions on Cybernetics, 2016, 46(3): 792-803 doi: 10.1109/TCYB.2015.2416274 [28] Li T, Li M, Gao Q, Xie D. F-norm distance metric based robust 2DPCA and face recognition. Neural Networks, 2017, 94: 204-211 doi: 10.1016/j.neunet.2017.07.011 [29] Wang Q, Gao Q, Gao X, Nie F. Optimal mean two-dimensional principal component analysis with F-norm minimization. Pattern Rrecognition, 2017, 68: 286-294 doi: 10.1016/j.patcog.2017.03.026 [30] Gao Q, Ma L, Liu Y, Gao X, Nie F. Angle 2DPCA: a new formulation for 2DPCA. IEEE Transactions on Cybernetics, 2018, 48(5): 1672-1678 doi: 10.1109/TCYB.2017.2712740 期刊类型引用(8)
1. 葛泉波，程惠茹，张明川，郑瑞娟，朱军龙，吴庆涛. 基于PCA和ICA模式融合的非高斯特征检测识别. 自动化学报. 2024(01): 169-180 . 
本站查看2. 邢楷，单耀，王克南，樊林林，邹阳. 耦合主成分分析下煤矿矿井突水水源识别方法. 能源与环保. 2024(04): 38-43 . 百度学术3. 李伟，雷鹏，黄天尘，张晓利，叶鸥. 基于MLP-DBN模型的构造煤分布预测策略分析. 能源与环保. 2024(04): 118-123 . 百度学术4. 吕弢，陈璟，薛善烨. 离散制造中基于多源多模数据的产品指标联合分析. 空天防御. 2024(05): 110-119+126 . 百度学术5. 毕鹏飞，胡志远，陈璇，杜雪. 双灵活度量自适应加权2DPCA在水下光学图像识别中的应用. 电子与信息学报. 2024(11): 4188-4197 . 百度学术6. 裴正中，赵旭俊. 基于同构化角度的离群检测方法. 计算机工程与设计. 2024(12): 3622-3630 . 百度学术7. 陈璇，毕鹏飞，胡志远. 任意三角形结构2DPCA在水下光学图像识别中的应用. 电子测量与仪器学报. 2024(12): 43-53 . 百度学术8. 吴磊，康英伟. 基于深度强化学习的湿法脱硫系统运行优化. 系统科学与数学. 2022(05): 1067-1087 . 百度学术其他类型引用(5)
-
下载:
