refactor:make powerflow work

rspower/lib/newtonpf.txt

@@ -48,13 +48,13 @@ fn newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt) {
     nb = length(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
-    pv_pq_i = vertcat(pv, pq);
+    j1 = 0;         j2 = npv;           // j1:j2 - V angle of pv buses
+    j3 = j2;    j4 = j2 + npq;      // j3:j4 - V angle of pq buses
+    j5 = j4;    j6 = j4 + npq;      // j5:j6 - V mag of pq buses
+    pv_pq_i = horzcat(pv, pq);
     // evaluate F(x0)
     mis = V .* conj(Ybus * V) - Sbus;
-    F = vertcat(real(get_multi(mis, pv_pq_i)), imag(get_mult(mis, pq)));
+    F = vertcat(real(get_multi(mis, pv_pq_i))', imag(get_multi(mis, pq))');
     // check tolerance
     normF = norm_max(F);
     // do Newton iterations
@@ -76,7 +76,7 @@ fn newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt) {
         );
         // compute update step
         // dx = linsolve(J, -F);
-        dx = dSbus_dVm \ -F;
+        dx = J \ -F;
         // update voltage vector
         if npv {
             Va = set2(Va, pv, slice(dx, [j1,j2], [0]));
@@ -90,9 +90,9 @@ fn newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt) {
         Va = angle(V);          // we wrapped around with a negative Vm
 
         mis = V .* conj(Ybus * V) - Sbus;
-        F = vertcat(real(get_multi(mis, pv_pq_i)), imag(get_mult(mis, pq)));
+        F = vertcat(real(get_multi(mis, pv_pq_i))', imag(get_multi(mis, pq))');
         // check for convergence
-         normF = norm_max(F);
+        normF = norm_max(F);
         if normF < tol {
             converged = 1;
         }

rspower/lib/runpf.txt

@@ -5,7 +5,8 @@ fn runpf() {
     // bus_gen_status(gen(:, GEN_BUS)) = gen(:, GEN_STATUS) > 0;
     nb = size(bus, 0);
     ng = size(gen, 0);
-    Cg = sparse(slice(gen, [0], [GEN_BUS-1,GEN_BUS]), range(0,ng-1), slice(gen, [0], [GEN_STATUS-1,GEN_STATUS]) > 0, nb, ng);  // gen connection matrix
+    g_i = slice(gen, [0], [GEN_BUS-1,GEN_BUS]) - 1; // generator bus indices
+    Cg = sparse(g_i, range(0,ng-1), 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