Matlab-Matpower制作IEEE14-电力虚假数据注入攻击FDIA数据集

文章目录

- - 1. 加载Matpower-IEEE14电力数据
  - 2. 导入原始数据集
  - - 两个数据集结合的意义
    - 潮流分析和状态估计的意义
  - 3. 初始化变量
  - 4. 分离有功和无功功率
  - 4. 潮流计算
  - 5. 生成测量向量
  - 6. 选择是否篡改数据
  - 7. 状态估计和雅可比矩阵
  - 8. 保存未篡改数据
  - 9. 篡改数据
  - 生成FDIA仿真数据集完整代码
  - IEEE14节点数据
  - 状态估计函数

MatpowerIEEE14_2">1. 加载Matpower-IEEE14电力数据

matlab">mpc = loadcase('case14');

通过matpower的loadcase函数加载 case14 文件（电力系统IEEE14节点数据案例），
mpc 是一个包含电力系统参数的结构体。该案例定义了一个14节点的电网结构，包括节点类型、发电机信息、线路信息、有功功率、无功功率等测试数据。这个案例可作为FDIA、电力潮流分析等研究的基础。

2. 导入原始数据集

matlab">rawDataset = readmatrix('RegroupedDataset\RegroupedData.csv');

读取存储在 CSV 文件中的原始数据集 RegroupedData.csv，该数据集源自该档案包含 2006 年 12 月至 2010 年 11 月（47 个月）期间在位于 Sceaux（法国巴黎 7 公里）的一所房屋中收集的 2075259 个测量值。

1.（global_active_power*1000/60 - sub_metering_1 - sub_metering_2 - sub_metering_3）表示家庭中未在子计量 1、2 和 3 中测量的电气设备每分钟消耗的有功电能（以瓦时为单位）。
2.数据集包含一些测量值缺失值（接近 1.25% 的行）。所有日历时间戳都存在于数据集中，但对于某些时间戳，测量值缺失：缺失值表示两个连续的分号属性分隔符之间没有值。例如，数据集显示 2007 年 4 月 28 日有缺失值。

数据集下载链接：个人家庭用电量公开数据集

主要包含不同时间步的活动功率和无功功率数据，将用于仿真电力系统中的负荷变化。

两个数据集结合的意义

外部数据集：提供时间序列中的负荷变化，确保仿真过程中的负荷是动态的。
IEEE 14节点案例：提供标准电力网络结构，使得仿真能够在一个现实的电网拓扑上。

外部数据集 (RegroupedData.csv)
个人家庭用电量CSV文件中的数据包含不同时间步（时刻）的负荷数据，即有功功率（active power）和无功功率（reactive power）的测量值。用于电力系统仿真分析，电力系统中的负荷是动态变化的，因为实际电力系统中的用电负荷并不是恒定的，而是随着时间波动的，和交通流量的早晚高峰类似。
IEEE 14节点案例 (case14)
IEEE 14节点案例是一个标准的电力系统测试模型，常用于电力系统仿真与分析。它定义了电网的拓扑结构，包括节点（bus）、发电机位置、输电线路等。这些参数相对静态，描述了电力系统在物理上如何连接、如何传输电力。

系统结构：case14 文件定义了14个节点的电力网络，包括发电机的位置（通常是节点7和8），负荷节点，线路阻抗等。这些信息在仿真过程中保持不变。
使用方式：IEEE 14节点系统中定义了电网的基础结构，根据外部数据集更新某些节点的负荷，并执行潮流计算来模拟不同时间步电力系统的状态。

动态负荷+静态电网结构：外部CSV文件提供了时间步变化的负荷数据，而EEE 14节点提供了一个静态的电力网络结构。
外部数据用于在每个时间步更新IEEE 14节点系统中的负荷，从而进行潮流计算。使得每个时间步的潮流分析基于在同一个电力网络结构（14节点系统）上，进行不同的负荷配置。

潮流分析和状态估计的意义

