前言
最近的课题需要做的是信息物理系统的攻击检测类和安全估计类算法。以经典的Kalman算法为例。动态系统模型如下:
X
k
+
1
=
A
X
k
+
B
U
k
+
ω
Y
k
=
C
X
k
+
v
\begin{matrix}X_{k+1}=AX_{k}+BU_{k}+\omega \\ Y_{k}=CX_{k}+v \end{matrix}
Xk+1=AXk+BUk+ωYk=CXk+v
A,B,C矩阵分别为状态转移矩阵、输入矩阵和量测矩阵。假设更简单的情况,没有系统输入的情况下,系统的模型如下:
X k + 1 = A X k + ω Y k = C X k + v \begin{matrix}X_{k+1}=AX_{k}+\omega \\ Y_{k}=CX_{k}+v \end{matrix} Xk+1=AXk+ωYk=CXk+v
前一段时间纠结了很久的利用标准IEEE电网模型去做这类算法的仿真,实际上在比较正式的IEEE期刊论文当中,进行仿真的也都是去利用IEEE的标准模型来做的。
一开始在网上收集到的是利用Matlab的Simulink工具箱对电网模型的展开。但是相对比较复杂。在咨询以及翻阅包括自己收集到的一写论文之后发现基本都是用基于Matlab的一个叫Matpower的工具包。这个工具包的用处和原理在前面的文章中给出。
Matpower软件简介和参数介绍
MATPOWER工具本质原理解析
MATPOWER用于动态安全估计的仿真
经过一段时间的查阅资料和思考猜测,思路大概有一下两种:
模型的动态由实际负载等参数决定
首先需要说明的是,我们可以确定Matpower进行潮流计算后产生的数据是作为量测值来使用的,而且我们设计的攻击也是叠加在此之上。
我们做的是动态估计,假设的系统模型中,系统的状态是随时间变化的(即使是准稳态下)。这种动态的方式是通过负载等参数的随时间变化来设计,经过Matpower的潮流计算后,会体现在节点参数,也就是量测值当中。
用系统函数来体现,即状态量X表征的是负载等参数,量测值Y表征的就是潮流计算之后的节点参数。状态转移矩阵A是预设的负载的动态变化,这种变化可以有两种方法来实现:
- 自己假定A的值,比如作为准静态的过程可以假设A为单位矩阵,所有的动态都是叠加上去的高斯白噪声;
- 直接调取某些现有的实际电网中的历史数据,这种情况下使用状态转移方程的原形式是无法获得状态转移矩阵A的,因此可以使用电力行业中经常使用的Holt‘s方法,利用历史数据来对状态的动态进行拟合。
通过以上可以解决状态转移矩阵A的动态问题。在这种假设之下,量测矩阵表征的实际上就是潮流运算,量测矩阵H隐含在Matpower的潮流运算过程中。因此即使解决了状态转移方程中状态转移矩阵A的问题,如何在Matpower的潮流运算中提取出量测矩阵C依然是严峻的问题。
就目前搜集到的资料,猜测:
- 与Matpower运算过程中的中间雅可比矩阵有关;
- 根据电网模型有特定的现有量测矩阵,且在行业内是基础知识。
模型的动态自己来设置
这种情况下,实际上是将Matpower进行潮流运算后获取到的节点参数既作为状态值也作为量测值来使用,因此这种方式下,系统的状态转移矩阵A由自己来设置,同时系统的量测矩阵C也常常设置为单位矩阵。
最简单的方式是,假设电网就处于准静态情况下,出现的动态都是因为故障或攻击。因此A矩阵为单位矩阵,状态值包含噪声;C矩阵同样为单位矩阵,实际包含攻击的系统模型如下:
X
k
+
1
=
X
k
+
ω
Y
k
=
X
k
+
v
+
a
\begin{matrix}X_{k+1}=X_{k}+\omega \\ Y_{k}=X_{k}+v+a \end{matrix}
Xk+1=Xk+ωYk=Xk+v+a
其中,a为精心设计的攻击矩阵。
具体仿真的两种形式
Matpower
上述
自建模型
也有的学位论文并没有直接采用IEEE的标准电网的Matpower工具,而是根据IEEE标准电网的模型,进行了非常专业的电力系统建模,然后通过建模实现了量测矩阵C的获取。但是这种的前提是:
- 电网不能太复杂,最多是IEEE9节点电网模型;
- 系统的动态过程依然需要Holt’s方法来获取。(电力系统本身就是一个庞大的非线性系统,获取A很困难)
最后
说到这里,实际上如果研究的是线性系统,基本上可以利用IEEE的标准电网和Matpower工具包做的,只能是自己设置动态,也就是
X
k
+
1
=
X
k
+
ω
Y
k
=
X
k
+
v
+
a
\begin{matrix}X_{k+1}=X_{k}+\omega \\ Y_{k}=X_{k}+v+a \end{matrix}
Xk+1=Xk+ωYk=Xk+v+a
这种形式。而且电网模型实际上对于整个CPS的研究来说并不是完全具备典型性。