Skip to content

S
sparrowzz
Commit ffb0e397 by Lian-runzhe
Options

状态估计及之前一些小错误修改

parent 6e5f2f04
正在显示 包含 1000 行增加 和 84 行删除
rspower/data/case3bus_P6_6.txt 0 → 100644
// CASE3BUS_P6_6 Case of 3 bus system.
// From Problem 6.6 in book 'Computational
// Methods for Electric Power Systems' by Mariesa Crow
// created by Rui Bo on 2007/11/12
// ----- Power Flow Data -----
// system MVA base
baseMVA = 1000;
// bus data
// bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin
bus = [
[1, 3, 350, 100, 0, 0, 1, 1, 0, 230, 1, 1.00, 1.00],
[2, 2, 400, 250, 0, 0, 1, 1, 0, 230, 1, 1.02, 1.02],
[3, 2, 250, 100, 0, 0, 1, 1, 0, 230, 1, 1.02, 1.02]
];
// generator data
// Note:
// 1 It's better of gen to be in number order, otherwise gen and genbid
// should be sorted to make the lp solution output clearly(in number order as well)
// 2 set Pmax to nonzero. set to 999 if no limit
// 3 If change the order of gen, then must change the order in genbid
// accordingly
// bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin
gen = [
[1, 182.18, 0, 999, -999, 1.00, 100, 1, 600, 0],
[2, 272.77, 0, 999, -999, 1.02, 100, 1, 400, 0],
[3, 545.05, 0, 999, -999, 1.02, 100, 1, 100, 0]
];
//gen = assign(gen, [[999],[999],[999]], [0], [8,9] ); % inactive the Pmax constraints
// branch data
// fbus tbus r x b rateA rateB rateC ratio angle status
branch = [
[1, 2, 0.01, 0.1, 0.050, 999, 100, 100, 0, 0, 1],
[1, 3, 0.05, 0.1, 0.025, 999, 100, 100, 0, 0, 1],
[2, 3, 0.05, 0.1, 0.025, 999, 100, 100, 0, 0, 1]
];
//----- OPF Data -----
// area data
areas = [
[1, 1]
];
// generator cost data
// 2 startup shutdown n c(n-1) ... c0
gencost = [
[2, 0, 0, 3, 1.5, 1, 0],
[2, 0, 0, 3, 1, 2, 0],
[2, 0, 0, 3, 0.5, 2.5, 0]
];
\ No newline at end of file
rspower/examples/test_dSbr_dV.txt 0 → 100644
#include ../data/case3bus_P6_6.txt
#include ../lib/idx_gen.txt
#include ../lib/idx_bus.txt
#include ../lib/idx_brch.txt
#include ../lib/make_y_bus.txt
#include ../lib/dSbr_dV.txt
V = [[1], [1.02], [1.02]];
[Ybus, Yf, Yt] = make_y_bus(baseMVA, bus, branch);
vcart = 0;
[dSf_dV1, dSf_dV2, dSt_dV1, dSt_dV2, Sf, St] = dSbr_dV(branch, Yf, Yt, V, vcart);
assert_relative_eq(dSf_dV1, [
[c(10.0990, -1.0099), c(-10.0990, 1.0099), c(0,0)],
[c(8.1600, -4.0800), c(0, 0), c(-8.1600, 4.0800)],
[c(0, 0), c(8.3232, -4.1616), c(-8.3232, 4.1616)]
],0.001);
assert_relative_eq(dSf_dV2, [
[c(0.9703, 9.6530), c(-0.9901, -9.9010), c(0,0)],
[c(3.9200, 7.8150), c(0, 0), c(-4.0000, -8.0000)],
[c(0, 0), c(4.0800, 8.1345), c(-4.0800, -8.1600)]
],0.001);
assert_relative_eq(dSt_dV1, [
[c(-10.0990, 1.0099), c(10.0990, -1.0099), c(0,0)],
[c(-8.1600, 4.0800), c(0, 0), c(8.1600, -4.0800)],
[c(0, 0), c(-8.3232, 4.1616), c(8.3232, -4.1616)]
],0.001);
assert_relative_eq(dSt_dV2, [
[c(-1.0099, -10.0990), c(1.0297, 10.2460), c(0,0)],
[c(-4.0800, -8.1600), c(0, 0), c(4.1600, 8.2945)],
[c(0, 0), c(-4.0800, -8.1600), c(4.0800, 8.1345)]
],0.001);
assert_relative_eq(Sf, [
[c(-0.0198, -0.2230)],
[c(-0.0800, -0.1725)],
[c(-4.4391e-16, -0.0130)]
],0.001);
assert_relative_eq(St, [
[c(0.0202,0.1760)],
[c(0.0816,0.1502)],
[c(4.4391e-16, -0.0130)]
],0.001);
rspower/examples/test_doSE.txt 0 → 100644
//同时测试isobservable()
#include ../data/case3bus_P6_6.txt
#include ../lib/idx_gen.txt
#include ../lib/idx_bus.txt
#include ../lib/idx_brch.txt
#include ../lib/isobservable.txt
#include ../lib/bustypes.txt
#include ../lib/make_y_bus.txt
#include ../lib/dSbr_dV.txt
#include ../lib/dSbus_dV.txt
#include ../lib/doSE.txt
[Ybus, Yf, Yt] = make_y_bus(baseMVA, bus, branch);
V0 = [[1], [1.02], [1.02]];
[ref,pv,pq] = bustypes(bus, gen);
//函数输入的自变量“measure, idx, sigma”matpower中为结构体类型，拆分的话变量太多了，故先不作为函数输入自变量
// which measurements are available
idx_zPF = [[1],[2]];
idx_zPT = [3];
idx_zPG = [[1],[2],[3]];
idx_zVa = [];
idx_zQF = [];
idx_zQT = [];
idx_zQG = [];
idx_zVm = [[2],[3]];
// specify measurements
measure_PF = [[0.12],[0.10]];
measure_PT = [-0.04];
measure_PG = [[0.58],[0.30],[0.14]];
measure_Va = [];
measure_QF = [];
measure_QT = [];
measure_QG = [];
measure_Vm = [[1.04],[0.98]];
// specify measurement variances
sigma_PF = 0.02;
sigma_PT = 0.02;
sigma_PG = 0.015;
sigma_Va = [];
sigma_QF = [];
sigma_QT = [];
sigma_QG = [];
sigma_Vm = 0.01;
[V, converged, iter_num, z, z_est, error_sqrsum] = doSE(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V0, ref, pv, pq);
assert_relative_eq(V, [
[c(1,0)],
[c(1.025603,-0.01745)],
[c(0.979010,0.000709)]
], 0.001);
assert_relative_eq(converged, 1);
assert_relative_eq(iter_num, 7);
assert_relative_eq(z, [
[0.1200],
[0.1000],
[-0.0400],
[0.5800],
[0.3000],
[0.1400],
[1.0400],
[0.9800]
]);
assert_relative_eq(z_est, [
[0.1474],
[0.0783],
[-0.0399],
[0.5757],
[0.3034],
[0.1336],
[1.0258],
[0.9790]
], 0.001);
assert_relative_eq(error_sqrsum, c(5.4179,-4.7371e-16), 0.001);
\ No newline at end of file
rspower/examples/test_getV0.txt 0 → 100644
#include ../data/case3bus_P6_6.txt
#include ../lib/idx_gen.txt
#include ../lib/idx_bus.txt
#include ../lib/idx_brch.txt
#include ../lib/getV0.txt
type_initialguess = 2;
V0_in = [[1],[2],[3]];
V0 = getV0(bus, gen, type_initialguess, V0_in);
assert_relative_eq(V0, [
[c(1,0)],
[c(1.02,0)],
[c(1.02,0)]
]);
\ No newline at end of file
rspower/examples/test_make_sdzip.txt
...@@ -20,4 +20,4 @@ assert_relative_eq(sdzip, [ ...@@ -20,4 +20,4 @@ assert_relative_eq(sdzip, [
[c(0,0), c(0,0), c(0.0610,0.0160)], [c(0,0), c(0,0), c(0.0610,0.0160)],
[c(0,0), c(0,0), c(0.1350,0.0580)], [c(0,0), c(0,0), c(0.1350,0.0580)],
[c(0,0), c(0,0), c(0.1490,0.0500)] [c(0,0), c(0,0), c(0.1490,0.0500)]
],4); ],0.001);
\ No newline at end of file \ No newline at end of file
rspower/examples/test_pfsoln.txt
...@@ -3,14 +3,14 @@ ...@@ -3,14 +3,14 @@
#include ../lib/idx_bus.txt #include ../lib/idx_bus.txt
#include ../lib/idx_brch.txt #include ../lib/idx_brch.txt
#include ../lib/make_y_bus.txt #include ../lib/make_y_bus.txt
#include ../lib/dsbus_dv.txt #include ../lib/dSbus_dV.txt
#include ../lib/make_sdzip.txt #include ../lib/make_sdzip.txt
#include ../lib/make_sbus.txt #include ../lib/make_sbus.txt
#include ../lib/newtonpf.txt #include ../lib/newtonpf.txt
#include ../lib/runpf.txt #include ../lib/runpf.txt
#include ../lib/pfsoln.txt #include ../lib/pfsoln.txt
#include ../lib/total_load.txt #include ../lib/total_load.txt
#include ../lib/bustypes.txt
info(`Running power flow on case14...`); info(`Running power flow on case14...`);
for i in 0..5 { for i in 0..5 {
...@@ -20,10 +20,8 @@ info(`Run power flow on case14 end`); ...@@ -20,10 +20,8 @@ info(`Run power flow on case14 end`);
[Ybus, Yf, Yt] = make_y_bus(baseMVA, bus, branch); [Ybus, Yf, Yt] = make_y_bus(baseMVA, bus, branch);
V = r; V = r;
bus_type = slice(bus, [0], [BUS_TYPE-1,BUS_TYPE]); [ref, pv, pq] = bustypes(bus, gen);
ref = find(bus_type == REF) + 1;
pv = find(bus_type == PV) + 1;
pq = find(bus_type == PQ) + 1;
[bus_output, gen_output, branch_output] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq, []); [bus_output, gen_output, branch_output] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq, []);
//bus结果 //bus结果
...@@ -42,7 +40,7 @@ assert_relative_eq(bus_output, [ ...@@ -42,7 +40,7 @@ assert_relative_eq(bus_output, [
[c(12,0), c(1,0), c(6.1,0.0), c(1.6,0.0), c(0.0,0.0), c(0.0,0.0), c(1.0,0.0), c(1.0552,0.0), c(-15.0756,0.0), c(0.0,0.0), c(1.0,0.0), c(1.06,0.0), c(0.94,0.0)], [c(12,0), c(1,0), c(6.1,0.0), c(1.6,0.0), c(0.0,0.0), c(0.0,0.0), c(1.0,0.0), c(1.0552,0.0), c(-15.0756,0.0), c(0.0,0.0), c(1.0,0.0), c(1.06,0.0), c(0.94,0.0)],
[c(13,0), c(1,0), c(13.5,0.0), c(5.8,0.0), c(0.0,0.0), c(0.0,0.0), c(1.0,0.0), c(1.0504,0.0), c(-15.1563,0.0), c(0.0,0.0), c(1.0,0.0), c(1.06,0.0), c(0.94,0.0)], [c(13,0), c(1,0), c(13.5,0.0), c(5.8,0.0), c(0.0,0.0), c(0.0,0.0), c(1.0,0.0), c(1.0504,0.0), c(-15.1563,0.0), c(0.0,0.0), c(1.0,0.0), c(1.06,0.0), c(0.94,0.0)],
[c(14,0), c(1,0), c(14.9,0.0), c(5.0,0.0), c(0.0,0.0), c(0.0,0.0), c(1.0,0.0), c(1.0355,0.0), c(-16.0336,0.0), c(0.0,0.0), c(1.0,0.0), c(1.06,0.0), c(0.94,0.0)] [c(14,0), c(1,0), c(14.9,0.0), c(5.0,0.0), c(0.0,0.0), c(0.0,0.0), c(1.0,0.0), c(1.0355,0.0), c(-16.0336,0.0), c(0.0,0.0), c(1.0,0.0), c(1.06,0.0), c(0.94,0.0)]
],4); ],0.01);
//gen结果 //gen结果
assert_relative_eq(gen_output, [ assert_relative_eq(gen_output, [
...@@ -51,7 +49,7 @@ assert_relative_eq(gen_output, [ ...@@ -51,7 +49,7 @@ assert_relative_eq(gen_output, [
[c(3.0,0.0), c(0.0,0.0), c(25.0753,0.0), c(40.0,0.0), c(0.0,0.0), c(1.0100,0.0), c(100.0,0.0), c(1.0,0.0), c(100.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0)], [c(3.0,0.0), c(0.0,0.0), c(25.0753,0.0), c(40.0,0.0), c(0.0,0.0), c(1.0100,0.0), c(100.0,0.0), c(1.0,0.0), c(100.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0)],
[c(6.0,0.0), c(0.0,0.0), c(12.7309,0.0), c(24.0,0.0), c(-6.0,-0.0), c(1.0700,0.0), c(100.0,0.0), c(1.0,0.0), c(100.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0)], [c(6.0,0.0), c(0.0,0.0), c(12.7309,0.0), c(24.0,0.0), c(-6.0,-0.0), c(1.0700,0.0), c(100.0,0.0), c(1.0,0.0), c(100.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0)],
[c(8.0,0.0), c(0.0,0.0), c(17.6235,0.0), c(24.0,0.0), c(-6.0,-0.0), c(1.0900,0.0), c(100.0,0.0), c(1.0,0.0), c(100.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0)] [c(8.0,0.0), c(0.0,0.0), c(17.6235,0.0), c(24.0,0.0), c(-6.0,-0.0), c(1.0900,0.0), c(100.0,0.0), c(1.0,0.0), c(100.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0), c(0.0,0.0)]
],4); ],0.01);
//branch结果 //branch结果
assert_relative_eq(branch_output, [ assert_relative_eq(branch_output, [
...@@ -75,4 +73,4 @@ assert_relative_eq(branch_output, [ ...@@ -75,4 +73,4 @@ assert_relative_eq(branch_output, [
[c(10.0, 0.0), c(11.0, 0.0), c(0.0821, 0.0), c(0.1921, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(1.0, 0.0), c(-360.0, -0.0), c(360.0, 0.0), c(-3.7853, 0.0), c(-1.6151, 0.0), c(3.7979, 0.0), c(1.6445, 0.0)], [c(10.0, 0.0), c(11.0, 0.0), c(0.0821, 0.0), c(0.1921, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(1.0, 0.0), c(-360.0, -0.0), c(360.0, 0.0), c(-3.7853, 0.0), c(-1.6151, 0.0), c(3.7979, 0.0), c(1.6445, 0.0)],
[c(12.0, 0.0), c(13.0, 0.0), c(0.2209, 0.0), c(0.1999, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(1.0, 0.0), c(-360.0, -0.0), c(360.0, 0.0), c(1.6143, 0.0), c(0.7540, 0.0), c(-1.6080, 0.0), c(-0.7483, 0.0)], [c(12.0, 0.0), c(13.0, 0.0), c(0.2209, 0.0), c(0.1999, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(1.0, 0.0), c(-360.0, -0.0), c(360.0, 0.0), c(1.6143, 0.0), c(0.7540, 0.0), c(-1.6080, 0.0), c(-0.7483, 0.0)],
[c(13.0, 0.0), c(14.0, 0.0), c(0.1709, 0.0), c(0.3480, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(1.0, 0.0), c(-360.0, -0.0), c(360.0, 0.0), c(5.6439, 0.0), c(1.7472, 0.0), c(-5.5898, 0.0), c(-1.6371, 0.0)] [c(13.0, 0.0), c(14.0, 0.0), c(0.1709, 0.0), c(0.3480, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(0.0, 0.0), c(1.0, 0.0), c(-360.0, -0.0), c(360.0, 0.0), c(5.6439, 0.0), c(1.7472, 0.0), c(-5.5898, 0.0), c(-1.6371, 0.0)]
],4); ],0.01);
rspower/examples/test_runpf.txt
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
#include ../lib/idx_bus.txt #include ../lib/idx_bus.txt
#include ../lib/idx_brch.txt #include ../lib/idx_brch.txt
#include ../lib/make_y_bus.txt #include ../lib/make_y_bus.txt
#include ../lib/dsbus_dv.txt #include ../lib/dSbus_dV.txt
#include ../lib/make_sdzip.txt #include ../lib/make_sdzip.txt
#include ../lib/make_sbus.txt #include ../lib/make_sbus.txt
#include ../lib/newtonpf.txt #include ../lib/newtonpf.txt
...@@ -30,4 +30,4 @@ assert_relative_eq(r, [ ...@@ -30,4 +30,4 @@ assert_relative_eq(r, [
[c(1.0189,-0.2744)], [c(1.0189,-0.2744)],
[c(1.0138,-0.2746)], [c(1.0138,-0.2746)],
[c(0.9952,-0.2860)] [c(0.9952,-0.2860)]
],4); ], 0.001);
\ No newline at end of file \ No newline at end of file
rspower/examples/test_runse.txt 0 → 100644
#include ../data/case3bus_P6_6.txt
#include ../lib/idx_gen.txt
#include ../lib/idx_bus.txt
#include ../lib/idx_brch.txt
#include ../lib/ext2int.txt
#include ../lib/int2ext.txt
#include ../lib/bustypes.txt
#include ../lib/make_y_bus.txt
#include ../lib/getV0.txt
#include ../lib/runse.txt
#include ../lib/dSbr_dV.txt
#include ../lib/dSbus_dV.txt
#include ../lib/doSE.txt
#include ../lib/isobservable.txt
#include ../lib/make_sdzip.txt
#include ../lib/make_sbus.txt
#include ../lib/pfsoln.txt
#include ../lib/total_load.txt
//输入观测值
// which measurements are available
idx_zPF = [[1],[2]];
idx_zPT = [3];
idx_zPG = [[1],[2],[3]];
idx_zVa = [];
idx_zQF = [];
idx_zQT = [];
idx_zQG = [];
idx_zVm = [[2],[3]];
// specify measurements
measure_PF = [[0.12],[0.10]];
measure_PT = [-0.04];
measure_PG = [[0.58],[0.30],[0.14]];
measure_Va = [];
measure_QF = [];
measure_QT = [];
measure_QG = [];
measure_Vm = [[1.04],[0.98]];
// specify measurement variances
sigma_PF = 0.02;
sigma_PT = 0.02;
sigma_PG = 0.015;
sigma_Va = [];
sigma_QF = [];
sigma_QT = [];
sigma_QG = [];
sigma_Vm = 0.01;
type_initialguess = 2;
V0 = [[1], [1.02], [1.02]];
[baseMVA, bus, gen, branch, success, z, z_est, error_sqrsum] = run_se(type_initialguess, V0);
assert_relative_eq(baseMVA, 1000);
assert_relative_eq(bus, [
[1, 3, 350, 100, 0, 0, 1, 1, 0, 230, 1, 1, 1],
[2, 2, 400, 250, 0, 0, 1, 1.0258, -0.9748, 230, 1, 1.02, 1.02],
[3, 2, 250, 100, 0, 0, 1, 0.9790, 0.0415, 230, 1, 1.02, 1.02]
], 0.001);
assert_relative_eq(gen, [
[1, 575.7215, -37.5179, 999, -999, 1, 100, 1, 600, 0],
[2, 303.4183, 946.8882, 999, -999, 1.02, 100, 1, 400, 0],
[3, 133.5777, -527.2488, 999, -999, 1.02, 100, 1, 100, 0]
], 0.001);
assert_relative_eq(branch, [
[1, 2, 0.01, 0.1, 0.05, 999, 100, 100, 0, 0, 1, 147.4384, -295.7713, -146.4878, 253.9727],
[1, 3, 0.05, 0.1, 0.025, 999, 100, 100, 0, 0, 1, 78.2831, 158.2534, -76.5189, -179.2057],
[2, 3, 0.05, 0.1, 0.025, 999, 100, 100, 0, 0, 1, 49.9061, 442.9156, -39.9035, -448.0432]
], 0.001);
assert_relative_eq(success, 1);
assert_relative_eq(z, [
[0.1200],
[0.1000],
[-0.0400],
[0.5800],
[0.3000],
[0.1400],
[1.0400],
[0.9800]
]);
assert_relative_eq(z_est, [
[0.1474],
[0.0783],
[-0.0399],
[0.5757],
[0.3034],
[0.1336],
[1.0258],
[0.9790]
], 0.001);
assert_relative_eq(error_sqrsum, c(5.4179,-4.7371e-16), 0.001);
\ No newline at end of file
rspower/lib/bustypes.txt
...@@ -14,7 +14,8 @@ ...@@ -14,7 +14,8 @@
// Covered by the 3-clause BSD License (see LICENSE file for details). // Covered by the 3-clause BSD License (see LICENSE file for details).
// See https://matpower.org for more info. // See https://matpower.org for more info.
fn bustypes_ref(bus, gen) { fn bustypes(bus, gen) {
//ref, pv, pq索引从'0'开始
// get generator status // get generator status
nb = size(bus, 0); nb = size(bus, 0);
ng = size(gen, 0); ng = size(gen, 0);
...@@ -25,45 +26,9 @@ fn bustypes_ref(bus, gen) { ...@@ -25,45 +26,9 @@ fn bustypes_ref(bus, gen) {
// form index lists for slack, PV, and PQ buses // form index lists for slack, PV, and PQ buses
bus_type = slice(bus, [0], [BUS_TYPE-1,BUS_TYPE]); // bus type column bus_type = slice(bus, [0], [BUS_TYPE-1,BUS_TYPE]); // bus type column
ref = find(bus_type == REF && bus_gen_status); // reference bus index ref = find(bus_type == REF && bus_gen_status)'; // reference bus index
pv = find(bus_type == PV && bus_gen_status); // PV bus indices pv = find(bus_type == PV && bus_gen_status)'; // PV bus indices
pq = find(bus_type == PQ || ~~bus_gen_status); // PQ bus indices pq = find(bus_type == PQ || ~~bus_gen_status)'; // PQ bus indices
return ref; return(ref, pv, pq);
}
fn bustypes_pv(bus, gen) {
// get generator status
nb = size(bus, 0);
ng = size(gen, 0);
g_i = slice(gen, [0], [GEN_BUS-1,GEN_BUS]) - 1; // generator bus indices
Cg = full(sparse(g_i, range(0, ng), slice(gen, [0], [GEN_STATUS-1,GEN_STATUS]) > 0, nb, ng)); // gen connection matrix
// element i, j is 1 if, generator j at bus i is ON
bus_gen_status = Cg * ones(ng, 1); // number of generators at each bus that are ON
// form index lists for slack, PV, and PQ buses
bus_type = slice(bus, [0], [BUS_TYPE-1,BUS_TYPE]); // bus type column
ref = find(bus_type == REF && bus_gen_status); // reference bus index
pv = find(bus_type == PV && bus_gen_status); // PV bus indices
pq = find(bus_type == PQ || ~~bus_gen_status); // PQ bus indices
return pv;
}
fn bustypes_pq(bus, gen) {
// get generator status
nb = size(bus, 0);
ng = size(gen, 0);
g_i = slice(gen, [0], [GEN_BUS-1,GEN_BUS]) - 1; // generator bus indices
Cg = full(sparse(g_i, range(0, ng), slice(gen, [0], [GEN_STATUS-1,GEN_STATUS]) > 0, nb, ng)); // gen connection matrix
// element i, j is 1 if, generator j at bus i is ON
bus_gen_status = Cg * ones(ng, 1); // number of generators at each bus that are ON
// form index lists for slack, PV, and PQ buses
bus_type = slice(bus, [0], [BUS_TYPE-1,BUS_TYPE]); // bus type column
ref = find(bus_type == REF && bus_gen_status); // reference bus index
pv = find(bus_type == PV && bus_gen_status); // PV bus indices
pq = find(bus_type == PQ || ~~bus_gen_status); // PQ bus indices
return pq;
} }
\ No newline at end of file
rspower/lib/dSbr_dV.txt 0 → 100644
//DSBR_DV Computes partial derivatives of branch power flows w.r.t. voltage.
//
// The derivatives can be take with respect to polar or cartesian coordinates
// of voltage, depending on the 5th argument.
//
// [DSF_DVA, DSF_DVM, DST_DVA, DST_DVM, SF, ST] = DSBR_DV(BRANCH, YF, YT, V)
// [DSF_DVA, DSF_DVM, DST_DVA, DST_DVM, SF, ST] = DSBR_DV(BRANCH, YF, YT, V, 0)
//
// Returns four matrices containing partial derivatives of the complex
// branch power flows at "from" and "to" ends of each branch w.r.t voltage
// magnitude and voltage angle, respectively (for all buses).
//
// [DSF_DVR, DSF_DVI, DST_DVR, DST_DVI, SF, ST] = DSBR_DV(BRANCH, YF, YT, V, 1)
//
// Returns four matrices containing partial derivatives of the complex
// branch power flows at "from" and "to" ends of each branch w.r.t real and
// imaginary parts of voltage, respectively (for all buses).
//
// If YF is a sparse matrix, the partial derivative matrices will be as well.
// Optionally returns vectors containing the power flows themselves. The
// following explains the expressions used to form the matrices:
//
// If = Yf * V;
// Sf = diag(Vf) * conj(If) = diag(conj(If)) * Vf
//
// Polar coordinates:
// Partials of V, Vf & If w.r.t. voltage angles
// dV/dVa = j * diag(V)
// dVf/dVa = sparse(1:nl, f, j * V(f)) = j * sparse(1:nl, f, V(f))
// dIf/dVa = Yf * dV/dVa = Yf * j * diag(V)
//
// Partials of V, Vf & If w.r.t. voltage magnitudes
// dV/dVm = diag(V./abs(V))
// dVf/dVm = sparse(1:nl, f, V(f)./abs(V(f))
// dIf/dVm = Yf * dV/dVm = Yf * diag(V./abs(V))
//
// Partials of Sf w.r.t. voltage angles
// dSf/dVa = diag(Vf) * conj(dIf/dVa)
// + diag(conj(If)) * dVf/dVa
// = diag(Vf) * conj(Yf * j * diag(V))
// + conj(diag(If)) * j * sparse(1:nl, f, V(f))
// = -j * diag(Vf) * conj(Yf * diag(V))
// + j * conj(diag(If)) * sparse(1:nl, f, V(f))
// = j * (conj(diag(If)) * sparse(1:nl, f, V(f))
// - diag(Vf) * conj(Yf * diag(V)))
//
// Partials of Sf w.r.t. voltage magnitudes
// dSf/dVm = diag(Vf) * conj(dIf/dVm)
// + diag(conj(If)) * dVf/dVm
// = diag(Vf) * conj(Yf * diag(V./abs(V)))
// + conj(diag(If)) * sparse(1:nl, f, V(f)./abs(V(f)))
//
// Cartesian coordinates:
// Partials of V, Vf & If w.r.t. real part of complex voltage
// dV/dVr = diag(ones(n,1))
// dVf/dVr = Cf
// dIf/dVr = Yf
// where Cf is the connection matrix for line & from buses
//
// Partials of V, Vf & If w.r.t. imaginary part of complex voltage
// dV/dVi = j * diag(ones(n,1))
// dVf/dVi = j * Cf
// dIf/dVi = j * Yf
//
// Partials of Sf w.r.t. real part of complex voltage
// dSf/dVr = conj(diag(If)) * Cf + diag(Vf) * conj(Yf)
//
// Partials of Sf w.r.t. imaginary part of complex voltage
// dSf/dVi = j * (conj(diag(If)) * Cf - diag(Vf) * conj(Yf))
//
// Derivations for "to" bus are similar.
//
// Examples:
// [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
// [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = ...
// dSbr_dV(branch, Yf, Yt, V);
// [dSf_dVr, dSf_dVi, dSt_dVr, dSt_dVi, Sf, St] = ...
// dSbr_dV(branch, Yf, Yt, V, 1);
// 计算支路功率流对电压的偏导数。
fn dSbr_dV(branch, Yf, Yt, V, vcart) {
// 默认输入参数
if is_empty(vcart) {
vcart = 0; // 默认为极坐标
}
// 定义
f = slice(branch, [0], [F_BUS - 1, F_BUS]) - 1; // "from" 母线列表(f索引从'0'开始)
t = slice(branch, [0], [T_BUS - 1, T_BUS]) - 1; // "to" 母线列表(t索引从'0'开始)
nl = size(f, 0);
nb = size(V, 0);
// 计算中间值
Yfc = conj(Yf);
Ytc = conj(Yt);
Vc = conj(V);
Ifc = Yfc * Vc; // "from" 电流的共轭
Itc = Ytc * Vc; // "to" 电流的共轭
diagVf = diag(get_multi(V, f));
diagVt = diag(get_multi(V, t));
diagIfc = diag(Ifc);
diagItc = diag(Itc);
if vcart == 0 { // 极坐标
Vnorm = V ./ abs(V);
diagVc = diag(Vc);
diagVnorm = diag(Vnorm);
CVf = full(sparse(range(0, nl), f, get_multi(V, f), nl, nb));
CVnf = full(sparse(range(0, nl), f, get_multi(Vnorm, f), nl, nb));
CVt = full(sparse(range(0, nl), t, get_multi(V, t), nl, nb));
CVnt = full(sparse(range(0, nl), t, get_multi(Vnorm, t), nl, nb));
}
if vcart { // 直角坐标
Cf = sparse(range(0, nl), f, ones(nl, 1), nl, nb); // 连接矩阵 for line & from buses
Ct = sparse(range(0, nl), t, ones(nl, 1), nl, nb); // 连接矩阵 for line & to buses
Af = diagIfc * Cf;
Bf = diagVf * Yfc;
At = diagItc * Ct;
Bt = diagVt * Ytc;
dSf_dV1 = Af + Bf; // dSf_dVr
dSf_dV2 = c(0, 1) * (Af - Bf); // dSf_dVi
dSt_dV1 = At + Bt; // dSt_dVr
dSt_dV2 = c(0, 1) * (At - Bt); // dSt_dVi
} else { // 极坐标
dSf_dV1 = c(0, 1) * (diagIfc * CVf - diagVf * Yfc * diagVc); // dSf_dVa
dSf_dV2 = diagVf * conj(Yf * diagVnorm) + diagIfc * CVnf; // dSf_dVm
dSt_dV1 = c(0, 1) * (diagItc * CVt - diagVt * Ytc * diagVc); // dSt_dVa
dSt_dV2 = diagVt * conj(Yt * diagVnorm) + diagItc * CVnt; // dSt_dVm
}
// 返回功率流本身
Sf = get_multi(V, f)' .* Ifc;
St = get_multi(V, t)' .* Itc;
/*if is_sparse(Yf) { // 稀疏矩阵版本 (如果Yf是稀疏矩阵)
diagVf = sparse(range(0, nl), range(0, nl), get_multi(V, f), nl, nl);
diagVt = sparse(range(0, nl), range(0, nl), get_multi(V, t), nl, nl);
diagIfc = sparse(range(0, nl), range(0, nl), Ifc, nl, nl);
diagItc = sparse(range(0, nl), range(0, nl), Itc, nl, nl);
if vcart == 0 { // 极坐标
Vnorm = V ./ abs(V);
diagVc = sparse(range(0, nb), range(0, nb), Vc, nb, nb);
diagVnorm = sparse(range(0, nb), range(0, nb), Vnorm, nb, nb);
CVf = sparse(range(0, nl), f, get_multi(V, f), nl, nb);
CVnf = sparse(range(0, nl), f, get_multi(get_multi(Vnorm, f), range(0, size(f, 0))), nl, nb);
CVt = sparse(range(0, nl), t, get_multi(V, t), nl, nb);
CVnt = sparse(range(0, nl), t, get_multi(get_multi(Vnorm, t), range(0, size(t, 0))), nl, nb);
}
} else { // 稠密矩阵版本
diagVf = diag(get_multi(V, f));
diagVt = diag(get_multi(V, t));
diagIfc = diag(Ifc);
diagItc = diag(Itc);
if vcart == 0 { // 极坐标
Vnorm = V ./ abs(V);
diagVc = diag(Vc);
diagVnorm = diag(Vnorm);
CVf = full(sparse(range(0, nl), f, get_multi(V, f), nl, nb));
CVnf = full(sparse(range(0, nl), f, get_multi(Vnorm, f), nl, nb));
CVt = full(sparse(range(0, nl), t, get_multi(V, t), nl, nb));
CVnt = full(sparse(range(0, nl), t, get_multi(Vnorm, t), nl, nb));
}
}
*/
return (dSf_dV1, dSf_dV2, dSt_dV1, dSt_dV2, Sf, St);
}
rspower/lib/doSE.txt 0 → 100644
// 状态估计主函数
fn doSE(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V0, ref, pv, pq) {
// 迭代参数设置
tol = 1e-5; // 收敛 tolerance
max_it = 100; // 最大迭代次数
verbose = 0; // 输出控制
// 初始化变量
converged = 0;
iter_num = 0;
V = V0;
Va = angle(V); // 电压角度（弧度）
Vm = abs(V); // 电压幅值
nb = size(Ybus, 0); // 母线数量
f = slice(branch, [0], [F_BUS - 1, F_BUS]) - 1; // 支路"从"母线索引(f索引从'0'开始)
t = slice(branch, [0], [T_BUS - 1, T_BUS]) - 1; // 支路"到"母线索引(t索引从'0'开始)
// 非平衡节点（PV+PQ）
nonref = vertcat(pv, pq); //(pv, pq索引从'0'开始)
// ---------------------- 构建测量向量 z ----------------------
// 测量类型顺序：PF -> PT -> PG -> Va -> QF -> QT -> QG -> Vm
z = vertcat(
measure_PF,
measure_PT,
measure_PG,
measure_Va,
measure_QF,
measure_QT,
measure_QG,
measure_Vm
);
// ---------------------- 构建权重矩阵 R_inv ----------------------
// 构建sigma向量（与z维度匹配）
sigma_vector = vertcat(
sigma_PF * ones(size(idx_zPF, 0), 1),
sigma_PT * ones(size(idx_zPT, 0), 1),
sigma_PG * ones(size(idx_zPG, 0), 1),
sigma_Va * ones(size(idx_zVa, 0), 1),
sigma_QF * ones(size(idx_zQF, 0), 1),
sigma_QT * ones(size(idx_zQT, 0), 1),
sigma_QG * ones(size(idx_zQG, 0), 1),
sigma_Vm * ones(size(idx_zVm, 0), 1)
);
sigma_square = sigma_vector .^ 2;
R_inv = diag(1.0 ./ sigma_square); // 权重矩阵（对角阵）
// ---------------------- 牛顿迭代主循环 ----------------------
while (~~converged && iter_num < max_it) {
iter_num = iter_num + 1;
// ---------------------- 计算估计测量值 z_est ----------------------
// 支路潮流计算
Sfe = get_multi(V, f)' .* conj(Yf * V); // 从端潮流
Ste = get_multi(V, t)' .* conj(Yt * V); // 到端潮流
// 发电机节点注入功率计算
gbus = slice(gen, [0], [GEN_BUS - 1, GEN_BUS]) - 1; // 发电机所在母线(索引从“0”开始)
Sgbus = get_multi(V, gbus)' .* conj(select(Ybus, gbus) * V); // 母线注入功率
// 发电机输出功率 = 注入功率 + 本地负荷（标幺值）
Sgen = Sgbus * baseMVA + (select(bus, gbus, PD-1) + c(0,1) * select(bus, gbus, QD-1));
Sgen = Sgen / baseMVA;
// ---------------------- 计算雅可比矩阵 H ----------------------
// 节点注入功率雅可比
[dSbus_dVa, dSbus_dVm] = dSbus_dV(Ybus, V, 0);
// 支路潮流雅可比
[dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V, 0);
genbus_row = gbus;
// 雅可比子矩阵提取（仅非平衡节点）
// 1. 有功潮流雅可比
dPF_dVa = real(dSf_dVa); // 从端有功对角度的偏导
dPF_dVm = real(dSf_dVm); // 从端有功对幅值的偏导
dPT_dVa = real(dSt_dVa); // 到端有功对角度的偏导
dPT_dVm = real(dSt_dVm); // 到端有功对幅值的偏导
// 2. 无功潮流雅可比
dQF_dVa = imag(dSf_dVa); // 从端无功对角度的偏导
dQF_dVm = imag(dSf_dVm); // 从端无功对幅值的偏导
dQT_dVa = imag(dSt_dVa); // 到端无功对角度的偏导
dQT_dVm = imag(dSt_dVm); // 到端无功对幅值的偏导
// 3. 发电机功率雅可比
dPG_dVa = real(select(dSbus_dVa, genbus_row)); // 有功对角度的偏导
dPG_dVm = real(select(dSbus_dVm, genbus_row)); // 有功对幅值的偏导
dQG_dVa = imag(select(dSbus_dVa, genbus_row)); // 无功对角度的偏导
dQG_dVm = imag(select(dSbus_dVm, genbus_row)); // 无功对幅值的偏导
// 4. 电压角度雅可比（单位矩阵）
dVa_dVa = eye(nb);
dVa_dVm = zeros(nb, nb);
// 5. 电压幅值雅可比（单位矩阵）
dVm_dVa = zeros(nb, nb);
dVm_dVm = eye(nb);
// ---------------------- 构建估计测量向量 z_est和雅可比矩阵 H ----------------------
//有功潮流（从端）
if ~~is_empty(idx_zPF){
z_est1 = real(select(Sfe, idx_zPF - 1));
H_1 = horzcat(select(dPF_dVa, idx_zPF - 1, nonref), select(dPF_dVm, idx_zPF - 1, nonref));
}else{
z_est1 = [];
H_1 = [];
}
//有功潮流（到端）
if ~~is_empty(idx_zPT){
z_est2 = real(select(Ste, idx_zPT - 1));
H_2 = horzcat(select(dPT_dVa, idx_zPT - 1, nonref), select(dPT_dVm, idx_zPT - 1, nonref));
}else{
z_est2 = [];
H_2 = [];
}
//发电机有功输出
if ~~is_empty(idx_zPG){
z_est3 = real(select(Sgen, idx_zPG - 1));
H_3 = horzcat(select(dPG_dVa, idx_zPG - 1, nonref), select(dPG_dVm, idx_zPG - 1, nonref));
}else{
z_est3 = [];
H_3 = [];
}
//母线电压角度
if ~~is_empty(idx_zVa){
z_est4 = select(angle(V), idx_zVa - 1);
H_4 = horzcat(select(dVa_dVa, idx_zVa - 1, nonref), select(dVa_dVm, idx_zVa - 1, nonref));
}else{
z_est4 = [];
H_4 = [];
}
//无功潮流（从端）
if ~~is_empty(idx_zQF){
z_est5 = imag(select(Sfe, idx_zQF - 1));
H_5 = horzcat(select(dQF_dVa, idx_zQF - 1, nonref), select(dQF_dVm, idx_zQF - 1, nonref));
}else{
z_est5 = [];
H_5 = [];
}
//无功潮流（到端）
if ~~is_empty(idx_zQT){
z_est6 = imag(select(Ste, idx_zQT - 1));
H_6 = horzcat(select(dQT_dVa, idx_zQT - 1, nonref), select(dQT_dVm, idx_zQT - 1, nonref));
}else{
z_est6 = [];
H_6 = [];
}
//发电机无功输出
if ~~is_empty(idx_zQG){
z_est7 = imag(select(Sgen, idx_zQG - 1));
H_7 = horzcat(select(dQG_dVa, idx_zQG - 1, nonref), select(dQG_dVm, idx_zQG - 1, nonref));
}else{
z_est7 = [];
H_7 = [];
}
//母线电压幅值
if ~~is_empty(idx_zVm){
z_est8 = select(abs(V), idx_zVm - 1);
H_8 = horzcat(select(dVm_dVa, idx_zVm - 1, nonref), select(dVm_dVm, idx_zVm - 1, nonref));
}else{
z_est8 = [];
H_8 = [];
}
// 构建估计测量向量 z_est
z_est = vertcat(z_est1,z_est2,z_est3,z_est4,z_est5,z_est6,z_est7,z_est8);
// 拼接完整雅可比矩阵 H
H = vertcat(H_1,H_2,H_3,H_4,H_5,H_6,H_7,H_8);
// ---------------------- 计算修正量 dx ----------------------
J = H' * R_inv * H; // 信息矩阵
mismatch = z - z_est; // 测量残差
F = H' * R_inv * mismatch; // 牛顿方向
// 可观测性检查
TorF = isobservable(H, pv, pq);
if ~~TorF {
println(`doSE: 系统不可观测`);
}
dx = J \ F; // 求解修正方程
// ---------------------- 收敛判断 ----------------------
normF = max(abs(F)); // 残差无穷范数
if verbose > 1 {
println("\niteration [%3d]\t\tnorm of mismatch: %10.3e", iter_num, normF);
}
if normF < tol {
converged = 1;
}
// ---------------------- 更新电压 ----------------------
// 角度更新（非平衡节点）
Va = set(Va, nonref, select(Va,nonref) + slice(dx, [0,size(nonref, 0)], [0]));
// 幅值更新（非平衡节点）
Vm = set(Vm, nonref, select(Vm,nonref) + slice(dx, [size(nonref, 0),2*size(nonref, 0)], [0]));
// 重构电压向量（角度为弧度）
V = Vm .* exp(c(0,1) * Va);
// 重新计算幅值和角度（避免负幅值问题）
Vm = abs(V);
Va = angle(V);
}
// ---------------------- 计算加权残差平方和 ----------------------
error_sqrsum = sum_all((z - z_est) .^ 2 ./ sigma_square);
// 返回结果（多返回值）
return (V, converged, iter_num, z, z_est, error_sqrsum);
}
\ No newline at end of file
rspower/lib/dsbus_dv.txt
...@@ -79,16 +79,12 @@ ...@@ -79,16 +79,12 @@
// [Online]. Available: https://matpower.org/docs/TN4-OPF-Derivatives-Cartesian.pdf // [Online]. Available: https://matpower.org/docs/TN4-OPF-Derivatives-Cartesian.pdf
// doi: 10.5281/zenodo.3237909 // doi: 10.5281/zenodo.3237909
// MATPOWER
// Copyright (c) 1996-2019, Power Systems Engineering Research Center (PSERC)
// by Ray Zimmerman, PSERC Cornell
// and Baljinnyam Sereeter, Delft University of Technology
//
// This file is part of MATPOWER.
// Covered by the 3-clause BSD License (see LICENSE file for details).
// See https://matpower.org for more info.
fn dSbus_dV(Ybus, V, vcart) { fn dSbus_dV(Ybus, V, vcart) {
// 默认输入参数
if is_empty(vcart) {
vcart = 0; // 默认为极坐标
}
n = length(V); n = length(V);
Ibus = Ybus * V; Ibus = Ybus * V;
diagV = diag(V); diagV = diag(V);
...@@ -102,5 +98,5 @@ fn dSbus_dV(Ybus, V, vcart) { ...@@ -102,5 +98,5 @@ fn dSbus_dV(Ybus, V, vcart) {
dSbus_dV1 = c(0,1) * diagV * conj(diagIbus - Ybus * diagV); // dSbus/dVa dSbus_dV1 = c(0,1) * diagV * conj(diagIbus - Ybus * diagV); // dSbus/dVa
dSbus_dV2 = diagV * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm; // dSbus/dVm dSbus_dV2 = diagV * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm; // dSbus/dVm
} }
return horzcat(dSbus_dV1, dSbus_dV2); return (dSbus_dV1, dSbus_dV2);
} }
\ No newline at end of file
rspower/lib/getV0.txt 0 → 100644
// GETV0 Get initial voltage profile for power flow calculation.
// Note: The pv bus voltage will remain at the given value even for
// flat start.
// type_initialguess: 1 - initial guess from case data
// 2 - flat start
// 3 - from input
// GETV0 获取潮流计算的初始电压分布
// 注意：即使是平启动，PV节点电压也会保持给定值
// type_initialguess: 1 - 从案例数据获取初始猜测值
// 2 - 平启动
// 3 - 从输入获取
fn getV0(bus, gen, type_initialguess, V0_in) {
// 发电机信息
on = find(slice(gen, [0], [GEN_STATUS - 1, GEN_STATUS]) > 0); //查找运行中的发电机索引（这里索引从0开始）
gbus = get_multi(slice(gen, [0], [GEN_BUS-1, GEN_BUS]), on); // 转换为一维张量（母线编号，注意是1基）
gbus_idx = gbus - 1; // 转换为0基索引
if type_initialguess == 1 {
// 使用案例数据中的先前值
// 注意：案例数据中的角度是度，而潮流求解器中是弧度，因此需要转换
vm_col = slice(bus, [0], [VM - 1, VM]); // 获取电压幅值列（二维张量）
vm = slice(vm_col, [0], 0); // 转换为一维张量
va_col = slice(bus, [0], [VA - 1, VA]); // 获取电压角度列（二维张量）
va_deg = slice(va_col, [0], 0); // 转换为一维张量（度）
va_rad = va_deg * pi / 180; // 转换为弧度
V0 = vm .* exp(c(0, 1) .* va_rad); // 计算初始电压（复数）
} else if type_initialguess == 2 {
// 平启动
nb = size(bus, 0); // 母线数量
V0 = ones(nb, 1); // 初始电压幅值为1，角度为0（复数形式为1+0i）
} else if type_initialguess == 3 {
// 使用给定的初始电压
V0 = V0_in;
} else {
//info(`Error: unknown 'type_initialguess'`);
return [];
}
// 将PV节点和参考节点的电压设置到初始猜测值中
gen_vg_col = select(gen, on, [VG - 1, VG]); // 获取运行中发电机的电压幅值设定值列（二维张量）
gen_vg = slice(gen_vg_col, [0], 0); // 转换为一维张量
v0_gbus = get_multi(V0, gbus_idx); // 获取初始猜测中发电机所在母线的电压
new_v0 = gen_vg .* v0_gbus ./ abs(v0_gbus); // 新的电压（设定幅值，保持原有角度）
V0 = set(V0, gbus_idx, new_v0); // 更新初始电压
return V0;
}
\ No newline at end of file
rspower/lib/int2ext.txt 0 → 100644
fn int2ext(i2e, bus, gen, branch, areas) {
//INT2EXT Converts internal to external bus numbering.
//
// This function has two forms, (1) the old form that operates on
// and returns individual matrices and (2) the new form that operates
// on and returns an entire MATPOWER case struct.
//
// 1. [BUS, GEN, BRANCH, AREAS] = INT2EXT(I2E, BUS, GEN, BRANCH, AREAS)
// [BUS, GEN, BRANCH] = INT2EXT(I2E, BUS, GEN, BRANCH)
//
// Converts from the consecutive internal bus numbers back to the originals
// using the mapping provided by the I2E vector returned from EXT2INT,
// where EXTERNAL_BUS_NUMBER = I2E(INTERNAL_BUS_NUMBER).
// AREAS is completely ignored and is only included here for backward
// compatibility of the API.
//
// Examples:
// [bus, gen, branch, areas] = int2ext(i2e, bus, gen, branch, areas);
// [bus, gen, branch] = int2ext(i2e, bus, gen, branch);
//
// 2. MPC = INT2EXT(MPC)
// MPC = INT2EXT(MPC, MPOPT)
//
// If the input is a single MATPOWER case struct, followed optionally
// by a MATOWER options struct, then it restores all buses, generators
// and branches that were removed because of being isolated or off-line,
// and reverts to the original generator ordering and original bus
// numbering. This requires that the 'order' field created by EXT2INT be
// in place.
//
// Examples:
// mpc = int2ext(mpc);
// mpc = int2ext(mpc, mpopt);
// 重新编号母线
index_i2e_bus = slice(bus, [0], [BUS_I-1, BUS_I]) - 1; // 索引值=内部母线编号-1
i2e_bus = get_multi(i2e, index_i2e_bus)'; // 获取对应的外部母线编号
output_bus = assign(bus, i2e_bus, [0], [BUS_I-1, BUS_I]); // 转换为外部母线编号
// 重新编号发电机连接的母线
index_i2e_gen = slice(gen, [0], [GEN_BUS-1, GEN_BUS]) - 1;
i2e_gen = get_multi(i2e, index_i2e_gen)';
output_gen = assign(gen, i2e_gen, [0], [GEN_BUS-1, GEN_BUS]);
// 重新编号支路的起始母线
index_i2e_fbus = slice(branch, [0], [F_BUS-1, F_BUS]) - 1;
i2e_fbus = get_multi(i2e, index_i2e_fbus)';
output_branch1 = assign(branch, i2e_fbus, [0], [F_BUS-1, F_BUS]);
// 重新编号支路的终止母线
index_i2e_tbus = slice(branch, [0], [T_BUS-1, T_BUS]) - 1;
i2e_tbus = get_multi(i2e, index_i2e_tbus)';
output_branch2 = assign(output_branch1, i2e_tbus, [0], [T_BUS-1, T_BUS]);
return (output_bus, output_gen, output_branch2, areas);
}
\ No newline at end of file
rspower/lib/isobservable.txt 0 → 100644
// ISOBSERVABLE Test for observability.
// returns 1 if the system is observable, 0 otherwise.
// 定义函数：判断系统是否可观测
fn isobservable(H, pv, pq) {
// 默认容差
tol = 1e-5;
// 是否检查不可观测的原因（0不检查，1检查。由于不支持字符串赋值自变量，故暂不支持检查。）
check_reason = 0;
// 获取H矩阵的行列数
m = size(H, 0);
n = size(H, 1);
// 计算H矩阵的秩
r = rank(H);
// 判断是否满秩（满秩则可观测）
if r < min(m, n) {
TorF = 0;
} else {
TorF = 1;
}
/*// 如果需要检查不可观测的原因且系统不可观测
if check_reason && !TorF {
// 存储零列索引和对应的变量名
idx_trivialColumns = [];
varNames = [];
// 遍历每一列检查是否为零列
for j in 0..n {
// 计算列的无穷范数
col_j = slice(H, [0], j);
normJ = norm_max(col_j);
if normJ < tol {
// 记录零列索引（转换为1基显示）
idx_trivialColumns = horzcat(idx_trivialColumns, [j + 1]);
// 获取变量名
varName = getVarName(j, pv, pq);
varNames = horzcat(varNames, [varName]);
}
}
// 如果存在零列（变量未被任何测量关联）
if ~~is_empty(idx_trivialColumns) {
println("Warning: The following variables are not observable since they are not related with any measurement!");
println("Variables:");
println(varNames);
println("Indices (1-based):");
println(idx_trivialColumns);
} else {
// 检查是否存在线性相关的列
for j in 0..n {
// 计算前j+1列的秩（因为RustScript是0基，j从0开始）
subH = slice(H, [0], [0, j+1]);
rr = rank(subH);
// 如果秩小于列数，说明存在线性相关
if rr < (j + 1) {
// 获取当前列
colJ = slice(H, [0], j);
varJName = getVarName(j, pv, pq);
// 检查当前列与前面所有列的线性相关性
for k in 0..j {
colK = slice(H, [0], k);
// 合并两列计算秩
combined = horzcat(colK, colJ);
if rank(combined) < 2 {
varKName = getVarName(k, pv, pq);
println("Warning: {}(th) column vector (w.r.t. {}) of H is linearly dependent of {}(th) column vector (w.r.t. {})!",
j+1, varJName, k+1, varKName);
return TorF;
}
}
}
}
println("Warning: No specific reason was found for system being not observable.");
}
}
*/
return TorF;
}
\ No newline at end of file
rspower/lib/make_jac.txt
...@@ -19,10 +19,7 @@ fn make_jac(baseMVA, bus, branch, gen, fullJac) { ...@@ -19,10 +19,7 @@ fn make_jac(baseMVA, bus, branch, gen, fullJac) {
V = set(V, index, gen_vg ./ abs(voltage_g) .* voltage_g); V = set(V, index, gen_vg ./ abs(voltage_g) .* voltage_g);
// build Jacobian // build Jacobian
ds_dv = dSbus_dV(Ybus, V, 0); [dSbus_dVa,dSbus_dVm] = dSbus_dV(Ybus, V, 0);
dSbus_dVa = slice(ds_dv, [0], [0,nb]);
dSbus_dVm = slice(ds_dv, [0], [nb,2*nb]);
// fulljac // fulljac
if fullJac { if fullJac {
...@@ -31,9 +28,8 @@ fn make_jac(baseMVA, bus, branch, gen, fullJac) { ...@@ -31,9 +28,8 @@ fn make_jac(baseMVA, bus, branch, gen, fullJac) {
j21 = imag(dSbus_dVa); j21 = imag(dSbus_dVa);
j22 = imag(dSbus_dVm); j22 = imag(dSbus_dVm);
} else { } else {
pv = bustypes_pv(bus, gen); [_,pv,pq] = bustypes(bus, gen);
pq = bustypes_pq(bus, gen); pv_pq_i = horzcat(pv', pq');
pv_pq_i = horzcat(pv, pq);
j11 = real(select(dSbus_dVa, pv_pq_i, pv_pq_i)); j11 = real(select(dSbus_dVa, pv_pq_i, pv_pq_i));
j12 = real(select(dSbus_dVm, pv_pq_i, pq)); j12 = real(select(dSbus_dVm, pv_pq_i, pq));
j21 = imag(select(dSbus_dVa, pq, pv_pq_i)); j21 = imag(select(dSbus_dVa, pq, pv_pq_i));
......
rspower/lib/newtonpf.txt
...@@ -67,9 +67,7 @@ fn newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt) { ...@@ -67,9 +67,7 @@ fn newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt) {
// update iteration counter // update iteration counter
i = i + 1; i = i + 1;
// evaluate Jacobian 因迭代中V变化，因此jac需要从dSbus_dV重新计算 // evaluate Jacobian 因迭代中V变化，因此jac需要从dSbus_dV重新计算
jac = dSbus_dV(Ybus, V, 0); [dSbus_dVa, dSbus_dVm] = dSbus_dV(Ybus, V, 0);
dSbus_dVa = slice(jac, [0], [0,nb]);
dSbus_dVm = slice(jac, [0], [nb,2*nb]);
j11 = real(select(dSbus_dVa, pv_pq_i, pv_pq_i)); j11 = real(select(dSbus_dVa, pv_pq_i, pv_pq_i));
j12 = real(select(dSbus_dVm, pv_pq_i, pq)); j12 = real(select(dSbus_dVm, pv_pq_i, pq));
......
rspower/lib/pfsoln.txt
...@@ -95,11 +95,11 @@ fn pfsoln(baseMVA, bus0, gen0, branch0, Ybus, Yf, Yt, V, ref, pv, pq, mpopt){ ...@@ -95,11 +95,11 @@ fn pfsoln(baseMVA, bus0, gen0, branch0, Ybus, Yf, Yt, V, ref, pv, pq, mpopt){
} // (terms are mult by 0 anyway) } // (terms are mult by 0 anyway)
// update Pg for slack gen(s) // update Pg for slack gen(s)
for k in 1..length(ref){ for k in 0..length(ref){
refgen = find(gbus == get_multi(get_multi(ref, [k-1]), 0) ); // which is(are) the reference gen(s)? refgen = find(gbus == get_multi(get_multi(ref, [k]), 0) + 1); // which is(are) the reference gen(s)?
PQ_refk = total_load(select(bus, get_multi(ref, [k-1] )), [], 1, [], [], mpopt); PQ_refk = total_load(select(bus, get_multi(ref, [k])), [], 1, [], [], mpopt);
Pd_refk = slice(PQ_refk, [0], [0,1]); Pd_refk = slice(PQ_refk, [0], [0,1]);
Pd_refk = [[c(0,0)]]; //Pd_refk = [[c(0,0)]];
gen_PG = set(slice(gen, [0], [PG-1, PG]), get_multi(on, [get_multi(refgen,0)]), real(get_multi(Sbus, [get_multi(refgen,0)])) * baseMVA + Pd_refk ); gen_PG = set(slice(gen, [0], [PG-1, PG]), get_multi(on, [get_multi(refgen,0)]), real(get_multi(Sbus, [get_multi(refgen,0)])) * baseMVA + Pd_refk );
gen = assign(gen, gen_PG, [0], [PG-1, PG] ); gen = assign(gen, gen_PG, [0], [PG-1, PG] );
...@@ -115,8 +115,8 @@ fn pfsoln(baseMVA, bus0, gen0, branch0, Ybus, Yf, Yt, V, ref, pv, pq, mpopt){ ...@@ -115,8 +115,8 @@ fn pfsoln(baseMVA, bus0, gen0, branch0, Ybus, Yf, Yt, V, ref, pv, pq, mpopt){
//----- update/compute branch power flows -----（计算支路潮流） //----- update/compute branch power flows -----（计算支路潮流）
out = find(slice(branch, [0], [BR_STATUS-1, BR_STATUS]) == 0); // out-of-service branches out = find(slice(branch, [0], [BR_STATUS-1, BR_STATUS]) == 0); // out-of-service branches
br = find(slice(branch, [0], [BR_STATUS-1, BR_STATUS] )); // in-service branches br = find(slice(branch, [0], [BR_STATUS-1, BR_STATUS] )); // in-service branches
Sf = get_multi(V, get_multi(br, slice(branch, [0], [F_BUS-1, F_BUS] ) ) - 1)' .* conj(select(Yf, br) * V) * baseMVA; // complex power at "from" bus(bus索引需-1) Sf = get_multi(V, get_multi(br, slice(branch, [0], [F_BUS-1, F_BUS]) - 1 ))' .* conj(select(Yf, br) * V) * baseMVA; // complex power at "from" bus(bus索引需-1)
St = get_multi(V, get_multi(br, slice(branch, [0], [T_BUS-1, T_BUS] ) ) - 1)' .* conj(select(Yt, br) * V) * baseMVA; // complex power injected at "to" bus St = get_multi(V, get_multi(br, slice(branch, [0], [T_BUS-1, T_BUS]) - 1 ))' .* conj(select(Yt, br) * V) * baseMVA; // complex power injected at "to" bus
branch_PF = zeros(size(branch,0), 1); branch_PF = zeros(size(branch,0), 1);
branch_PF = set(branch_PF, br, real(Sf)); branch_PF = set(branch_PF, br, real(Sf));
branch = horzcat(branch, branch_PF); branch = horzcat(branch, branch_PF);
...@@ -131,5 +131,4 @@ fn pfsoln(baseMVA, bus0, gen0, branch0, Ybus, Yf, Yt, V, ref, pv, pq, mpopt){ ...@@ -131,5 +131,4 @@ fn pfsoln(baseMVA, bus0, gen0, branch0, Ybus, Yf, Yt, V, ref, pv, pq, mpopt){
branch = horzcat(branch, branch_QT); branch = horzcat(branch, branch_QT);
return(bus, gen, branch); return(bus, gen, branch);
} }
\ No newline at end of file
rspower/lib/runse.txt 0 → 100644
//RUN_SE Run state estimation.
// [INPUT PARAMETERS]
// measure: measurements
// idx: measurement indices
// sigma: measurement variances
// [OUTPUT PARAMETERS]
// z: Measurement Vector. In the order of PF, PT, PG, Va, QF, QT, QG, Vm (if
// applicable), so it has ordered differently from original measurements
// z_est: Estimated Vector. In the order of PF, PT, PG, Va, QF, QT, QG, Vm
// (if applicable)
// error_sqrsum: Weighted sum of error squares
fn run_se(type_initialguess, V0) {
// 转换为内部母线编号
[se_i2e, se_bus, se_gen, se_branch] = ext2int(bus, gen, branch); //索引从“1”开始
// 获取母线类型
[ref, pv, pq] = bustypes(se_bus, se_gen); //索引从“0”开始
// 构建导纳矩阵
[Ybus, Yf, Yt] = make_y_bus(baseMVA, se_bus, se_branch);
Ybus = full(Ybus); // 转为稠密矩阵
Yf = full(Yf);
Yt = full(Yt);
// 初始电压猜测
V0 = getV0(se_bus, se_gen, type_initialguess, V0);
// 执行状态估计
[V, success, iterNum, z, z_est, error_sqrsum] = doSE(baseMVA, se_bus, se_gen, se_branch, Ybus, Yf, Yt, V0, ref, pv, pq);
// 更新发电机有功功率（PG列）
idx_zPG_len = length(idx_zPG);
if idx_zPG_len > 0 {
cnt = length(idx_zPF) + length(idx_zPT);
// 提取z_est中对应PG的部分并转换为标幺值
z_est_PG = slice(z_est, [cnt, cnt + idx_zPG_len], [0]);
gen_PG = z_est_PG * baseMVA;
// 更新gen矩阵的PG列
se_gen_PG = set(slice(se_gen, [0], [PG - 1, PG]), idx_zPG - 1, gen_PG);
se_gen = assign(se_gen, se_gen_PG, [0], [PG - 1, PG]);
}
// 更新发电机无功功率（QG列）
idx_zQG_len = length(idx_zQG);
if idx_zQG_len > 0 {
cnt = length(idx_zPF) + length(idx_zPT) + length(idx_zPG) + length(idx_zVa) + length(idx_zQF) + length(idx_zQT);
// 提取z_est中对应QG的部分并转换为标幺值
z_est_QG = slice(z_est, [cnt, cnt + idx_zQG_len], [0]);
gen_QG = z_est_QG * baseMVA;
// 更新gen矩阵的QG列
se_gen_QG = set(slice(se_gen, [0], [QG - 1, QG]), idx_zQG - 1, gen_QG);
se_gen = assign(se_gen, se_gen_QG, [0], [QG - 1, QG]);
}
// 更新潮流解
[new_bus, new_gen, new_branch] = pfsoln(baseMVA, se_bus, se_gen, se_branch, Ybus, Yf, Yt, V, ref, pv, pq, []);
// 转换回原始母线编号
[output_bus, output_gen, output_branch] = int2ext(se_i2e, new_bus, new_gen, new_branch, areas);
// 输出结果
//outputpfsoln(baseMVA, bus, gen, branch, success, et, 1, iterNum);
//outputsesoln(idx, sigma, z, z_est, error_sqrsum);
return (baseMVA, output_bus, output_gen, output_branch, success, z, z_est, error_sqrsum);
}
\ No newline at end of file
rustscript/RustScript手册.md
# RustScript 语言规范 # RustScript 语言规范
...@@ -1251,6 +1251,16 @@ for i in 0..b { ...@@ -1251,6 +1251,16 @@ for i in 0..b {
``` ```
**注意**：循环序列需“0..10”格式，不能跟tensor。 **注意**：循环序列需“0..10”格式，不能跟tensor。
```rustscript
//while语句
max_num = 10;
a = 0;
converged = 0;
while (!converged && a < max_num) {
a = a + 1;
}
```
## 9. 函数定义 ## 9. 函数定义
### 9.1 基本函数语法 ### 9.1 基本函数语法
......
Markdown 格式
0% 或
您添加了 0 人 到此讨论。请谨慎行事。
请先完成此评论的编辑！
请 注册 或者 后发表评论