FDIA攻击的核心目标是通过篡改电力系统的测量数据，干扰状态估计的结果，从而使得系统误判当前的状态。为达到这一目的，需要通过对IEEE节点系统和负荷变化的数据集进行不同时间步的潮流计算和状态估计。

潮流分析和状态估计：每次读取外部数据集中的一组负荷数据后，都会根据这些负荷在IEEE 14节点系统上运行潮流分析，并生成测量向量用于状态估计。这些测量数据根据一定的选择概率被伪造（通过添加攻击向量），生成不同的FDIA测试数据集。

3. 初始化变量

matlab">untamperedVectors = [];
tamperedVectors = [];
labels = [];
actives = [];
reactives = [];
data_len = size(rawDataset,1);

初始化存储未篡改、篡改的测量向量及其对应标签的数组。
actives 和 reactives 用于存储每个时间步的有功和无功功率测量值。
data_len 是数据集中时间步的总数。

4. 分离有功和无功功率

matlab">for j = 1:data_len % j is the Time step%--------------------------------------------%% Split active and reactive power measurementstempActives = [];tempReactives = [];for i = 1:22if mod(i,2) ~= 0tempActives = [tempActives, rawDataset(j,i)];elsetempReactives = [tempReactives, rawDataset(j,i)];endendactives = [actives; tempActives];reactives = [reactives; tempReactives];

遍历数据集中每个时间步的22个数据点，将奇数索引的数据视为有功功率，偶数索引的数据视为无功功率，分别存储到 tempActives 和 tempReactives 中。

4. 潮流计算

matlab">savecase('CustomCase14.m', mpc);
results = rundcpf(mpc);

保存更新后的电力系统模型 mpc，然后通过 rundcpf 运行潮流分析，计算出每个节点的功率注入和每条线路的功率流。

5. 生成测量向量

matlab">realPowerInjections = results.bus(:,3);
realPowerFlows = results.branch(:,14);
measurementVector = [realPowerInjections; realPowerFlows];

从潮流分析结果中提取每个节点的有功功率注入（realPowerInjections）和每条线路的有功功率流（realPowerFlows），将它们组合成测量向量。

6. 选择是否篡改数据

matlab">r = randi([0 100], 1, 1);
if(r < 20) && (j > floor(data_len*0.7))decision = true;
elsedecision = false;
end
labels = [labels; decision];

使用随机数生成器决定当前时间步是否篡改测量向量。每个时间步有20%的几率被篡改，且仅在数据集的最后30%时间步上进行篡改。decision 变量记录是否篡改，并将其保存到 labels 标签数组中。

7. 状态估计和雅可比矩阵

matlab">ser = StateEstimation(results);
H = makeJac(mpc);

使用 StateEstimation 函数进行状态估计，通过 makeJac 计算雅可比矩阵 H，它表示测量与状态变量之间的线性关系。

8. 保存未篡改数据

matlab">untamperedVectors = [untamperedVectors, measurementVector];

9. 篡改数据

如果决定篡改，则通过以下步骤生成伪造的测量向量：

1. 计算误差向量
1. 生成攻击向量
1. 生成篡改的测量向量

matlab">za = measurementVector + (H*c) + error;
status = strcat(int2str(j), '/', int2str(data_len), 'simulations complete');  % 进度显示
disp(status)
writematrix(untamperedVectors, 'VectorDataset2\untamperedVectorData.csv')
writematrix(tamperedVectors, 'VectorDataset2\tamperedVectorData.csv')
writematrix(labels, 'VectorDataset2\labelData.csv')

最后，保存未篡改和篡改后的测量向量，以及标签数据到 CSV 文件中。

FDIA_138">生成FDIA仿真数据集完整代码

