5.2 电力系统状态估计

1、案例概况

本案例介绍利用mems实现matpower中的状态估计计算，包括从scada数据源获取量测数据、网络设备信息查询、状态估计计算、计算结果写入状态估计se数据源。

matpower状态估计功能包含以下主要函数：

• run_se: 运行状态估计
• doSE: 求解状态估计问题
• checkDataIntegrity：校验输入数据完整性
• isobservable：校验系统是否可观测

上述函数均已通过RustScript语言实现，并在github中上传，可在mems中直接调用，URL如下：

• run_se: https://shufengdong.github.io/sparrowzz/rspower/lib/runse.txt
• doSE: https://shufengdong.github.io/sparrowzz/rspower/lib/doSE.txt
• checkDataIntegrity：https://shufengdong.github.io/sparrowzz/rspower/lib/checkDataIntegrity.txt
• isobservable：https://shufengdong.github.io/sparrowzz/rspower/lib/isobservable.txt

状态估计计算输入数据中，除潮流计算的输入数据以外，还包括以下8类量测数据：

• PF：支路首端有功功率量测（由首端流向尾端）
• QF：支路首端无功功率量测（由首端流向尾端）
• PT：支路尾端有功功率量测（由首端流向尾端）
• QT：支路尾端无功功率量测（由首端流向尾端）
• Vm：母线电压幅值量测
• Va：母线电压相角量测
• PG：发电机有功出力量测
• QG：发电机无功出力量测

这8类量测可分别设置其测量误差。

在输入数据中通过位置索引变量确定量测与设备的对应关系，例如idx.idx_zPF = [1;2]、measure.PF = [0.12;0.10]中，idx_zPF指定了PF向量中的2个量测分别为branch矩阵中第1条和第2条支路的量测。

2、量测配置

状态估计案例中，bus、branch、gen设备信息的导入与潮流计算类似，除此之外，还需要配置量测信息。由第1节可知，每条支路包含PF、QF、PG、QG四个量测，每个母线包含Vm、Va两个量测，每台发电机包含PG、QG两个量测。对于每个量测，在mems界面设置项-测点中导入对应的测点。其中，测点单位需正确配置，后续需通过单位区分同一元件的不同量测物理量。对于存在量测值的测点，设置其默认值为量测值，对于不存在量测值的测点，设置其默认值为NaN。

此外，系统基准容量和8类量测的误差参数也作为测点数据进行配置，后续计算从测点中获取。

然后，在设备-设备维护-测点导入设备的量测信息，建立设备与测点的关联关系。量测信息配置文件如下图所示：

其中，ID列为量测id，Point ID列为对应测点点号，Dev ID为对应设备ID，Terminal ID为对应端子ID，Phase为量测相别。注意导入量测后需要提交并应用版本。

3、状态估计计算配置

在设置项-报表配置实现状态估计计算的DataFrame流，包括获取量测数据、获取网络设备信息、执行状态估计计算、输出结果等，以下说明具体配置步骤。

（1）获取测点id和单位信息

报表配置文件中，在节点1中获取测点id和单位信息，如下图所示。

节点1类型为SOURCE_POINT，表示获取测点表信息，配置参数为如下sql语句：

select id,data_unit from points order by id

其中，points为测点表，id为测点id列名，data_unit为测点单位列名。

（2）获取scada数据源实时量测信息

在节点2中获取scada数据源实时量测信息，如下图所示。

节点2类型为SOURCE_MEAS，表示获取量测表信息。参数1为0，表示获取数据源为0（即scada数据源）的量测数据，参数2为如下sql语句：

select id,value from meas order by id id

其中，meas为量测表，id为量测对应测点id列名，value为量测值列名。

（3）获取母线量测信息

通过节点3、4、5、6和边1;4、2;4、3;4、4;5、4;6获取母线电压幅值和电压相角量测，如下图所示。

节点3类型为SOURCE_DEV，表示获取设备表信息，参数为如下sql语句：

select t2.dev,t2.terminal,t3.point from Bus t1 
left join topo t2 on t1.id=t2.dev 
left join measdef t3 on t2.terminal = t3.terminal 
order by t1.id

