refactor: change pf alg

rspower/lib/make_sbus.txt

@@ -28,12 +28,21 @@
 
 // define named indices into bus, gen matrices
 
-fn make_sbus(baseMVA, bus, gen, mpopt, Vm, Sg) {
+fn make_sbus(baseMVA, bus, gen) {
     nb = size(bus, 0);
     Sd = make_sdzip(baseMVA, bus, mpopt);
     on = find(slice(gen, [0], [GEN_STATUS-1, GEN_STATUS]) > 0);  // which generators are on?
     gbus = get_multi(slice(gen, [0], [GEN_BUS-1, GEN_BUS]), on) - 1; // what buses are they at?
     ngon = size(on, 0);
     Cg = sparse(gbus, range(0, ngon-1), 1, nb, ngon);
-    return Cg;
+    Gp = slice(gen, [0], [PG-1, PG]);
+    Gq = slice(gen, [0], [PG-1, PG]);
+    Sbusg = Cg * (get_multi(Gp, on)  + c(0,1) * get_multi(Gq, on)) / baseMVA;
+    z = slice(Sd, [0], [0,1]);
+    i = slice(Sd, [0], [1,2]);
+    p = slice(Sd, [0], [2,3]);
+    Vm = ones(nb, 1);
+    Sbusd = p + i .* Vm + z .* Vm.^2;
+    // complex power injection at each bus
+    return Sbusg' - Sbusd;
 }
\ No newline at end of file

rspower/lib/make_sdzip.txt

@@ -25,5 +25,5 @@ fn make_sdzip(baseMVA, bus, mpopt) {
     z = (pd * get(pw,2)  + c(0,1) * qd * get(qw,2)) / baseMVA;
     i = (pd * get(pw,1)  + c(0,1) * qd * get(qw,1)) / baseMVA;
     p = (pd * get(pw,0)  + c(0,1) * qd * get(qw,0)) / baseMVA;
-    return p;
+    return horzcat(z,i,p);
 }
\ No newline at end of file

rspower/lib/newtonpf.txt

@@ -53,8 +53,7 @@ fn newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt) {
     j5 = j4 + 1;    j6 = j4 + npq;      // j5:j6 - V mag of pq buses
 
     // evaluate F(x0)
-    mis = V .* conj(Ybus * V);
-    F = mis;
+    F = V .* conj(Ybus * V) - Sbus;
     // check tolerance
     normF = norm_max(F);
     // do Newton iterations
@@ -64,10 +63,17 @@ fn newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt) {
         // evaluate Jacobian
         J = dSbus_dV(Ybus, V);
         // compute update step
-        dx = linsolve(J, -F);
+        // dx = linsolve(J, -F);
+        dx = J \ -F;
+        // update voltage vector
+        V = V + dx';
+
+        F = V .* conj(Ybus * V) - Sbus;
+        // check for convergence
+         normF = norm_max(F);
         if normF < tol {
             converged = 1;
         }
     }
-    return dx;
+    return V;
 }
\ No newline at end of file

rspower/lib/runpf.txt

@@ -36,7 +36,7 @@ fn runpf() {
     qlim = 0;
     loop {
         // compute bus power injections
-        Sbus = make_sbus(baseMVA, bus, gen, 1, 1, 1);
+        Sbus = make_sbus(baseMVA, bus, gen);
         V = newtonpf(Ybus, Sbus, V_init, ref, pv, pq, [1e-6, 100]);
         if ~~is_empty(V) && qlim {
             // if V is empty, it means the power flow did not converge