matlab">function Simulation%------------------------------------------------------------------------------%
% One-off Declarations
mpc = loadcase('case14');%------------------------------------------------------------------------------%
% Import raw dataset
rawDataset = readmatrix('RegroupedDataset\RegroupedData.csv');%------------------------------------------------------------------------------%
% Power Flow Analysis
untamperedVectors = [];
tamperedVectors = [];
labels = [];actives = [];
reactives = [];
data_len = size(rawDataset,1);for j = 1:data_len % j is the Time step%--------------------------------------------%% Split active and reactive power measurementstempActives = [];tempReactives = [];for i = 1:22if mod(i,2) ~= 0tempActives = [tempActives, rawDataset(j,i)];elsetempReactives = [tempReactives, rawDataset(j,i)];endendactives = [actives; tempActives];reactives = [reactives; tempReactives];%--------------------------------------------%% Set the case file parameterscell = 1;index = j;for i = 2:14if (i ~= 7) && (i ~= 8)mpc.bus(i,2) = 1;  % Set as load busmpc.bus(i,3) = actives(index,cell);mpc.bus(i,4) = reactives(index,cell);cell = cell + 1;endend% Set as generator busesmpc.bus(7,2) = 2;mpc.bus(8,2) = 2;% Save the casefile and run the power flow analysissavecase('CustomCase14.m', mpc);results = rundcpf(mpc);%--------------------------------------------%% Create the measurement vectorrealPowerInjections = results.bus(:,3);realPowerFlows = results.branch(:,14);measurementVector = [realPowerInjections; realPowerFlows];%--------------------------------------------%% Decide if measurement vector is to be falsifiedr = randi([0 100], 1, 1); % Generate a random number from 0 to 100if(r < 20) && (j > floor(data_len*0.7)) % Do not falsify the first 2k measurements of the testing datasetdecision = true; % 20 percent chance of falsificationelsedecision = false;endlabels = [labels; decision];%--------------------------------------------%% Perform state estimationser = StateEstimation(results);%--------------------------------------------%% Compute the Jacobian matrix for the bus systemH = makeJac(mpc);jacLen = length(H);%--------------------------------------------%% Save measurement vector to untamperedVectorsmeasurementVector = measurementVector(1:jacLen);untamperedVectors = [untamperedVectors, measurementVector];%--------------------------------------------%% Falsification if decidedif j > floor(data_len*0.7) % Falsify only the testing portion (last 10%)if ~decision% Append the untampered measurement vectortamperedVectors = [tamperedVectors, measurementVector];else% Compute error vectorerror = [];for i = 1:14error = [error; 0];endser.z = ser.z(1:20);ser.z_est = ser.z_est(1:20);flowsError = ser.z - ser.z_est;error = [error; flowsError];error = error(1:jacLen);% Compute attack vectorr = randi([1 10], 1, 1); % Generate a random number from 0 to 100c = [];for i = 1:jacLenc = [c, r];endc = transpose(c);% Compute the falsified measurement vectorza = measurementVector + (H*c) + error;% Append the tampered measurement vectortamperedVectors = [tamperedVectors, za];endendstatus = strcat(int2str(j), '/', int2str(data_len), 'simulations complete');disp(status)
enduntamperedVectors = transpose(untamperedVectors);
tamperedVectors = transpose(tamperedVectors);%------------------------------------------------------------------------------%
% Save results
writematrix(untamperedVectors, 'VectorDataset2\untamperedVectorData.csv')
writematrix(tamperedVectors, 'VectorDataset2\tamperedVectorData.csv')
writematrix(labels, 'VectorDataset2\labelData.csv')

代码中的 if(r < 20) && (j > floor(data_len*0.7)) 条件决定了在一部分时间步（后30%）中，有20%的几率进行数据篡改，而前70%的数据保持真实。

这样做有两个目的：

用于训练和测试：在使用时间序列模型如LSTM、GRU做攻击检测的过程中，通常使用正常数据集进行模型训练，而测试时使用含有异常数据的数据集进行差异性对比从而实现攻击检测。

即异常检测（Anomaly Detection），通过学习正常数据的统计特性（如均值、标准差等），当某些数据偏离正常范围时，标记为异常。

避免篡改数据过早被发现：如果一开始就篡改数据，可能会导致防御系统过早发现问题，降低攻击的隐蔽性。