其中，Bus为母线表；topo为拓扑信息表，包含设备id及其对应端子id信息；measdef为量测定义表，包含测点点号及其对应端子id对应关系。

节点4类型为TensorEval，参数2为4，表示进行DataFrame操作。在该节点中对查询到的测点信息、量测值、设备信息进行组合操作，参数1配置如下：

bus_points_unit = join(bus_points, points, `point`, `id`, inner);
bus_meas_all = join(bus_points_unit, meas, `point`, `id`, inner);

在该节点中先对设备表信息与测点表信息按设备表的顺序进行拼接，再与量测表信息进行拼接，均按照测点id字段进行拼接。

边4;5和4;6按照量测单位对母线2个物理量进行区分，分别将电压幅值量测和电压相角量测信息传给节点5和节点6。以边4;5为例，其动作类型为Eval，动作参数如下：

filter(col(dev), col(terminal), col(point), col(value), col(data_unit)==`UnitOne`);

其中，filter函数用于对DataFrame进行过滤筛选，最后一个参数为过滤条件表达式，UnitOne为电压幅值标幺值单位，之前的参数为保留的列。

（4）获取支路首端量测信息

通过节点7、8、9、10和边1;8、2;8、7;8、8;9、8;10获取支路首端有功和无功潮流量测，如下图所示。

获取支路测点信息整体过程与（3）中母线操作相同，不同之处是还需要区分支路首端和尾端量测，为此在节点7配置参数如下：

SELECT 
    t2.dev,
    t2.terminal,
    t3.point
FROM Branch t1 
LEFT JOIN (
    SELECT 
        dev,
        terminal
    FROM (
        SELECT 
            dev,
            terminal,
            ROW_NUMBER() OVER (PARTITION BY dev ORDER BY terminal) AS rn
        FROM topo
    ) tmp
    WHERE rn = 1
) t2 ON t1.id = t2.dev
LEFT JOIN measdef t3 ON t2.terminal = t3.terminal
ORDER BY t1.id

其中，Branch为支路表。该案例中根据端子id大小区分支路首端和尾端，端子id小的表示首端，端子id大的表示尾端。通过子查询t2获取每个支路较小的端子id，从而获取支路首端测点信息，具体逻辑为：

• 从topo表中，按dev即设备id分组
• 每组内按terminal排序，用ROW_NUMBER()生成编号
• 取每组第1行（rn=1），即每个设备对应的最小端子id值
• 然后与Branch表通过t1.id=t2.dev左连接。

（5）获取支路尾端量测信息

通过节点11、12、13、14和边1;12、2;12、11;12、12;13、12;14获取支路首端有功和无功潮流量测，过程与获取支路首端量测信息类似，如下图所示。

（6）获取发电机量测信息

通过节点15、16、17、18和边1;16、2;16、15;16、16;17、16;18获取发电机有功和无功出力量测，如下图所示。

（7）获取量测误差信息

从测点表中获取各类量测误差信息，如下图所示：

（8）获取设备参数和系统基准容量

从设备表中获取各类设备参数，对应matpower的bus、branch、gen、baseMVA输入参数，如下图所示：

（9）输入数据传给状态估计计算节点

通过边将上述数据传给状态估计计算节点，配置如下图所示：

其中对于量测需要滤除NaN值，以边5;27过滤节点电压幅值中的NaN值为例，边类型为Eval，配置参数如下

drop_nans();

该函数用于删除DataFrame中包含NaN值的行。

（10）进行状态估计计算

在节点27中进行状态估计计算，该节点类型为TensorEval，参数2为3，表示进行复数张量运算，参数1配置及说明如下：

