小修改

rspower/examples/test_newtonpf.txt

@@ -8,6 +8,31 @@
 #include ../lib/make_sdzip.txt
 #include ../lib/make_sdzip.txt
 #include ../lib/make_sbus.txt
 #include ../lib/make_sbus.txt
 
 
+// get bus index lists of each type of bus
+// get generator status
+// bus_gen_status = zeros(size(bus, 1), 1);
+// bus_gen_status(gen(:, GEN_BUS)) = gen(:, GEN_STATUS) > 0;
+nb = size(bus, 1);
+ng = size(gen, 1);
+Cg = sparse(slice(gen, [0], [GEN_BUS-1,GEN_BUS]), (1: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
+
+
+// initial state
+// V0    = ones(size(bus, 1), 1);            // flat start
+V0  = slice(bus, [0], [VM-1,VM]) .* exp(c(0,1) * pi/180 * slice(bus, [0], [VA-1,VA]));
+vcb = ones(size(V0));           // create mask of voltage-controlled buses
+vcb(pq) = 0;                    // exclude PQ buses
+k = find(vcb(gbus));            // in-service gens at v-c buses
+V0(gbus(k)) = gen(on(k), VG) ./ abs(V0(gbus(k))).* V0(gbus(k));
+
 Ybus = make_y_bus(baseMVA, bus, branch);
 Ybus = make_y_bus(baseMVA, bus, branch);
 r = make_sbus(baseMVA, bus, gen, 1, 1, 1);
 r = make_sbus(baseMVA, bus, gen, 1, 1, 1);
 return r;
 return r;
\ No newline at end of file

rspower/lib/newtonpf.txt

@@ -45,4 +45,14 @@ fn newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt) {
     V = V0;
     V = V0;
     Va = angle(V);
     Va = angle(V);
     Vm = abs(V);
     Vm = abs(V);
+    
+    npv = length(pv);
+    npq = length(pq);
+    j1 = 1;         j2 = npv;           // j1:j2 - V angle of pv buses
+    j3 = j2 + 1;    j4 = j2 + npq;      // j3:j4 - V angle of pq buses
+    j5 = j4 + 1;    j6 = j4 + npq;      // j5:j6 - V mag of pq buses
+
+    // evaluate F(x0)
+    mis = V .* conj(Ybus * V) - Sbus(Vm);
+   ·
 }
 }
\ No newline at end of file