IEEE14节点数据

matlab">function mpc = CustomCase14
%CUSTOMCASE14%% MATPOWER Case Format : Version 2
mpc.version = '2';%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = 100;%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [1	3	0	0	0	0	1	1.06	0	0	1	1.06	0.94;2	1	0.346	0.298	0	0	1	1.045	-4.98	0	1	1.06	0.94;3	1	0.3	0.252	0	0	1	1.01	-12.72	0	1	1.06	0.94;4	1	0.716	0.202	0	0	1	1.019	-10.33	0	1	1.06	0.94;5	1	0.606	0.14	0	0	1	1.02	-8.78	0	1	1.06	0.94;6	1	1.298	0.116	0	0	1	1.07	-14.22	0	1	1.06	0.94;7	2	0	0	0	0	1	1.062	-13.37	0	1	1.06	0.94;8	2	0	0	0	0	1	1.09	-13.36	0	1	1.06	0.94;9	1	3.692	0	0	19	1	1.056	-14.94	0	1	1.06	0.94;10	1	0.412	0.112	0	0	1	1.051	-15.1	0	1	1.06	0.94;11	1	1.298	0.046	0	0	1	1.057	-14.79	0	1	1.06	0.94;12	1	0.67	0	0	0	1	1.055	-15.07	0	1	1.06	0.94;13	1	0.49	0.142	0	0	1	1.05	-15.16	0	1	1.06	0.94;14	1	0.938	0	0	0	1	1.036	-16.04	0	1	1.06	0.94;
];%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
mpc.gen = [1	232.4	-16.9	10	0	1.06	100	1	332.4	0	0	0	0	0	0	0	0	0	0	0	0;2	40	42.4	50	-40	1.045	100	1	140	0	0	0	0	0	0	0	0	0	0	0	0;3	0	23.4	40	0	1.01	100	1	100	0	0	0	0	0	0	0	0	0	0	0	0;6	0	12.2	24	-6	1.07	100	1	100	0	0	0	0	0	0	0	0	0	0	0	0;8	0	17.4	24	-6	1.09	100	1	100	0	0	0	0	0	0	0	0	0	0	0	0;
];%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [1	2	0.01938	0.05917	0.0528	0	0	0	0	0	1	-360	360;1	5	0.05403	0.22304	0.0492	0	0	0	0	0	1	-360	360;2	3	0.04699	0.19797	0.0438	0	0	0	0	0	1	-360	360;2	4	0.05811	0.17632	0.034	0	0	0	0	0	1	-360	360;2	5	0.05695	0.17388	0.0346	0	0	0	0	0	1	-360	360;3	4	0.06701	0.17103	0.0128	0	0	0	0	0	1	-360	360;4	5	0.01335	0.04211	0	0	0	0	0	0	1	-360	360;4	7	0	0.20912	0	0	0	0	0.978	0	1	-360	360;4	9	0	0.55618	0	0	0	0	0.969	0	1	-360	360;5	6	0	0.25202	0	0	0	0	0.932	0	1	-360	360;6	11	0.09498	0.1989	0	0	0	0	0	0	1	-360	360;6	12	0.12291	0.25581	0	0	0	0	0	0	1	-360	360;6	13	0.06615	0.13027	0	0	0	0	0	0	1	-360	360;7	8	0	0.17615	0	0	0	0	0	0	1	-360	360;7	9	0	0.11001	0	0	0	0	0	0	1	-360	360;9	10	0.03181	0.0845	0	0	0	0	0	0	1	-360	360;9	14	0.12711	0.27038	0	0	0	0	0	0	1	-360	360;10	11	0.08205	0.19207	0	0	0	0	0	0	1	-360	360;12	13	0.22092	0.19988	0	0	0	0	0	0	1	-360	360;13	14	0.17093	0.34802	0	0	0	0	0	0	1	-360	360;
];%%-----  OPF Data  -----%%
%% generator cost data
%	1	startup	shutdown	n	x1	y1	...	xn	yn
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [2	0	0	3	0.0430292599	20	0;2	0	0	3	0.25	20	0;2	0	0	3	0.01	40	0;2	0	0	3	0.01	40	0;2	0	0	3	0.01	40	0;
];%% bus names
mpc.bus_name = {'Bus 1     HV';'Bus 2     HV';'Bus 3     HV';'Bus 4     HV';'Bus 5     HV';'Bus 6     LV';'Bus 7     ZV';'Bus 8     TV';'Bus 9     LV';'Bus 10    LV';'Bus 11    LV';'Bus 12    LV';'Bus 13    LV';'Bus 14    LV';
};

