pfsoln函数开发以及make_ybus函数修改

rspower/examples/test_pfsoln.txt

+#include ../data/case14.txt
+#include ../lib/idx_gen.txt
+#include ../lib/idx_bus.txt
+#include ../lib/idx_brch.txt
+#include ../lib/make_y_bus.txt
+#include ../lib/make_y_bus_Yf.txt
+#include ../lib/make_y_bus_Yt.txt
+#include ../lib/dsbus_dv.txt
+#include ../lib/make_sdzip.txt
+#include ../lib/make_sbus.txt
+#include ../lib/newtonpf.txt
+#include ../lib/runpf.txt
+#include ../lib/pfsoln.txt
+
+
+info("Running power flow on case14...");
+for i in 0..5 {
+    r = runpf();
+}
+info("Run power flow on case14 end");
+
+Ybus = make_y_bus(baseMVA, bus, branch);
+Yf = make_y_bus_Yf(baseMVA, bus, branch);
+Yt = make_y_bus_Yt(baseMVA, bus, branch);
+V = r;
+bus_type = slice(bus, [0], [BUS_TYPE-1,BUS_TYPE]);
+ref = find(bus_type == REF) + 1;
+pv  = find(bus_type == PV) + 1;
+pq  = find(bus_type == PQ) + 1;
+mpopt = 12;
+bus_output = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt);
+
+return bus_output;
\ No newline at end of file

rspower/lib/make_y_bus.txt