idx_zVm = select(bus_meas_Vm, [], 0)-10000; //母线电压幅值量测对应设备位置索引
measure_Vm = select(bus_meas_Vm, [], 3);  //母线电压幅值量测值
idx_zVa = select(bus_meas_Va, [], 0)-10000; //母线电压相角量测对应设备位置索引
measure_Va = select(bus_meas_Va, [], 3);  //母线电压相角量测值
idx_zPF = select(brf_meas_Pf, [], 0)-30000; //支路首端有功量测对应设备位置索引
measure_PF = select(brf_meas_Pf, [], 3);  //支路首端有功量测值
idx_zQF = select(brf_meas_Qf, [], 0)-30000; //支路首端无功量测对应设备位置索引
measure_QF = select(brf_meas_Qf, [], 3);  //支路首端无功量测值
idx_zPT = select(brt_meas_Pt, [], 0)-30000; //支路尾端有功量测对应设备位置索引
measure_PT = select(brt_meas_Pt, [], 3);  //支路尾端有功量测值
idx_zQT = select(brt_meas_Qt, [], 0)-30000; //支路尾端无功量测对应设备位置索引
measure_QT = select(brt_meas_Qt, [], 3);  //支路尾端无功量测值
idx_zPG = select(gen_meas_Pg, [], 0)-20000; //发电机有功出力量测对应设备位置索引
measure_PG = select(gen_meas_Pg, [], 3);  //发电机有功出力量测值
idx_zQG = select(gen_meas_Qg, [], 0)-20000; //发电机无功出力量测对应设备位置索引
measure_QG = select(gen_meas_Qg, [], 3);  //发电机无功出力量测值
sigma_PF=get(sigma_PF_mat,0,0); //支路首端有功潮流量测误差
sigma_PT=get(sigma_PT_mat,0,0); //支路尾端有功潮流量测误差
sigma_PG=get(sigma_PG_mat,0,0); //发电机有功出力量测误差
sigma_Vm=get(sigma_Vm_mat,0,0); //母线电压幅值量测误差
sigma_Va = []; //母线电压相角量测误差
sigma_QF = []; //支路首端无功潮流量测误差
sigma_QT = []; //支路尾端无功潮流量测误差
sigma_QG = []; //发电机无功出力量测误差
baseMVA=get(baseMVA_mat,0,0); //基准容量
areas = [1,1];  //所属区域
//通过URL导入状态估计计算所需函数
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/idx_gen.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/idx_bus.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/idx_brch.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/ext2int.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/int2ext.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/bustypes.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/make_y_bus.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/getV0.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/runse.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/dSbr_dV.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/dsbus_dv.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/doSE.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/isobservable.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/make_sdzip.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/make_sbus.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/pfsoln.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/total_load.txt
#include_url https://shufengdong.github.io/sparrowzz/rspower/lib/checkDataIntegrity.txt
type_initialguess = 2;  //平启动计算
[baseMVA, bus, gen, branch, success, z, z_est, error_sqrsum] = run_se(type_initialguess, []); //执行状态估计计算
return (bus, gen, branch);  //返回计算结果

该案例中通过设备id减去其编号基数获取设备位置索引。由于该节点计算结果包含多个返回参数，需要分别传给后续的边，还需要在界面对该节点填写列名为返回参数，并勾选自定义名称和Multi-returns，如下图所示：

（11）获取状态估计计算结果

通过边获取状态估计计算结果的bus、gen、branch返回参数，其中动作名称需与返回参数变量名称一致，如下图所示：

（12）状态估计结果写入测点的状态估计数据源

计算完成后，将状态估计写入测点的状态估计se数据源。以将电压幅值计算结果写入se数据源的节点和边配置如下图所示：

其中WriteFile类型的边用于将DataFrame数据写入量测数据源，DataFrame需包含point和analog/discrete2列，动作逻辑为：如果参数为meas://1, 表示用DataFrame更新数据源id为1的量测值，DataFrame包含point的列表示测点号，如果有discrete的列表示更新离散量，如果有analog的列表示更新模拟量，注意discrete或analog列不能同时出现。

上图在边5;31中选取point列，在边28;31选取计算结果列并将列名重命名为analog。在节点31中将从边获取的2个DataFrame拼接为1个DataFrame，此外，由于写入测点的值需为实数，通过字符串替换replace函数将计算值中的+0i虚数部分替换为空格，保留实部。

将状态估计其他结果写入se数据源的配置方法类似，如下图所示。

4. 状态估计结果查看

运行状态估计的DataFrame流，在事件-量测页面选择se数据源，可看到测点中写入的计算结果，如下图所示：

以下附件中提供了上述状态估计DataFrame流的配置文件，以及matpower状态估计3节点系统案例对应的的测点、设备、量测配置文件，可参考使用。

附件：电力系统状态估计案例.zip