基本系统参数
mpc.version: 系统的MATPOWER格式版本，通常是 ‘2’。
mpc.baseMVA: 系统的基准功率，单位是MVA。在这个例子中，基准功率是100 MVA。
Bus Data (节点数据)
每个节点的参数描述了该节点的电力负荷、电压、区域等信息。它们的列意义如下：

bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin

1、bus_i: 节点编号
2、type: 节点类型

1：负荷节点（PQ 节点）
2：发电机节点（PV 节点）
3：平衡节点（Slack 节点）

3、Pd: 有功负荷，单位为 MW（兆瓦）
4、Qd: 无功负荷，单位为 MVar（兆乏）
5、Gs: 并联到节点的电导，单位为 S（西门子），通常为零
6、Bs: 并联到节点的电纳，单位为 S（西门子），通常为零
7、area: 区域编号，通常为 1，表示系统中节点所在的区域
8、Vm: 节点电压幅值，单位为标幺（pu），无量纲
9、Va: 节点电压相角，单位为度（degree）
10、baseKV: 节点的基准电压，单位为千伏（kV）
11、zone: 区域编号，通常为 1
12、Vmax: 节点允许的最大电压，单位为标幺
13、Vmin: 节点允许的最小电压，单位为标幺

例如，第一个节点的描述如下：

1	3	0	0	0	0	1	1.06	0	0	1	1.06	0.94

这个节点是平衡节点（Slack Node），有功负荷和无功负荷都为 0，电压为 1.06 pu，相角为 0 度。

Generator Data (发电机数据)
发电机的数据描述了每个发电机的功率输出、电压控制和相关约束条件。每列的含义如下：

bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf

1、bus: 发电机所在的节点编号
2、Pg: 发电机的有功功率输出，单位为 MW（兆瓦）
3、Qg: 发电机的无功功率输出，单位为 MVar（兆乏）
4、Qmax: 无功功率的上限，单位为 MVar
5、Qmin: 无功功率的下限，单位为 MVar
6、Vg: 发电机控制的电压设定值，单位为标幺
7、mBase: 发电机的基准容量，通常为 100 MVA
8、status: 发电机的状态，1 表示在线，0 表示离线
9、Pmax: 发电机的有功功率上限，单位为 MW
10、Pmin: 发电机的有功功率下限，单位为 MW
11、Pc1, Pc2: 启停功率控制参数，通常为 0
12、Qc1min, Qc1max, Qc2min, Qc2max: 启停无功功率控制参数，通常为 0
13、ramp_agc, ramp_10, ramp_30, ramp_q: 发电机的爬坡率，单位分别为 MW/min 或 MVar/min
14、apf: 发电机的参与因子，通常为 0

例如：发电机 1 的数据为：

1	232.4	-16.9	10	0	1.06	100	1	332.4	0	0	0	0	0	0	0	0	0	0	0	0

这表示发电机 1 在节点 1 上输出 232.4 MW 的有功功率和 -16.9 MVar 的无功功率，最大无功功率限制为 10 MVar，电压设定为 1.06 pu.

Branch Data (支路数据)
支路数据描述了系统中各条线路（包括传输线和变压器）的阻抗、传输能力等信息。每列的含义如下：

fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax

