Analysis and System Design of Multi-convex Hull Stabilization Domain for Double-layered Model Predictive Control System
-
摘要: 针对双层模型预测控制（Model predictive control，MPC）中出现的由于系统状态在动态控制（Dynamic control，DC）过程中超出约束集，导致下层优化不可行的问题，本文在综合控制方法的基础上提出一种新的动态控制策略，引入多包镇定域（Multi-convex hull stabilization domain，MHSD）的概念.通过离线计算多包镇定域，并根据系统每一时刻的实测状态值，在线决定（Dynamic control）层的镇定域以及相应的控制时域，结合变约束思想，保证动态控制过程递归可行，从而有效控制在大范围内变化的系统状态.另外，本文通过设计非线性反馈控制器，扩大了终端不变集和多包镇定域的范围，提高了DC层对稳态目标值的跟踪效果.本文的控制算法可以使得DC层在目标跟踪过程中保证递归可行性，并最大程度地实现无静差跟踪.仿真算例验证了本文算法对稳定系统和不稳定系统都有效.Abstract: In order to solve the feasibility problem in double-layered model predictive control (MPC) caused by some states in dynamic control (DC) that violate the constraints, this paper proposes a new control strategy based on the overall control solution and introduces a new definition multi-convex hull stabilization domain (MHSD). This strategy designs the MHSD off line and then chooses a proper stabilization range online according to the real-time system state. The control horizon of the DC layer can be calculated at the same time. What is more, this paper enlarges the invariant sets and the stabilization domain through designing a nonlinear feedback controller so that the system states varying in a wide range can be controlled and the tracking effect is significantly improved. By using the above algorithm, the control process is recursively feasible and the optimal targets can be tracked precisely. The effectiveness of this method is verified in both stable and unstable systems through two examples.
-
工业控制系统(Industrial control system)广泛应用于国家基础设施, 工控系统一旦遭到破坏可能会造成难以估量的经济损失甚至人员伤亡.一个典型的工业控制网路可分为企业网络层、控制网络层和现场网络层, 如图 1所示.
近年来, 由于IT系统的软件和硬件技术不断集成到工控领域, IT领域的一些漏洞和后门也出现在工控系统上.以智能制造为核心的工业4.0, 在推动工业转型的同时也使得工控系统面临更多来自互联网的威胁[1]. 2010年, 世界首个网络超级破坏武器"震网"病毒[2]被检测出来, 证明了工控系统是可被攻击并被利用的.此后又相继爆发了"毒区"、"火焰"等病毒和一连串的工控系统入侵事件.工控信息安全受到越来越多的重视, 针对工控系统入侵检测的研究也成为目前的一个热点.
由于工控保护机制的存在, 大的破坏性攻击容易被检测出来, 因此对入侵检测的研究越来越集中于隐蔽攻击[3-5].当前, 工控系统入侵检测的研究方法主要包含三类: 1)基于概率统计的方法[6-7]; 2)基于机器学习的方法[8-12]; 3)基于状态的方法[13-16].
基于状态的方法因为更好地考虑了工控系统的物理特性, 检测准确度高, 目前大部分研究都是基于此类方法[13].基于状态的方法通过历史数据和专家经验定义系统的临界状态, 实时监控系统当前状态与临界状态的距离以判断系统是否处于危险之中.Fovino等[15]通过监控系统状态的变化来检测复杂攻击, 创立了系统的虚拟镜像作为系统的内部表示, 并用规则语言描述系统的临界状态, 但是仅考虑了离散输入输出数据, 没有考虑连续数据. Carcano等[17]进一步考虑了连续输入和输出, 但在输入输出数据量庞大的情况下, 很难定义系统的临界状态.针对这一问题, Khalili等[13]提出了一个系统化的解决方案SysDetect, 该方案基于Apriori算法, 通过专家经验的判定可显著减少算法下一次迭代产生的候选临界状态数.但是该方法依然不能很好地处理高维数据, 且在一定程度上依赖人工判定.
针对当前工控入侵检测系统存在的问题, 本文提出一种新的基于状态迁移图的异常检测方法.利用工控系统与物理世界交互的特性, 采用数据驱动的方法进行建模, 即利用系统运行时系统变量的数据来建立检测模型.但考虑到数据维度可能过高, 不直接以系统变量数据来描述系统运行状态, 而是以相邻数据向量间的余弦相似度和欧氏距离来表征, 因而可以处理高维数据.状态迁移图刻画系统运行过程中的正常模型, 根据正常历史数据样本训练得出, 所以不需要依赖专家预先定义系统的临界状态.
1. 恶意数据攻击建模
一个工业控制过程可简单地用图 2来表示.其中传感器和执行器易成为攻击者的攻击目标, 因为它们直接与物理过程相连, 一旦遭到攻击将会产生不可估量的后果, 例如"震网"病毒恶意篡改了伊朗核工厂离心机的转速, 使核工厂被迫停工.
图 2 工业控制过程:传感器和执行器易成为攻击目标Fig. 2 Industrial control process: sensors and actuators are vulnerable targets由于工控系统的物理特性, 恶意数据注入成为一种简单且收益高的方式.攻击者通过攻击执行器和传感器, 注入恶意数据, 修改变量的值, 对系统造成破坏.
系统变量分为传感变量和操纵变量, 传感变量的值由传感器测量取得, 而操纵变量则由控制器通过系统参数和传感变量值计算得出.记系统含$l$个测量变量, 用集合${\mathit{\boldsymbol{y}}} = \{ {{y_k}|k = 1, 2, \cdots , l} \}$表示, 含$p$个操纵变量, 用集合${\mathit{\boldsymbol{u}}} = \{ {{u_k}|k = 1, 2, \cdots , p} \}$表示.则系统变量集合可用${\mathit{\boldsymbol{s}}} = \{ {{s_j}|j = 1, 2, \cdots , l + p} \}$表示.根据攻击者试图破坏系统速度的快慢, 将恶意数据注入攻击分为快速注入和慢速注入.
快速注入使得系统变量在短时间内发生较大改变, 以最快速度造成破坏, 表现为最大/最小值攻击.
$ {\tilde s_i}\left( t \right) = \begin{cases} {s_i}\left( t \right)\text{, }{\kern 20pt} t \notin {T_a}\\ {s_i}^{\max }\text{, }{\kern 20pt} t \in {T_a} \end{cases} $
(1) $ {\tilde s_i}\left( t \right) = \begin{cases} {s_i}\left( t \right)\text{, }{\kern 20pt} t \notin {T_a}\\ {s_i}^{\min }\text{, }{\kern 22pt} t \in {T_a} \end{cases} $
(2) 其中, ${T_a}$为攻击时段.
慢速注入缓慢改变系统变量的值, 以达到潜伏的目的, 可以表现为偏置注入和几何注入.偏置注入连续注入多个较小的常量.
$ {\tilde s_i}\left( t \right) = \begin{cases} {s_i}\left( t \right) + {c_i}\text{, }{\kern 20pt} t \in {T_a}\\ {s_i}\left( t \right)\text{, } t \notin {T_a} \end{cases} $
(3) 其中, 表示常量, 通常取值较小.几何注入逐渐增大系统变量的改变幅度, 具有一定的潜伏性, 同时不失破坏性, 可表现为指数形式的注入.
$ {\tilde s_i}\left( t \right) = \begin{cases} {s_i}\left( t \right) + a{b^{t - {t_0}}}\text{, }{\kern 10pt} t \in {T_a}\\ {s_i}\left( t \right)\text{, }{\kern 49pt} t \notin {T_a} \end{cases} $
(4) 其中, $a$, $b$可以是常数, 也可以是变量, ${t_0}$为常数.
以上只是列举了三种可能的攻击形式, 实际上恶意数据注入攻击可以表现为更多的形式.
2. 基于状态迁移图的异常检测方法
由于工控系统的诸多限制, 例如有限的内存、有限的计算能力、较高的实时性要求等, 工控入侵检测系统必须是轻量级的, 检测规则或检测模型应设置得相对简单.本文采用一个状态迁移图的检测模型, 状态迁移图的状态数与数据向量的维度无关, 因而可以处理高维的数据向量.
2.1 状态表示
实际的工控系统中, 状态变量可能有很多个, 因此数据向量的维度可能很高, 而现有的基于状态的入侵检测算法无法很好地处理数据维度较高的情况.考虑控制系统的运行模式, 控制器根据前一时刻传感变量的值, 计算出操纵变量的值, 操纵变量的值再作用于物理系统使之发生变化, 传感变量随着系统的运动而变化, 在下一个采样周期将传感变量的值传给控制器, 如此循环往复.通过分析这个过程可以发现, 系统变量在当前时刻的取值很大程度上取决于上一时刻的取值, 因此正常情况下, 相邻两个数据向量间的变化不会太大.
相邻数据向量之间的变化反应了系统变量运行的时间特性, 本文以余弦相似度和欧氏距离来衡量这种变化.余弦相似度反应数据向量方向的变化, 常用于文档聚类和信息提取.欧氏距离衡量两个数据向量绝对距离的远近.
选择作为系统状态表征的特征应该是相互关联的, 即相邻数据向量间的余弦相似度和欧氏距离应该是相互关联的.这样当余弦相似度在某个区间取值时, 欧氏距离只位于某些特定区间, 而不会在所有区间取值.如图 3所示, 假设余弦相似度取值区间为$A$, $B$, $C$, 欧氏距离取值区间为$D$, $E$, $F$, 二者在二维空间所有可能的区域为1 $\sim$ 9, 实际取值区域为1, 4, 5, 9.当余弦相似度取值区间为$B$时, 欧氏距离只在区间$D, E$取值, 而不会在$F$区间取值, 若欧氏距离在$F$区间取值就可以认为发生了异常.这样可以更加细致地刻画数据的轮廓, 而不是将所有可能的区域为1 $\sim$ 9都视为可接受的区域.余弦相似度和欧氏距离的关联性如图 4所示, 其中余弦相似度和欧氏距离采用Normal数据集计算得出(经过归一化处理), 图 4展现了前30组数据.从图 4可以看出, 余弦相似度和欧氏距离的变化基本相反, 当余弦相似度在高位取值时, 欧氏距离总在低位取值, 证明了系统变量数据向量间的余弦相似度和欧氏距离是相互关联的.
图 3 余弦相似度和欧氏距离取值示意Fig. 3 Possible taking-value for cosine similarity and Euclidian distance本文以数据向量间的余弦相似度和欧氏距离组成的二元组来表征系统运行状态, 避免了直接对数据进行建模, 将高维的数据向量转化为二维数据, 有效降低了计算复杂度.同时, 考虑数据的时间特性, 更加符合工控系统的动态特性.
2.2 状态矩阵求解
为了求解出状态迁移图, 需要知道图的顶点和各顶点之间的转换关系.状态迁移图的顶点表示了系统可能的状态, 本文用一个状态矩阵求解状态迁移图, 需要知道图的顶点和各顶点之间的转换关系.状态迁移图的顶点表示了系统可能的状态, 本文以状态矩阵$STATE$表示可能的状态, $STATE$的每一行称为一个状态点, 一个状态点表示一个状态, 对应状态迁移图中的一个顶点.
2.2.1 余弦相似度和欧氏距离的求解
在某一时刻所有系统变量取值的集合称为一样本点, 例如${{\mathit{\boldsymbol{m}}}_i} = \{ {{s_{1, i}}, {s_{2, i}}, \cdots , {s_{l + p, i}}} \}$表示系统在$i$时刻的样本点, ${s_{k, t}} = {s_k}( t )$为第$k$个系统变量在$t$时刻的取值.为了计算方便, 用行向量的形式来表示样本点, 即${{\mathit{\boldsymbol{m}}}_i} = [ {{s_{1, i}}, {s_{2, i}}, \cdots , {s_{l + p, i}}} ]$.假设训练样本点数为$n$, 相邻样本点(经过z-score标准化的样本点)间的余弦相似度由式(5)计算得出.
$ cossi{m_i} = \frac{{{{\mathit{\boldsymbol{m}}}_i} \times {{\mathit{\boldsymbol{m}}}}^{\rm T}_{i + 1}}}{{\left\| {{{\mathit{\boldsymbol{m}}}_i}} \right\| \times \left\| {{{\mathit{\boldsymbol{m}}}_{i + 1}}} \right\|}}, \quad i = 1, 2, \cdots , n - 1 $
(5) 令
${\mathit{\boldsymbol{cossim}}}$是所有训练样本点任意相邻两点间余弦相似度组成的列向量.相邻样本点间的欧氏距离由式(6)计算得出.
$ e{d_i} = \sqrt {{{{{\mathit{\boldsymbol{m}}}}}_i} \times {{\mathit{\boldsymbol{m}}}}^{\rm T}_{i + 1}} , \quad i = 1, 2, \cdots , n{\rm{ - }}1 $
(6) 令${\mathit{\boldsymbol{ed}}} = {[ {e{d_1}, e{d_2}, \cdots , e{d_{n - 1}}} ]^{\rm{T}}}$, ${\mathit{\boldsymbol{ed}}}$是所有训练样本点任意相邻两点间欧氏距离组成的列向量.
2.2.2 cossim和ed量化
由于${\mathit{\boldsymbol{cossim}}}$和${\mathit{\boldsymbol{ed}}}$中的元素都是实数, 这样通过二者确定的状态图中的状态数目将是无限个, 所以需要对它们进行量化处理, 使得${\mathit{\boldsymbol{cossim}}}$和${\mathit{\boldsymbol{ed}}}$中的元素只能取一些固定值.可以用一些量化阈值将${\mathit{\boldsymbol{cossim}}}$和${\mathit{\boldsymbol{ed}}}$的取值分成数个量化区间, 同一量化区间的点用相同的值来替换.本文用$intervals$表示${\mathit{\boldsymbol{cossim}}}$和${\mathit{\boldsymbol{ed}}}$量化区间的数量.
为便于处理, 首先对${\mathit{\boldsymbol{cossim}}}$和${\mathit{\boldsymbol{ed}}}$进行归一化处理.将${\mathit{\boldsymbol{cossim}}}$归一化到$[-1, 1]$区间, 得到$\pmb{comsim}'$; ${\mathit{\boldsymbol{ed}}}$归一化到$[0, 1]$区间, 得到${\mathit{\boldsymbol{ed}}}'$.令$\pmb{comsim}'$的量化阈值集为
$ {\mathit{\boldsymbol{r}}} = \left\{ {{r_0}, {r_1}, {r_2}, \cdots , {r_{\rm intervals}}} \right\} $
(7) ${\mathit{\boldsymbol{ed}}}'$的量化阈值集为
$ {\mathit{\boldsymbol{d}}} = \left\{ {{d_0}, {d_1}, {d_2}, \cdots , {d_{\rm intervals}}} \right\} $
(8) 按照各量化区间点数相等的原则计算${\mathit{\boldsymbol{r}}}$和${\mathit{\boldsymbol{d}}}$.定义函数$k:{\boldsymbol{\rm{R}}} \to {\boldsymbol{\rm{R}}}$, 表示某实数集合${{P}}$中元素的个数.定义函数
$ f(x) = k(\{ j|j \le x, j \in {\mathit{\boldsymbol{cossim}}}'\} ) $
$ g(x) = k(\{ j|j \le x, j \in {\mathit{\boldsymbol{ed}}}'\} ) $
${\mathit{\boldsymbol{r}}}$和${\mathit{\boldsymbol{d}}}$的计算如下:
$ {r_k} = \begin{cases} - 1\text{, }\quad\ k = 0\\ \arg\min \limits_{x \in \left[ { -1, 1} \right]} f\left( x \right) - f\left( {{r_{k - 1}}} \right) = \dfrac{n}{{intervals}}\text{, }\\{\kern 54pt}\quad\ \, k = 1, 2, \cdots , intervals - 1\\ 1\text{, }\quad\ \, \;\;\;\;\;\; k = intervals \end{cases} $
(9) $ {d_k} = \begin{cases} 0\text{, }\quad\ k = 0\\ \arg\min\limits_{x \in \left[ { 0, 1} \right]} g\left( x \right) - g\left( {{d_{k-1}}} \right) = \dfrac{n}{{intervals}}\text{, }\\{\kern 54pt}\quad\ \, k = 1, 2, \cdots , intervals - 1\\ 1\text{, }\quad\ \;\;\;\;\;\; k = intervals \end{cases} $
(10) 2.2.3 状态矩阵生成
确定量化阈值后, 根据量化阈值对$\pmb{comsim}'$和${\mathit{\boldsymbol{ed}}}'$进行量化处理.
$ \begin{align} & qcossi{m_i} = \arg\min\limits_{r \in {{\mathit{\boldsymbol{r}}}}} {\kern 1pt} cossi{m_i} \le r, \nonumber\\&\qquad\qquad\qquad\qquad\qquad\quad i = 1, 2, \cdots , n - 1 \end{align} $
(11) $\pmb{qcomsim}'\!\!=\!\!{[ {qcossi{m_1}|qcossi{m_2}| \cdots |qcossi{m_{n-1}}} ]^{\rm{T}}}$为量化后的$\pmb{comsim}'$.
$ qe{d_i} = \arg\min\limits_{d \in {\mathit{\boldsymbol{d}}}} e{d_i} \le d, \quad i = 1, 2, \cdots , n - 1 $
(12) ${\mathit{\boldsymbol{qed}}}' = [ {qe{d_1}|qe{d_2}| \cdots |qe{d_{n - 1}}} ]$为量化后的${\mathit{\boldsymbol{ed}}}'$.最终的状态矩阵为$STATE = [ {\pmb{qcossim}'|{\mathit{\boldsymbol{qed}}}'} ]$.
2.3 状态迁移图生成
状态迁移图是一个有向图, 作为异常检测模型, 其顶点由$STATE$确定, 边用邻接矩阵表示.
2.3.1 顶点确定
由于对相邻样本点间的余弦相似度和欧氏距离进行了量化, 状态矩阵$STATE$中含有相同的状态点.相同的状态点对应状态迁移图中的同一个顶点, 应将其合并, 具体做法为:
步骤 1. 用自然数对$STATE$中的每个状态点进行标号, 称为状态号, 并用向量${\mathit{\boldsymbol{sign}}}$来表示, 例如${\mathit{\boldsymbol{sign}}}\left( k \right)$表示$STATE$中的第$k$个状态点对应的状态号.状态号唯一标明了状态迁移图中的状态(即顶点).
步骤 2. 令$STATE$第1个状态点的状态号为1, 即${\mathit{\boldsymbol{sign}}}(1) = 1$, 令状态数$sn = 1$.
步骤 3. 从第2个状态点开始, 对第$k$个状态点遍历$STATE$中该状态点之前的状态点, 将该状态点与之前的状态点进行比较, 若找到相同的状态点$j$则停止遍历, 将该点的状态号设为第$j$个状态点的状态号, 即${\mathit{\boldsymbol{sign}}}(k) = {\mathit{\boldsymbol{sign}}}(j)$; 若没有找到相同的状态点, 则令, 将该点的状态号设为$sn$, 即${\mathit{\boldsymbol{sign}}}(k) = sn$.
经过步骤1 $\sim$ 3后, $STATE$中相同的状态点具有相同的状态号, 最终的状态数为$sn$.状态迁移图的顶点对应$STATE$中状态号为1 $\sim$ $sn$的状态点.
2.3.2 邻接矩阵确定
用$NM$表示邻接矩阵, 邻接矩阵反映了状态迁移图中的状态转移关系, 若${NM_{i, j}} = 1$, 则状态迁移图中的第$j$个顶点到第$j$个顶点可达; 若${NM_{i, j}}$ $=$ $0$, 则第$i$个顶点到第$i$个顶点不可达.由于样本点是随时间增长逐渐采样而来的, 因此前一样本点到当前样本点总是可达的.相应地, $STATE$中前一状态点到当前状态点也总是可达的.根据这一原则可以求解出邻接矩阵$NM$, 具体做法为:
步骤 1. 将$NM$置为全零.
步骤 2. 根据前一样本点(状态点)到当前样本点(状态点)可达的原则, 令$NM({\mathit{\boldsymbol{sign}}}( i )$, ${\mathit{\boldsymbol{sign}}}( i$ $+$ $1 ))$ $= 1$, 由此计算出$NM$.
确定了顶点和邻接矩阵后, 就能生成状态迁移图.
2.4 基于状态迁移图的在线检测
状态迁移图生成后, 可以据此进行在线异常检测.基于状态迁移图的在线检测流程如图所示, 具体步骤如下:
步骤 1. 获取当前样本点.
步骤 2. 计算当前样本点与上一样本的余弦相似度$t\_cossi{m_i}$和欧氏距离$t\_e{d_i}$.
步骤 3. 按照训练样本的标准对$t\_cossi{m_i}$和$t\_e{d_i}$进行归一化处理, 得到$qt\_cossi{m_i}$和$qt\_e{d_i}$.
步骤 4. 遍历$STATE$, 判断$STATE$中是否有和当前点相同的点.若存在, 根据${\mathit{\boldsymbol{sign}}}$确定当前状态点对应的状态号.
步骤 5. 若当前状态点对应的状态不在状态迁移图内, 产生告警, 记录异常类别为"异常1";若当前点对应的状态位于状态图内, 即$STATE$中有和当前点相同的点, 则根据状态迁移图判断上一节点对应状态到当前节点对应状态是否可达.若不可达, 则产生告警, 记录异常类别为"异常2";若可达, 继续下一个样本点的检测.
3. 实验与测试
目前工控领域没有一个标准的测试集, 而实际工控系统在攻击下的数据样本很难获取.针对这一情况, 本文采用田纳西-伊斯曼(Tennessee-eastman, TE)过程[18]的MATLAB模型作为仿真平台.根据第1节建立的恶意数据攻击模型, 选定斜坡注入和偏置注入两种攻击形式, 在特定时刻对TE过程实施攻击, 生成系统正常运行数据和攻击情形下的数据, 并用第2.4节描述的检测方法进行检测.
3.1 TE仿真模型
TE过程的原型由Downs和Vogel在1993年根据一个真实的化工过程创建. TE过程主要部件有冷凝器、反应器、循环压缩机、汽提塔和气液分离器等.该过程包含两个放热反应和两个副反应, 生成产品$G$, $H$和副产品$F.$四种化学反应如下所示:
$ {{A}}\left( { {g}} \right) + {{C}}\left( { {g}} \right) + {{D}}\left( { {g}} \right) \to {{G}}\left( {{ {liq}}} \right) $
(13) $ {{A}}\left( { {g}} \right) + {{C}}\left( { {g}} \right) + {{E}}\left( { {g}} \right) \to {{H}}\left( {{ {liq}}} \right) $
(14) $ {{A}}\left( { {g}} \right) + {{E}}\left( { {g}} \right) \to {{F}}\left( {{ {liq}}} \right) $
(15) $ {{D}}\left( { {g}} \right) \to 2{{F}}\left( {{ {liq}}} \right) $
(16) TE过程包含41个测量变量和12个操纵变量, 是一个多变量, 多数据, 复杂的非线性控制系统, 常用于多变量控制、非线性控制和故障诊断等领域. TE过程的MATALB模型由Ricker[19]给出, 该模型采用文献[20]中的控制策略, 仿真时间为48 h, 采样时间为3 min, 每次仿真共产生960组数据, 每组数据中都包含高斯白噪声. TE过程的工艺流程如图 6所示.
3.2 测试数据
反应器是TE过程的核心部件, 是发生化学反应的地方, 对温度的要求很高.因此本文选定攻击对象为反应器的温度传感器, 并根据第1节所述攻击模型, 选定偏置注入和斜坡注入两种攻击方式.
$ {\tilde s_i}\left( t \right) = \begin{cases} {s_i}\left( t \right) + c\text{, }& t \in {T_a}\\ {s_i}\left( t \right)\text{, }& t \notin {T_a} \end{cases} $
(17) $ {\tilde s_i}\left( t \right) = \begin{cases} {s_i}\left( t \right) + k(t - {t_0})\text{, }& t \in {T_a}\\ {s_i}\left( t \right)\text{, }& t \notin {T_a} \end{cases} $
(18) 其中, 偏置注入中$c$为常数, 而斜坡注入在${t_0}$时注入斜率为$k$的斜坡信号.反应器温度控制回路如图 7所示, 两个输入分别代表反应器温度初始设定值(用输入点9表示)和反应器温度(用xmeas 9表示), 输出为反应器冷却水流量(用xmv10表示).该回路通过反应器温度测量值和反应器温度初始设定值计算冷却水流量, 以此来控制反应器温度.为达到攻击反应器温度的目的, 在反应器温度输入端口增加Add模块, 将其与攻击信号叠加, 攻击信号通过延时模块在特定时刻发挥作用.图 8为加入攻击信号后的反应器温度控制回路, 延时模块时间设置为20 h, 攻击信号为常数, 代表偏置注入信号.攻击信号可以替换为其他信号, 以此设计更多的攻击形式.
图 8 加入攻击信号的反应器温度控制回路Fig. 8 Control loop for reactor temperature added with attack signal采集系统正常运行时各变量的数据, 生成数据集合Normal.选定式(17)中$c$为0.01, 0.1和1, 在系统仿真20 h加入注入攻击, 生成偏置注入下的数据集Dataset1, Dataset2, Dataset3.选定式(18)中$k$为0.01, 0.1和1, 在仿真20 h实施斜坡攻击, 生成斜坡注入下的数据集Dataset4, Dataset5, Dataset6.各数据集及其对应参数见表 1.这些数据集将用来测试本文所提的方法.
表 1 测试数据集及其参数Table 1 The test data set and the corresponding parameters实验数据 $c$ $k$ Normal 0 0 Dataset1 0.01 0 Dataset2 0.1 0 Dataset3 1 0 Dataset4 0 0.01 Dataset5 0 0.1 Dataset6 0 1 为直观感受加入攻击信号后对反应器温度的影响, 用MATLAB采集生成了反应器温度在各种情况下的温度-时间变化曲线, 分别如图 9 $\sim$ 11所示.图 9为正常工况条件下的温度-时间变化曲线, 可以看出反应器温度稳定在123℃左右. 图 10为斜坡注入$k=0.01$情况下的温度变化情况, 可以看到在20 h后反应器温度呈线性下降趋势. 图 11为偏置注入$c=0.1$情况下的温度-时间变化曲线, 反应器温度在20 h突然下降, 之后则稳定下来.
图 9 正常工况下反应器温度随时间变化情况Fig. 9 Reactor temperature varies with time under normal condition
图 10 斜坡注入工况($k=0.01$)下, 反应器温度随时间变化情况Fig. 10 Reactor temperature varies with time under ramp signal injection with $k$ set at 0.01
图 11 偏置注入工况($c=0.1$)下, 反应器温度随时间变化情况Fig. 11 Reactor temperature varies with time under bias signal injection with $c$ set at 0.13.3 状态迁移图生成测试
由于状态迁移图通过正常数据训练得出, 所以用表 1中的Normal数据集来进行训练.根据第2.2节的分析, 生成的状态迁移图与余弦相似度和欧氏距离的量化区间数$intervals$相关, 因此本文不断改变$intervals$的取值, 进行了多组状态图测试.首先按照第2.2节和第2.3节的方式对Normal数据集进行处理, 然后设定量化区间数$intervals$, 计算状态矩阵, 最后生成状态迁移图.状态迁移图的顶点数和边数随$intervals$的变化情况如表 2所示.
表 2 状态迁移图的顶点数和边数随$intervals$的变化Table 2 The nodes and triangles number of state transition graph varies with $intervals$$intervals$ 顶点数 边数 3 9 48 5 20 147 8 44 315 10 73 677 12 97 782 15 143 871 从表 2可以看出, 随着$intervals$的增大, 状态迁移图的顶点数和边数也不断增大.而$intervals$越大, 意味着将余弦相似度和欧氏距离的量化区间划分得越小, 由于状态迁移图的检测由顶点和边来决定, 所以顶点数和边数越多, 意味着检测规则的粒度越细, 相应地对异常会更加敏感, 此外还会额外增加资源开销. 图 12和图 13分别显示了$intervals$为5和8时用MATLAB生成的状态迁移图.
图 12 $intervals$为5时的状态迁移图Fig. 12 The state transition graph when $intervals$ is equal to 5
图 13 $intervals$为8时的状态迁移图Fig. 13 The state transition graph when $intervals$ is equal to 83.4 异常检测性能测试
按照第2.4节的在线检测步骤, 用表 1的测试数据对状态迁移图检测模型的检测效果进行测试, 选定$intervals$为5, 8, 10, 15, 检测结果分别见表 3 $\sim$ 6.本文选用从攻击到正确检测到异常的时间以及误报率作为评价的标准.
表 3 $intervals=5$时的检测结果Table 3 Detection results when $intervals$ is equal to 5Dataset1 Dataset2 Dataset3 Dataset4 Dataset5 ataset6 正确检测到异常的样本点数 409 401 401 436 407 402 异常类别 异常2 异常2 异常2 异常2 异常2 异常2 误报率(%) 0.63 0.21 0.42 0.21 0.83 0.42 表 4 $intervals=8$时的检测结果Table 4 Detection results when $intervals$ is equal to 8Dataset1 Dataset2 Dataset3 Dataset4 Dataset5 Dataset6 正确检测到异常的样本点数 401 401 401 411 401 401 异常类别 异常2 异常2 异常2 异常2 异常2 异常2 误报率(%) 5.62 4.38 5.21 4.17 5.83 5.42 表 5 $intervals=10$时的检测结果Table 5 Detection results when $intervals$ is equal to 10Dataset1 Dataset2 Dataset3 Dataset4 Dataset5 Dataset6 正确检测到异常的样本点数 401 401 401 401 401 401 异常类别 异常2 异常2 异常2 异常2 异常2 异常2 误报率(%) 7.50 7.08 8.12 8.21 6.67 7.71 表 6 $intervals=15$时的检测结果Table 6 Detection results when $intervals$ is equal to 15Dataset1 Dataset2 Dataset3 Dataset4 Dataset5 Dataset6 正确检测到异常的样本点数 401 401 401 401 401 401 异常类别 异常2 异常1 异常1 异常1 异常1 异常1 误报率(%) 24.37 26.25 22.92 30.21 28.96 26.04 从检测时间来看, 在各种攻击情况下检测模型检测到异常的时间随着$intervals$的变化而变化.从检测最坏的结果对应的数据集来看(对应的数据集为Dataset4, 此时$c = 0$, $k = 0.01$), 当$intervals$为5和8时候, 分别在第436和第411个样本点检测到异常(在第20小时进行注入攻击, 对应的第一个异常样本点为第401个).而当$intervals$为10和15时, 检测模型在第401个样本点即检测到异常, 说明当$intervals$增大时, 检测模型对异常更加敏感, 检测异常的速度有所提升.
从误报率来看, 当$intervals$逐渐增大时, 相应地误报率也随之增大.当$intervals$为5时误报率均不超过1 %.当$intervals$增大到8时, 误报率迅速增加到5 %左右.而当$intervals$取15时, 误报率超过20 %, 这时可认为状态迁移图不能正确检测异常.
3.5 主元分析法测试
由于工控领域缺乏标准的数据集, 研究人员都用各自的模型和数据进行工控系统异常的研究, 很难复现, 难以对比各种方法的好坏.为了进一步评估状态迁移图检测的性能, 将其与控制系统常用的异常检测方法---主元分析法(Principal component analysis, PCA)进行比较. PCA采用${{\rm{T}}^2}$统计量和$\rm{SPE}$统计量作为检测异常的指标, 当样本的${{\rm{T}}^2}$统计量和$\rm{SPE}$统计量超过各自的控制限时, 检测到异常.选定状态迁移图模型检测结果最差的数据集Dataset4作为PCA方法的测试集, 其测试结果如图 14和图 15所示.从图 14和图 15可以看出, PCA虽然最终能检测到异常, 但是在第500个样本点才检测出来.相比于状态迁移图检测模型, 检测到异常的速度较慢.
3.6 讨论
本节从时空资源消耗以及检测性能两方面对状态迁移图检测模型进行讨论, 由于训练过程是通过历史数据离线进行的, 因此时空消耗只考虑检测时的情况.
1) 从时间消耗上来看, 对新样本点的检测需要判断新样本点是否位于状态迁移图内, 即需要遍历一次状态图, 若状态迁移图含$n$个顶点, 则时间复杂度为O$(n)$.实际检测过程中, 顶点数与设置的量化区间数$intervals$有关.通过第3.4节的实验结果不难发现, 当$intervals$取5或8时已有较好的检测结果, 继续增大$intervals$只会增大资源消耗, 同时增加误报率.且$intervals$越大, 需要越多的训练样本进行训练.以$intervals$取10为例, 顶点数$n$最多为100, 而实际会小于100, 因此这一步骤需要的计算量很小.第二个步骤判断上一状态到当前状态是否可达, 只需查询邻接矩阵中对应的数值即可, 时间复杂度为O$(1)$.
2) 从空间消耗上来看, 本文的异常检测方法需要存储量化阈值集合${\mathit{\boldsymbol{r}}}$和${\mathit{\boldsymbol{d}}}$, 状态矩阵$STATE$, 邻接矩阵$NM$, 其中邻接矩阵$NM$占用的存储空间为$n$的二阶次, ${\mathit{\boldsymbol{r}}}$和${\mathit{\boldsymbol{d}}}$, 状态矩阵$STATE$存储消耗为$n$的一阶次.但由于$n$取值很小(不超过100), 且邻接矩阵存储单位为比特, 当$intervals$取10时, 估计整体消耗的存储空间不超过1 MB.
3) 从检测性能来看, 即便对于微小的攻击(对应数据集Dataset1和Dataset4), 检测模型也能较快检测出异常, 相比于常规的异常检测算法, 例如第3.4节的PCA算法, 状态迁移图检测模型能够更好地刻画系统运行的动态变化, 更快地检测出异常.但是当$intervals$较大(取10以上)时, 检测误报率较高.而当$intervals$较小(例如, 取5时)时, 不能快速地对微小攻击做出反应.因此需要谨慎选择$intervals$以平衡误报率和检测速度之间的关系.
4. 总结
本文提出了一种新的基于状态迁移图的异常检测方法, 详细描述了其状态表示、状态图生成和在线检测过程, 建立了工控系统的恶意数据攻击模型, 并基于TE仿真模型进行了各项实验, 最后从时空资源消耗和检测性能两方面对方法进行了讨论.该方法具有较小的时空消耗, 且考虑了系统运行的动态特性, 能够快速检测出异常, 因而适用于资源受限和对实时性要求较高的工业控制系统.
-
图 3 控制时刻$k$从31到60对应的纸机系统控制过程
Fig. 3 Control process of the paper system with the control moment $k$ from 31 to 60
图 4 控制时刻$k$从61到90对应的纸机系统控制过程
Fig. 4 Control process of the paper system with the control moment $k$ from 61 to 90
图 5 控制时刻$k$从31到60对应的双积分系统控制过程
Fig. 5 Control process of the double-integrator system with the control moment $k$ from 31 to 60
表 1 本文符号及其含义
Table 1 The meanings of the notations in this paper
符号 含义 $x^*$ $x$的最优值 ${\bf R}^n$ $n$维欧氏空间 $k$ 离散采样间隔 $x$ 系统状态, $x \in {\bf R}^{n_x}$ $u$ 系统输入, $u \in {\bf R}^{n_u}$ $x_s(u_s)$ 稳态状态(输入) $x_t(u_t)$ 期望稳态状态(输入) $\bar{x}(\bar{u})$ 状态(输入)上界 $I_n $ $n$维单位矩阵 $Q_s, R_s$ 适维权重矩阵 $N_i$ 第$i$个镇定域所对应的控制时域 ${\| x\|}_{Q_s}^2$ $x^{\rm T}$$Q_s$$x$ $x(k+i|k)$ $k$时刻对未来状态的预测值 $u(k+i|k) $ $k$时刻对未来输入的预测值 
下载: 导出CSV
表 2 纸机系统的稳态目标计算结果
Table 2 The results of the SSTC in the AS DPS system
$k$ $u_{s, 1}$ $u_{s, 2}$ $x_{s, 1}$ $x_{s, 2}$ $1\sim90$ -0.39 -0.41 -0.32 -0.33
下载: 导出CSV
表 3 双积分器系统的稳态目标计算结果
Table 3 The results of SSTC in the double-integrator system
k $u_t$ $x_t$ $u_s$ $x_s$ $1\sim30$ (0, 0) (2, -2) (0, 0) (2, 0.5) $31\sim60$ (0, 0) (0, -2) (-0.1, 0.2) (-0.38, 0.2998) $61\sim90$ (0, 0) (0, 0) (0, 0) (0, 0)
下载: 导出CSV
-
[1] 席裕庚, 李德伟.预测控制定性综合理论的基本思想和研究现状.自动化学报, 2008, 34 (10):1225-1234 http://www.aas.net.cn/CN/abstract/abstract17992.shtmlXi Yu-Geng, Li De-Wei. Fundamental philosophy and status of qualitative synthesis of model predictive control. Acta Automatica Sinica, 2008, 34 (10):1225-1234 http://www.aas.net.cn/CN/abstract/abstract17992.shtml [2] 席裕庚, 李德伟, 林姝.模型预测控制-现状与挑战.自动化学报, 2013, 39 (3):222-236 http://www.aas.net.cn/CN/abstract/abstract17874.shtmlXi Yu-Geng, Li De-Wei, Lin Shu. Model predictive control-status and challenges. Acta Automatica Sinica, 2013, 39 (3):222-236 http://www.aas.net.cn/CN/abstract/abstract17874.shtml [3] 席裕庚.预测控制.第2版.北京:国防工业出版社, 2013. 15-67Xi Yu-Geng. Predictive Control (Second Edition). Beijing:National Defense Industry Press, 2013. 15-67 [4] Sildir H, Arkun Y, Cakal B, Gokce D, Kuzu E. Plant-wide hierarchical optimization and control of an industrial hydrocracking process. Journal of Process Control, 2013, 23 (9):1229-1240 doi: 10.1016/j.jprocont.2013.07.007 [5] Darby M L, Nikolaou M, Jones J, Nicholson D. RTO:an overview and assessment of current practice. Journal of Process Control, 2011, 21 (6):874-884 doi: 10.1016/j.jprocont.2011.03.009 [6] Scattolini R. Architectures for distributed and hierarchical model predictive control-a review. Journal of Process Control, 2009, 19 (5):723-731 doi: 10.1016/j.jprocont.2009.02.003 [7] 李世卿, 丁宝苍.基于动态矩阵控制的双层结构预测控制的整体解决方案.自动化学报, 2015, 41 (11):1857-1866 http://www.aas.net.cn/CN/abstract/abstract18761.shtmlLi Shi-Qing, Ding Bao-Cang. An overall solution to double-layered model predictive control based on dynamic matrix control. Acta Automatica Sinica, 2015, 41 (11):1857-1866 http://www.aas.net.cn/CN/abstract/abstract18761.shtml [8] 李世卿, 丁宝苍, 孙耀.双层预测控制中基于操作变量增量的多优先级稳态目标计算.控制理论与应用, 2015, 32 (2):239-245 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=kzllyyy201502014Li Shi-Qing, Ding Bao-Cang, Sun Yao. Multi-priority rank steady-state target calculation in double-layered model predictive control by optimizing increments of manipulated variables. Control Theory & Applications, 2015, 32 (2):239-245 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=kzllyyy201502014 [9] Mayne D Q, Rawlings J B, Rao C V, Scokaert P O M. Constrained model predictive control:stability and optimality. Automatica, 2000, 36 (6):789-814 doi: 10.1016/S0005-1098(99)00214-9 [10] Gutman P O, Cwikel M. An algorithm to find maximal state constraint sets for discrete-time linear dynamical systems with bounded controls and states. IEEE Transactions on Automatic Control, 1987, 32 (3):251-254 doi: 10.1109/TAC.1987.1104567 [11] Gilbert E G, Tan K T. Linear systems with state and control constraints:the theory and application of maximal output admissible sets. IEEE Transactions on Automatic Control, 1991, 36 (9):1008-1020 doi: 10.1109/9.83532 [12] 潘红光, 高海南, 孙耀, 张英, 丁宝苍.基于多优先级稳态优化的双层结构预测控制算法及软件实现.自动化学报, 2014, 40 (3):405-414 http://www.aas.net.cn/CN/abstract/abstract18305.shtmlPan Hong-Guang, Gao Hai-Nan, Sun Yao, Zhang Ying, Ding Bao-Cang. The algorithm and software implementation for double-layered model predictive control based on multi-priority rank steady-state optimization. Acta Automatica Sinica, 2014, 40 (3):405-414 http://www.aas.net.cn/CN/abstract/abstract18305.shtml [13] Blanchini F. Set invariance in control. Automatica, 1999, 35 (11):1747-1767 doi: 10.1016/S0005-1098(99)00113-2 [14] Li Z J, Tan W, Nian S C, Liu J Z. A stabilizing model predictive control for linear systems with input saturation. In: Proceedings of the 2006 International Conference on Machine Learning and Cybernetics. Dalian, China: IEEE, 2006. 671-675 [15] Li Z J, Liu J Z, Tan W. Multi-model H_∞ loop shaping controller for nonlinear system based on gap metric. In: Proceedings of the ICARCV 8th Control, Automation, Robotics and Vision Conference. Kunming, China: IEEE, 2004, 3: 1940-1944 [16] Kosut R. Design of linear systems with saturating linear control and bounded states. IEEE Transactions on Automatic Control, 1983, 28 (1):121-124 doi: 10.1109/TAC.1983.1103127 [17] 李志军. 约束模型预测控制的稳定性与鲁棒性研究[博士学位论文], 华北电力大学(北京), 中国, 2005.Li Zhi-Jun. Research on Stability and Robustness of Constrained Model Predictive Control[Ph. D. dissertation], North China Electric Power University (Beijing), China, 2005. [18] Limon D, Alvarado I, Alamo T, Camacho E F. MPC for tracking piecewise constant references for constrained linear systems. Automatica, 2008, 44 (9):2382-2387 doi: 10.1016/j.automatica.2008.01.023 [19] Chen H, Allgöwer F. A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability. Automatica, 1998, 34 (10):1205-1217 doi: 10.1016/S0005-1098(98)00073-9 [20] 曹永岩, 毛维杰, 孙优贤, 冯旭.现代控制理论的工程应用.杭州:浙江大学出版社, 2000.Cao Yong-Yan, Mao Wei-Jie, Sun You-Xian, Feng Xu. Engineering Application of Modern Control Theory. Hangzhou:Zhejiang University Press, 2000. [21] Pluymers B, Rossiter J A, Suykens J A K, Moor B D. The efficient computation of polyhedral invariant sets for linear systems with polytopic uncertainty. In: Proceedings of the 2005 American Control Conference. Portland, OR, USA: IEEE, 2005, 2: 804-809 期刊类型引用(20)
1. 巩彬，安爱民，石耀科，杜先君. 基于IMODA自适应深度信念网络的复杂模拟电路故障诊断方法. 电子科技大学学报. 2024(03): 327-344 . 
百度学术2. 金亮，闫银刚，杨庆新，刘素贞，张闯. 小样本条件下永磁同步电机深度迁移学习性能预测方法. 电工技术学报. 2023(18): 4921-4931 . 百度学术3. 朱小勇，陈胜. 基于ResNet-ViT的海战多目标态势感知. 信息与控制. 2023(05): 638-647 . 百度学术4. 吴智怀，王志伟. 基于SDAE的液压泵故障诊断方法研究. 冶金设备. 2023(05): 15-22 . 百度学术5. 李锦键，王兴贵，杨维满，赵玲霞. 基于改进递归深度信念网络的CSP电站短期出力预测. 太阳能学报. 2022(07): 225-232 . 百度学术6. 沈浩宇，江先志，冯涛. 二轮平衡车的神经元控制算法及仿真. 机床与液压. 2022(19): 159-166 . 百度学术7. 王功明，乔俊飞，关丽娜，贾庆山. 深度信念网络研究现状与展望. 自动化学报. 2021(01): 35-49 . 本站查看8. 石媛媛，裴志利，姜明洋. 基于深度信念网络的文本分类研究综述. 内蒙古民族大学学报(自然科学版). 2021(01): 44-50 . 百度学术9. 李娜娜，胡坚剑，顾军华，张亚娟. 深度置信网络优化模型在人才评价中的应用. 计算机工程. 2020(02): 80-87+102 . 百度学术10. 唐魏，郑源，潘虹，徐晶珺. 引入动态调节学习率的SAE轴承故障诊断研究. 计算机工程与应用. 2020(20): 264-269 . 百度学术11. 李少波，李传江，胡建军，张安思，杨静. 基于深度置信网络的机械设备故障诊断研究综述. 现代制造工程. 2020(10): 156-162+142 . 百度学术12. 史科，陆阳，刘广亮，毕翔，王辉. 基于多隐层Gibbs采样的深度信念网络训练方法. 自动化学报. 2019(05): 975-984 . 本站查看13. 高磊，范冰冰，黄穗. 基于残差的改进卷积神经网络图像分类算法. 计算机系统应用. 2019(07): 139-144 . 百度学术14. 张德，李国璋，王怀光，张峻宁. 位姿估计自适应学习率的改进. 电子测量与仪器学报. 2019(06): 51-58 . 百度学术15. 邢海霞，程乐. 一种基于强化学习的深度信念网络设计方法. 控制工程. 2019(11): 2115-2120 . 百度学术16. 索明何，程乐. 基于PLSR的深度信念网输出权值确定方法. 控制工程. 2018(04): 668-676 . 百度学术17. 王功明，乔俊飞，王磊. 一种能量函数意义下的生成式对抗网络. 自动化学报. 2018(05): 793-803 . 本站查看18. 牟善仲，徐天赐，符奥，王萌，白茹. 基于自适应深度学习模型的变压器故障诊断方法. 南方电网技术. 2018(10): 14-19 . 百度学术19. 刘冰尧，庞腾，雷菊阳. 基于深度信念网络的变压器故障诊断. 化工自动化及仪表. 2018(11): 873-878 . 百度学术20. 王献锋，张传雷，张善文，朱义海. 基于自适应判别深度置信网络的棉花病虫害预测. 农业工程学报. 2018(14): 157-164 . 百度学术其他类型引用(31)
-
下载:
计量
- 文章访问数: 2017
- HTML全文浏览量: 246
- PDF下载量: 876
- 被引次数: 51
