refactor: add dsbus_dv

rspower/examples/test_make_jac.txt

@@ -3,6 +3,7 @@
 #include ../lib/idx_bus.txt
 #include ../lib/idx_bus.txt
 #include ../lib/idx_brch.txt
 #include ../lib/idx_brch.txt
 #include ../lib/make_y_bus.txt
 #include ../lib/make_y_bus.txt
+#include ../lib/dsbus_dv.txt
 #include ../lib/make_jac.txt
 #include ../lib/make_jac.txt
 
 
 jac = make_jac(baseMVA, bus, branch, gen);
 jac = make_jac(baseMVA, bus, branch, gen);

rspower/lib/dsbus_dv.txt

+//DSBUS_DV   Computes partial derivatives of power injection w.r.t. voltage.
+//
+//   The derivatives can be take with respect to polar or cartesian coordinates
+//   of voltage, depending on the 3rd argument.
+//
+//   [DSBUS_DVA, DSBUS_DVM] = DSBUS_DV(YBUS, V)
+//   [DSBUS_DVA, DSBUS_DVM] = DSBUS_DV(YBUS, V, 0)
+//
+//   Returns two matrices containing partial derivatives of the complex bus
+//   power injections w.r.t voltage angle and voltage magnitude, respectively
+//   (for all buses).
+//
+//   [DSBUS_DVR, DSBUS_DVI] = DSBUS_DV(YBUS, V, 1)
+//
+//   Returns two matrices containing partial derivatives of the complex bus
+//   power injections w.r.t the real and imaginary parts of voltage,
+//   respectively (for all buses).
+//
+//   If YBUS is a sparse matrix, the return values will be also. The following
+//   explains the expressions used to form the matrices:
+//
+//   S = diag(V) * conj(Ibus) = diag(conj(Ibus)) * V
+//
+//   Polar coordinates:
+//     Partials of V & Ibus w.r.t. voltage magnitudes
+//       dV/dVm = diag(V./abs(V))
+//       dI/dVm = Ybus * dV/dVm = Ybus * diag(V./abs(V))
+//
+//     Partials of V & Ibus w.r.t. voltage angles
+//       dV/dVa = j * diag(V)
+//       dI/dVa = Ybus * dV/dVa = Ybus * j * diag(V)
+//
+//     Partials of S w.r.t. voltage magnitudes
+//       dS/dVm = diag(V) * conj(dI/dVm) + diag(conj(Ibus)) * dV/dVm
+//              = diag(V) * conj(Ybus * diag(V./abs(V)))
+//                                       + conj(diag(Ibus)) * diag(V./abs(V))
+//
+//     Partials of S w.r.t. voltage angles
+//       dS/dVa = diag(V) * conj(dI/dVa) + diag(conj(Ibus)) * dV/dVa
+//              = diag(V) * conj(Ybus * j * diag(V))
+//                                       + conj(diag(Ibus)) * j * diag(V)
+//              = -j * diag(V) * conj(Ybus * diag(V))
+//                                       + conj(diag(Ibus)) * j * diag(V)
+//              = j * diag(V) * conj(diag(Ibus) - Ybus * diag(V))
+//
+//   Cartesian coordinates:
+//     Partials of V & Ibus w.r.t. real part of complex voltage
+//       dV/dVr = diag(ones(n,1))
+//       dI/dVr = Ybus * dV/dVr = Ybus
+//
+//     Partials of V & Ibus w.r.t. imaginary part of complex voltage
+//       dV/dVi = j * diag(ones(n,1))
+//       dI/dVi = Ybus * dV/dVi = Ybus * j
+//
+//     Partials of S w.r.t. real part of complex voltage
+//       dS/dVr = diag(V) * conj(dI/dVr) + diag(conj(Ibus)) * dV/dVr
+//              = diag(V) * conj(Ybus) + conj(diag(Ibus))
+//
+//     Partials of S w.r.t. imaginary part of complex voltage
+//       dS/dVi = diag(V) * conj(dI/dVi) + diag(conj(Ibus)) * dV/dVi
+//              = j * (conj(diag(Ibus)) - diag(V) conj(Ybus))
+//
+//   Examples:
+//       [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
+//       [dSbus_dVa, dSbus_dVm] = dSbus_dV(Ybus, V);
+//       [dSbus_dVr, dSbus_dVi] = dSbus_dV(Ybus, V, 1);
+//
+//   For more details on the derivations behind the derivative code used
+//   in MATPOWER information, see:
+//
+//   [TN2]  R. D. Zimmerman, "AC Power Flows, Generalized OPF Costs and
+//          their Derivatives using Complex Matrix Notation", MATPOWER
+//          Technical Note 2, February 2010. [Online]. Available:
+//          https://matpower.org/docs/TN2-OPF-Derivatives.pdf
+//          doi: 10.5281/zenodo.3237866
+//   [TN4]  B. Sereeter and R. D. Zimmerman, "AC Power Flows and their
+//          Derivatives using Complex Matrix Notation and Cartesian
+//          Coordinate Voltages," MATPOWER Technical Note 4, April 2018.
+//          [Online]. Available: https://matpower.org/docs/TN4-OPF-Derivatives-Cartesian.pdf
+//          doi: 10.5281/zenodo.3237909
+
+//   MATPOWER
+//   Copyright (c) 1996-2019, Power Systems Engineering Research Center (PSERC)
+//   by Ray Zimmerman, PSERC Cornell
+//   and Baljinnyam Sereeter, Delft University of Technology
+//
+//   This file is part of MATPOWER.
+//   Covered by the 3-clause BSD License (see LICENSE file for details).
+//   See https://matpower.org for more info.
+
+fn dSbus_dV(Ybus, V, vcart) {
+    n = length(V);
+    Ibus = Ybus * V;
+    diagV       = diag(V);
+    diagIbus    = diag(Ibus);
+    // RustScript use ~~ means not
+    if ~~vcart {
+        diagVnorm = diag(V./abs(V));
+    }
+    
+    if vcart {
+        dSbus_dV1 = conj(diagIbus) + diagV * conj(Ybus);        // dSbus/dVr
+        dSbus_dV2 = c(0,1) * (conj(diagIbus) - diagV * conj(Ybus)); // dSbus/dVi
+    } else {
+        dSbus_dV1 = c(0,1) * diagV * conj(diagIbus - Ybus * diagV);                     // dSbus/dVa
+        dSbus_dV2 = diagV * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm;    // dSbus/dVm
+    }
+    return dSbus_dV1;
+}
\ No newline at end of file

rspower/lib/make_jac.txt

@@ -8,5 +8,7 @@ fn make_jac(baseMVA, bus, branch, gen) {
     // make sure we use generator setpoint voltage for PV and slack buses
     // make sure we use generator setpoint voltage for PV and slack buses
     on = find(slice(gen, [0], [GEN_STATUS-1, GEN_STATUS]) > 0);  // which generators are on?
     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); // what buses are they at?
     gbus = get_multi(slice(gen, [0], [GEN_BUS-1, GEN_BUS]), on); // what buses are they at?
-    return gbus;
+
+    r = dSbus_dV(Ybus, V, 0.);
+    return r;
 }
 }
\ No newline at end of file