@@ -39,8 +39,10 @@ fn make_y_bus(baseMVA, bus, branch) {
         // small case
         // build Yf and Yt such that Yf * V is the vector of complex branch currents injected
         // at each branch's "from" bus, and Yt is the same for the "to" bus end
-        i = horzcat(range(0, nl),range(0, nl)) - 1;                           // double set of row indices
+        i = horzcat(range(0, nl),range(0, nl));                           // double set of row indices
         j = range(0,nb);
+        upper = horzcat(Yff', Yft');
+        lower = horzcat(Ytf', Ytt');
         Yf = full(sparse(i, horzcat(f, t), upper, nl, nb));
         Yt = full(sparse(i, horzcat(f, t), lower, nl, nb));
 

rspower/lib/make_y_bus_Yf.txt

+fn make_y_bus_Yf(baseMVA, bus, branch) {
+    // constants
+    nb = size(bus, 0);          // number of buses
+    nl = size(branch, 0);       // number of lines
+
+    // for each branch, compute the elements of the branch admittance matrix where
+    //
+    //      | If |   | Yff  Yft |   | Vf |
+    //      |    | = |          | * |    |
+    //      | It |   | Ytf  Ytt |   | Vt |
+    //
+    stat = slice(branch, [0], [BR_STATUS-1,BR_STATUS]);     // ones at in-service branches
+    Ys = stat ./ (slice(branch, [0], [BR_R-1,BR_R])
+            + c(0,1) * slice(branch, [0], [BR_X-1,BR_X]));  // series admittance
+    Bc = stat .* slice(branch, [0], [BR_B-1,BR_B]);         // line charging susceptance
+    tap_col = slice(branch, [0], [TAP-1,TAP]); // tap ratio column
+
+    index = find(tap_col);      // indices of non-zero tap ratios
+    tap_init = set(ones(nl, 1), index, get_multi(tap_col, index)); // assign non-zero tap ratios
+    tap = tap_init .* exp(c(0,1)*pi/180 * slice(branch, [0], [SHIFT-1,SHIFT])); // add phase shifters
+    Ytt = Ys + c(0,1) * Bc/2;
+    Yff = Ytt ./ (tap .* conj(tap));
+    Yft = - Ys ./ conj(tap);
+    Ytf = - Ys ./ tap;
+
+    // compute shunt admittance
+    // if Psh is the real power consumed by the shunt at V = 1.0 p.u.
+    // and Qsh is the reactive power injected by the shunt at V = 1.0 p.u.
+    // then Psh - j Qsh = V * conj(Ysh * V) = conj(Ysh) = Gs - j Bs,
+    // i.e. Ysh = Psh + j Qsh, so ...
+
+    // vector of shunt admittances
+    Ysh = (slice(bus, [0], [GS-1,GS]) + c(0,1) * slice(bus, [0], [BS-1,BS])) / baseMVA;
+
+    // bus indices
+    f = slice(branch, [0], F_BUS-1) - 1; // list of "from" buses
+    t = slice(branch, [0], T_BUS-1) - 1; // list of "to" buses
+    if nb < 300 {
+        // small case
+        // build Yf and Yt such that Yf * V is the vector of complex branch currents injected
+        // at each branch's "from" bus, and Yt is the same for the "to" bus end
+        i = horzcat(range(0, nl),range(0, nl));                           // double set of row indices
+        j = range(0,nb);
+        upper = horzcat(Yff', Yft');
+        lower = horzcat(Ytf', Ytt');
+        Yf = full(sparse(i, horzcat(f, t), upper, nl, nb));
+        Yt = full(sparse(i, horzcat(f, t), lower, nl, nb));
+
+        // build Ybus
+        // branch admittances + shunt admittance
+        Ybus = full(sparse(horzcat(f,f,t,t), horzcat(f,t,f,t), vertcat(Yff,Yft,Ytf,Ytt), nb, nb)) +
+               full(sparse(j, j, Ysh, nb, nb));
+    } else {
+        // large case running on MATLAB
+        // build connection matrices
+        i = range(0, nl);
+        j = range(0, nb);
+        Cf = full(sparse(i, f, ones(nl, 1), nl, nb));      // connection matrix for line & from buses
+        Ct = full(sparse(i, t, ones(nl, 1), nl, nb));      // connection matrix for line & to buses
+
+        // build Yf and Yt such that Yf * V is the vector of complex branch currents injected
+        // at each branch's "from" bus, and Yt is the same for the "to" bus end
+
+        Yf = full(sparse(i, i, Yff, nl, nl)) * Cf + full(sparse(i, i, Yft, nl, nl)) * Ct;
+        Yt = full(sparse(i, i, Ytf, nl, nl)) * Cf + full(sparse(i, i, Ytt, nl, nl)) * Ct;
+
+        // build Ybus
+        // branch admittances + shunt admittance
+        Ybus = Cf' * Yf + Ct' * Yt +
+               full(sparse(j, j, Ysh, nb, nb));
+    }
+
+    return Yf;
+}
\ No newline at end of file

rspower/lib/make_y_bus_Yt.txt

+fn make_y_bus_Yt(baseMVA, bus, branch) {
+    // constants
+    nb = size(bus, 0);          // number of buses
+    nl = size(branch, 0);       // number of lines
+
+    // for each branch, compute the elements of the branch admittance matrix where
+    //
+    //      | If |   | Yff  Yft |   | Vf |
+    //      |    | = |          | * |    |
+    //      | It |   | Ytf  Ytt |   | Vt |
+    //
+    stat = slice(branch, [0], [BR_STATUS-1,BR_STATUS]);     // ones at in-service branches
+    Ys = stat ./ (slice(branch, [0], [BR_R-1,BR_R])
+            + c(0,1) * slice(branch, [0], [BR_X-1,BR_X]));  // series admittance
+    Bc = stat .* slice(branch, [0], [BR_B-1,BR_B]);         // line charging susceptance
+    tap_col = slice(branch, [0], [TAP-1,TAP]); // tap ratio column
+
+    index = find(tap_col);      // indices of non-zero tap ratios
+    tap_init = set(ones(nl, 1), index, get_multi(tap_col, index)); // assign non-zero tap ratios
+    tap = tap_init .* exp(c(0,1)*pi/180 * slice(branch, [0], [SHIFT-1,SHIFT])); // add phase shifters
+    Ytt = Ys + c(0,1) * Bc/2;
+    Yff = Ytt ./ (tap .* conj(tap));
+    Yft = - Ys ./ conj(tap);
+    Ytf = - Ys ./ tap;
+
+    // compute shunt admittance
+    // if Psh is the real power consumed by the shunt at V = 1.0 p.u.
+    // and Qsh is the reactive power injected by the shunt at V = 1.0 p.u.
+    // then Psh - j Qsh = V * conj(Ysh * V) = conj(Ysh) = Gs - j Bs,
+    // i.e. Ysh = Psh + j Qsh, so ...
+
+    // vector of shunt admittances
+    Ysh = (slice(bus, [0], [GS-1,GS]) + c(0,1) * slice(bus, [0], [BS-1,BS])) / baseMVA;
+
+    // bus indices
+    f = slice(branch, [0], F_BUS-1) - 1; // list of "from" buses
+    t = slice(branch, [0], T_BUS-1) - 1; // list of "to" buses
+    if nb < 300 {
+        // small case
+        // build Yf and Yt such that Yf * V is the vector of complex branch currents injected
+        // at each branch's "from" bus, and Yt is the same for the "to" bus end
+        i = horzcat(range(0, nl),range(0, nl));                           // double set of row indices
+        j = range(0,nb);
+        upper = horzcat(Yff', Yft');
+        lower = horzcat(Ytf', Ytt');
+        Yf = full(sparse(i, horzcat(f, t), upper, nl, nb));
+        Yt = full(sparse(i, horzcat(f, t), lower, nl, nb));
+
+        // build Ybus
+        // branch admittances + shunt admittance
+        Ybus = full(sparse(horzcat(f,f,t,t), horzcat(f,t,f,t), vertcat(Yff,Yft,Ytf,Ytt), nb, nb)) +
+               full(sparse(j, j, Ysh, nb, nb));
+    } else {
+        // large case running on MATLAB
+        // build connection matrices
+        i = range(0, nl);
+        j = range(0, nb);
+        Cf = full(sparse(i, f, ones(nl, 1), nl, nb));      // connection matrix for line & from buses
+        Ct = full(sparse(i, t, ones(nl, 1), nl, nb));      // connection matrix for line & to buses
+
+        // build Yf and Yt such that Yf * V is the vector of complex branch currents injected
+        // at each branch's "from" bus, and Yt is the same for the "to" bus end
+
+        Yf = full(sparse(i, i, Yff, nl, nl)) * Cf + full(sparse(i, i, Yft, nl, nl)) * Ct;
+        Yt = full(sparse(i, i, Ytf, nl, nl)) * Cf + full(sparse(i, i, Ytt, nl, nl)) * Ct;
+
+        // build Ybus
+        // branch admittances + shunt admittance
+        Ybus = Cf' * Yf + Ct' * Yt +
+               full(sparse(j, j, Ysh, nb, nb));
+    }
+
+    return Yt;
+}
\ No newline at end of file

rspower/lib/pfsoln.txt

+//PFSOLN  Updates bus, gen, branch data structures to match power flow soln.
+//   [BUS, GEN, BRANCH] = PFSOLN(BASEMVA, BUS0, GEN0, BRANCH0, ...
+//                                   YBUS, YF, YT, V, REF, PV, PQ, MPOPT)
+//  调用本函数时，需要调用如下文件
+//  #include ../lib/idx_gen.txt
+//  #include ../lib/idx_bus.txt
+//  #include ../lib/idx_brch.txt
+//  #include total_load.txt
+
+fn pfsoln(baseMVA, bus0, gen0, branch0, Ybus, Yf, Yt, V, ref, pv, pq, mpopt){
+
+    // initialize return values（创建原始数据的副本）
+    bus     = bus0;
+    gen     = gen0;
+    branch  = branch0;
+
+    //----- update bus voltages -----(将求解得到的复数电压转换为幅值和角度)
+    bus = assign(bus, abs(V), [0], [VM-1, VM] );
+    bus = assign(bus, angle(V) * 180 / pi, [0], [VA-1, VA] );
+
+    //----- update Qg for gens at PV/slack buses and Pg for slack bus(es) -----
+    // generator info
+    //取非PQ且GEN_STATUS=1的generators
+    on_temp = find(slice(gen, [0], [GEN_STATUS-1, GEN_STATUS]) > 0);    //which generators are on?
+    index_gbus_type = get_multi(slice(gen, [0], [GEN_BUS-1, GEN_BUS]), on_temp) - 1;    //得到GEN_STATUS=1的发电机母线节点索引
+    on = find(get_multi(slice(bus, [0], [BUS_TYPE-1, BUS_TYPE]), index_gbus_type) != PQ)'; // which generators are not at PQ buses
+    off = find(slice(gen, [0], [GEN_STATUS-1, GEN_STATUS]) <= 0)';       // which generators are off?
+    gbus = get_multi(slice(gen, [0], [GEN_BUS-1, GEN_BUS]), on)';         // what buses are they at?(注意索引需-1)
+
+    // compute total injected bus powers(计算节点注入功率)
+    Sbus = get_multi(V, gbus-1)' .* conj(select(Ybus, gbus-1) * V );
+
+    // update Qg for generators at PV/slack buses
+    gen_QG = set(slice(gen, [0], [QG-1, QG]), off, zeros(length(off), 1));  // zero out off-line Qg(离线发电机QG=0)
+    gen = assign(gen, gen_QG, [0], [QG-1, QG] );
+    // don't touch the ones at PQ buses
+    [Pd_gbus, Qd_gbus] = total_load(bus(gbus, :), [], 'bus', [], mpopt);    //total_load()函数现在没有，要另外开发(gbus索引注意是否-1)
+    gen_QG = set(slice(gen, [0], [QG-1, QG]), on, imag(Sbus) * baseMVA + Qd_gbus);  // inj Q + local Qd
+    gen = assign(gen, gen_QG, [0], [QG-1, QG] );
+
+    // ... at this point any buses with more than one generator will have
+    // the total Q dispatch for the bus assigned to each generator. This
+    // must be split between them. We do it first equally, then in proportion
+    // to the reactive range of the generator.（处理多发电机母线的无功分配）
+
+    if length(on) > 1 {
+        // build connection matrix, element i, j is 1 if gen on(i) at bus j is ON
+        nb = size(bus, 0);
+        ngon = size(on, 0);
+        Cg = full(sparse(range(0,ngon)', (gbus-1)', ones(ngon, 1), ngon, nb));
+
+        // divide Qg by number of generators at the bus to distribute equally(平均分配)
+        ngg = Cg * sum(Cg)';    // ngon x 1, number of gens at this gen's bus
+        gen_QG = set(slice(gen, [0], [QG-1, QG]), on, get_multi(slice(gen, [0], [QG-1, QG]), on)' ./ ngg);
+        gen = assign(gen, gen_QG, [0], [QG-1, QG] );
+
+        // set finite proxy M for infinite limits (for ~ proportional splitting)(处理无限限值)
+        // equal to sum over all gens at bus of abs(Qg) plus any finite Q limits
+        Qmin = select(gen, on, QMIN-1);
+        Qmax = select(gen, on, QMAX-1);
+        M = abs(select(gen, on, QG-1));
+        index_Qmax_inf = find( ~~is_inf(Qmax) );
+        M = set(M, index_Qmax_inf, get_multi(M, index_Qmax_inf) + abs(get_multi(Qmax, index_Qmax_inf)) );
+        index_Qmin_inf = find( ~~is_inf(Qmin) );
+        M = set(M, index_Qmin_inf, get_multi(M, index_Qmin_inf) + abs(get_multi(Qmin, index_Qmin_inf)) );
+        M = Cg * Cg' * M;   // each gen gets sum over all gens at same bus
+        // replace +/- Inf limits with proxy +/- M
+        Qmin = set(Qmin, find(Qmin == INF), get_multi(M, find(Qmin == INF) ) );
+        Qmin = set(Qmin, find(Qmin == -INF), -get_multi(M, find(Qmin == -INF) ) );
+        Qmax = set(Qmax, find(Qmax == INF), get_multi(M, find(Qmax == INF) ) );
+        Qmax = set(Qmax, find(Qmax == -INF), -get_multi(M, find(Qmax == -INF) ) );
+
+        // divide proportionally（按比例分配）
+        Cmin = full(sparse(range(0,ngon)', (gbus-1)', Qmin, ngon, nb) );
+        Cmax = full(sparse(range(0,ngon)', (gbus-1)', Qmax, ngon, nb) );
+        Qg_tot = Cg' * select(gen, on, QG-1);     // nb x 1 vector of total Qg at each bus(母线总无功)
+        Qg_min = sum(Cmin)';            // nb x 1 vector of min total Qg at each bus
+        Qg_max = sum(Cmax)';            // nb x 1 vector of max total Qg at each bus
+        gen_QG = set(slice(gen, [0], [QG-1, QG]), on, Qmin + (Cg * ((Qg_tot - Qg_min)./(Qg_max - Qg_min+eps))) .* (Qmax - Qmin) );
+        gen = assign(gen, gen_QG, [0], [QG-1, QG] );
+
+        // fix gens at buses with Qg range = 0 (use equal violation for all)（处理零范围特殊情况）
+        ig = find(abs(Cg * (Qg_min - Qg_max)) < 10*eps);  // gens at buses with Qg range = 0
+        if ~~(length(ig)==0){   //matlab中是isempty()函数，这里用length==0代替
+            ib = find(sum(select(Cg,ig), 0)')';   // buses with Qg range = 0
+            // total mismatch @ bus div by number of gens
+            mis_value = (get_multi(Qg_tot, ib)' - get_multi(Qg_min, ib)' ) ./ sum(select(Cg, [], ib)' , 1)';  
+            mis = full(sparse(range(0,nb)', zeros(nb,1), set(zeros(nb,1), ib, mis_value), nb, 1) );
+            gen_QG = set(slice(gen, [0], [QG-1, QG]), get_multi(on, ig), get_multi(Qmin, ig)' + select(Cg, ig) * mis);
+            gen = assign(gen, gen_QG, [0], [QG-1, QG] );
+        }
+    }                                            // (terms are mult by 0 anyway)
+
+    // update Pg for slack gen(s)
+    for k in 1..length(ref){
+        refgen = find(gbus == get_multi(get_multi(ref, [k-1]), 0) );              // which is(are) the reference gen(s)?
+        Pd_refk = total_load(bus(ref(k-1), :), [], 'bus', [], mpopt);     //total_load()函数现在没有，要另外开发
+        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] );
+
+        if length(refgen) > 1{          // more than one generator at this ref bus
+            // subtract off what is generated by other gens at this bus
+            refgen_PG_value = get_multi(slice(gen, [0], [PG-1, PG]), get_multi(on, [get_multi(refgen,0)]) ) 
+                                - sum(get_multi(slice(gen, [0], [PG-1, PG]), get_multi(on, get_multi(refgen,range(1,length(refgen)))) ));
+            gen_PG = set(slice(gen, [0], [PG-1, PG]), get_multi(on, [get_multi(refgen,0)]),  refgen_PG_value);
+            gen = assign(gen, gen_PG, [0], [PG-1, PG] );
+        }
+    }
+
+    //----- update/compute branch power flows -----（计算支路潮流）
+    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
+    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
+    branch_PF = zeros(size(branch,0), 1);
+    branch_PF = set(branch_PF, br, real(Sf));
+    branch = horzcat(branch, branch_PF);
+    branch_QF = zeros(size(branch,0), 1);
+    branch_QF = set(branch_QF, br, imag(Sf));
+    branch = horzcat(branch, branch_QF);
+    branch_PT = zeros(size(branch,0), 1);
+    branch_PT = set(branch_PT, br, real(St));
+    branch = horzcat(branch, branch_PT);
+    branch_QT = zeros(size(branch,0), 1);
+    branch_QT = set(branch_QT, br, imag(St));
+    branch = horzcat(branch, branch_QT);
+
+    //return bus, gen, branch;
+    return bus;
+}
\ No newline at end of file

rspower/lib/total_load.txt

+//TOTAL_LOAD Returns vector of total load in each load zone.
+//   PD = TOTAL_LOAD(MPC)
+//   PD = TOTAL_LOAD(MPC, LOAD_ZONE)
+//   PD = TOTAL_LOAD(MPC, LOAD_ZONE, OPT)
+//   PD = TOTAL_LOAD(MPC, LOAD_ZONE, OPT, MPOPT)
+//   PD = TOTAL_LOAD(BUS)
+//   PD = TOTAL_LOAD(BUS, GEN)
+//   PD = TOTAL_LOAD(BUS, GEN, LOAD_ZONE)
+//   PD = TOTAL_LOAD(BUS, GEN, LOAD_ZONE, OPT)
+//   PD = TOTAL_LOAD(BUS, GEN, LOAD_ZONE, OPT, MPOPT)
+//   [PD, QD] = TOTAL_LOAD(...) returns both active and reative power
+//   demand for each zone.
+//
+//   MPC - standard MATPOWER case struct
+//
+//   BUS - standard BUS matrix with nb rows, where the fixed active
+//       and reactive loads are specified in columns PD and QD
+//
+//   GEN - (optional) standard GEN matrix with ng rows, where the
+//       dispatchable loads are specified by columns PG, QG, PMIN,
+//       QMIN and QMAX (in rows for which ISLOAD(GEN) returns true).
+//       If GEN is empty, it assumes there are no dispatchable loads.
+//
+//   LOAD_ZONE - (optional) nb element vector where the value of
+//       each element is either zero or the index of the load zone
+//       to which the corresponding bus belongs. If LOAD_ZONE(b) = k
+//       then the loads at bus b will added to the values of PD(k) and
+//       QD(k). If LOAD_ZONE is empty, the default is defined as the areas
+//       specified in the BUS matrix, i.e. LOAD_ZONE = BUS(:, BUS_AREA)
+//       and load will have dimension = MAX(BUS(:, BUS_AREA)). LOAD_ZONE
+//       can also take the following string values:
+//           'all'  - use a single zone for the entire system (return scalar)
+//           'area' - use LOAD_ZONE = BUS(:, BUS_AREA), same as default
+//           'bus'  - use a different zone for each bus (i.e. to compute
+//               final values of bus-wise loads, including voltage dependent
+//               fixed loads and or dispatchable loads)
+//
+//   OPT - (optional) option struct, with the following fields:
+//           'type'  -  string specifying types of loads to include, default
+//                      is 'BOTH' if GEN is provided, otherwise 'FIXED'
+//               'FIXED'        : sum only fixed loads
+//               'DISPATCHABLE' : sum only dispatchable loads
+//               'BOTH'         : sum both fixed and dispatchable loads
+//           'nominal' -  1 : use nominal load for dispatchable loads
+//                        0 : (default) use actual realized load for
+//                             dispatchable loads
+//
+//       For backward compatibility with MATPOWER 4.x, OPT can also
+//       take the form of a string, with the same options as OPT.type above.
+//       In this case, again for backward compatibility, it is the "nominal"
+//       load that is computed for dispatchable loads, not the actual
+//       realized load. Using a string for OPT is deprecated and
+//       will be removed in a future version.
+//
+//   MPOPT - (optional) MATPOWER options struct, which may specify
+//       a voltage dependent (ZIP) load model for fixed loads
+//
+//   Examples:
+//       Return the total active load for each area as defined in BUS_AREA.
+//
+//       Pd = total_load(bus);
+//
+//       Return total active and reactive load, fixed and dispatchable, for
+//       entire system.
+//
+//       [Pd, Qd] = total_load(bus, gen, 'all');
+//
+//       Return the total of the nominal dispatchable loads at buses 10-20.
+//
+//       load_zone = zeros(nb, 1);
+//       load_zone(10:20) = 1;
+//       opt = struct('type', 'DISPATCHABLE', 'nominal', 1);
+//       Pd = total_load(mpc, load_zone, opt)
+//
+//   See also SCALE_LOAD.
+
+fn  total_load(bus, gen, load_zone, opt, mpopt){
+
+    //return Pd, Qd;
+    return Pd;
+}
\ No newline at end of file