refactor:change pf alg

rspower/lib/dsbus_dv.txt

@@ -102,5 +102,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 dSbus_dV1+dSbus_dV2;
+    return horzcat(dSbus_dV1, dSbus_dV2);
 }
 }
\ No newline at end of file

rspower/lib/make_jac.txt

@@ -16,6 +16,6 @@ fn make_jac(baseMVA, bus, branch, gen) {
     on_k = get_multi(on, k);
     on_k = get_multi(on, k);
     gen_vg = get_multi(slice(gen, [0], [VG-1,VG]), on_k);
     gen_vg = get_multi(slice(gen, [0], [VG-1,VG]), on_k);
     V = set(V, index, gen_vg ./ abs(voltage_g) .* voltage_g);
     V = set(V, index, gen_vg ./ abs(voltage_g) .* voltage_g);
-    dsbus_dv = dSbus_dV(Ybus, V, 0.);
+    dsbus_dv = dSbus_dV(Ybus, V, 0);
     return dsbus_dv;
     return dsbus_dv;
 }
 }
\ No newline at end of file

rspower/lib/newtonpf.txt

@@ -45,7 +45,7 @@ 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);
-    
+    nb = length(V);
     npv = length(pv);
     npv = length(pv);
     npq = length(pq);
     npq = length(pq);
     j1 = 1;         j2 = npv;           // j1:j2 - V angle of pv buses
     j1 = 1;         j2 = npv;           // j1:j2 - V angle of pv buses
@@ -61,10 +61,12 @@ fn newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt) {
         // update iteration counter
         // update iteration counter
         i = i + 1;
         i = i + 1;
         // evaluate Jacobian
         // evaluate Jacobian
-        J = dSbus_dV(Ybus, V);
+        J = dSbus_dV(Ybus, V, 0);
+        dSbus_dVa = slice(J, [0], [0,nb]);
+        dSbus_dVm = slice(J, [0], [nb,2*nb]);
         // compute update step
         // compute update step
         // dx = linsolve(J, -F);
         // dx = linsolve(J, -F);
-        dx = J \ -F;
+        dx = dSbus_dVm \ -F;
         // update voltage vector
         // update voltage vector
         V = V + dx;
         V = V + dx;