1、fbus: 起始节点编号
2、tbus: 终止节点编号
3、r: 支路电阻，单位为标幺（pu）
4、x: 支路电抗，单位为标幺（pu）
5、b: 支路的电纳，单位为标幺（pu）
6、rateA, rateB, rateC: 支路的热容量，单位为 MVA（不同的额定值通常表示不同的运行条件）
7、ratio: 变压器的变比，通常为 1
8、angle: 变压器的相位角，通常为 0
9、status: 支路的状态，1 表示在线，0 表示离线
10、angmin, angmax: 支路允许的最小和最大相角差，单位为度

例如，支路 1-2 的数据为：

1	2	0.01938	0.05917	0.0528	0	0	0	0	0	1	-360	360

这表示从节点 1 到节点 2 的支路电阻为 0.01938 pu，电抗为 0.05917 pu，电纳为 0.0528 pu.

Generator Cost Data (发电机成本数据)
发电机成本数据用于描述发电机的发电成本函数。常见的形式是二次成本函数，其表达式为：

$C(P_g) = a P_g^2 + b P_g + c$

1	startup	shutdown	n	x1	y1	...	xn	yn
2	startup	shutdown	n	c(n-1)	...	c0

startup, shutdown: 启动和关闭成本，通常为 0
n: 成本函数的阶数
c(n-1) … c0: 成本函数的系数

例如，第一个发电机的成本数据为：

2	0	0	3	0.0430292599	20	0

这表示其成本函数为 $0.0430292599 \cdot P_g^2 + 20 \cdot P_g+0$

Bus Names (节点名称)
节点名称为每个节点提供了文字描述，通常是为了方便识别。每个节点都对应一个名称，比如：

'Bus 1     HV';

表示节点 1 的名称是 “Bus 1 HV”。

状态估计函数

matlab">function ser = StateEstimation(results)%% which measurements are available
idx.idx_zPF = [1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20];
idx.idx_zPT = [4;5;7;11];
idx.idx_zPG = [1;2;3;4;5];
idx.idx_zVa = [];
idx.idx_zQF = [1;3;8;9;10;13;15;19];
idx.idx_zQT = [4;5;7;11];
idx.idx_zQG = [1;2];
idx.idx_zVm = [2;3;6;8;10;14];%% specify measurements
measure.PF = results.branch(:,14);measure.PT = [];
for i = 1: length(idx.idx_zPT)measure.PT = [measure.PT; results.branch(idx.idx_zPT(i),16);];
endmeasure.PG = [];
for i = 1: length(idx.idx_zPG)measure.PG = [measure.PG; results.gen(idx.idx_zPG(i),2);];
endmeasure.Va = [];measure.QF = [];
for i = 1: length(idx.idx_zQF)measure.QF = [measure.QF; results.branch(idx.idx_zQF(i),15);];
endmeasure.QT = [];
for i = 1: length(idx.idx_zQT)measure.QT = [measure.QT; results.branch(idx.idx_zQT(i),17);];
endmeasure.QG = [];
for i = 1: length(idx.idx_zQG)measure.QG = [measure.QG; results.branch(idx.idx_zQG(i),3);];
endmeasure.Vm = [1;1;1;1;1;1];%% specify measurement variances
sigma.sigma_PF = 0.02;
sigma.sigma_PT = 0.02;
sigma.sigma_PG = 0.015;
sigma.sigma_Va = [];
sigma.sigma_QF = 0.02;
sigma.sigma_QT = 0.02;
sigma.sigma_QG = 0.015;
sigma.sigma_Vm = 0.01;%% check input data integrity
nbus = 14;
[success, measure, idx, sigma] = checkDataIntegrity(measure, idx, sigma, nbus);
if ~successerror('State Estimation input data are not complete or sufficient!');
end%% run state estimation
casename = 'CustomCase14.m';
type_initialguess = 2; % flat start
[baseMVA, bus, gen, branch, success, et, z, z_est, error_sqrsum] = run_se(casename, measure, idx, sigma, type_initialguess);
ser.z = z;
ser.z_est = z_est;