Improved the reactive power calculation for the generators in the pow…

…er flow simulation

src/GridCalEngine/Compilers/circuit_to_data.py

@@ -189,6 +189,12 @@ def get_load_data(circuit: MultiCircuit,
             data.active[ii] = elm.active
             data.cost[ii] = elm.Cost
 
+        # reactive power sharing data
+        if data.active[ii]:
+            bus_data.q_fixed[i] -= data.S[ii].imag
+            bus_data.ii_fixed[i] -= data.I[ii].imag
+            bus_data.b_fixed[i] -= data.Y[ii].imag
+
         data.C_bus_elm[i, ii] = 1
         ii += 1
 
@@ -209,6 +215,12 @@ def get_load_data(circuit: MultiCircuit,
             data.active[ii] = elm.active
             data.cost[ii] = elm.Cost
 
+        # reactive power sharing data
+        if data.active[ii]:
+            bus_data.q_fixed[i] += data.S[ii].imag
+            # bus_data.ii_fixed[i] += data.I[ii].imag
+            # bus_data.b_fixed[i] += data.Y[ii].imag
+
         data.C_bus_elm[i, ii] = 1
         ii += 1
 
@@ -243,6 +255,12 @@ def get_load_data(circuit: MultiCircuit,
             data.S[ii] += complex(elm.P, elm.Q)
             data.active[ii] = elm.active
 
+        # reactive power sharing data
+        if data.active[ii]:
+            bus_data.q_fixed[i] += data.S[ii].imag
+            bus_data.ii_fixed[i] += data.I[ii].imag
+            bus_data.b_fixed[i] += data.Y[ii].imag
+
         data.C_bus_elm[i, ii] = 1
         ii += 1
 
@@ -266,6 +284,12 @@ def get_load_data(circuit: MultiCircuit,
             data.active[ii] = elm.active
             data.cost[ii] = elm.Cost
 
+        # reactive power sharing data
+        if data.active[ii]:
+            bus_data.q_fixed[i] += data.S[ii].imag
+            bus_data.ii_fixed[i] += data.I[ii].imag
+            bus_data.b_fixed[i] += data.Y[ii].imag
+
         data.C_bus_elm[i, ii] = 1
         ii += 1
 
@@ -293,15 +317,15 @@ def get_shunt_data(
     :param t_idx:
     :param time_series:
     :param use_stored_guess:
+    :param control_remote_voltage:
     :return:
     """
-    devices = circuit.get_shunts()
 
     data = ds.ShuntData(nelm=circuit.get_shunt_like_device_number(),
                         nbus=circuit.get_bus_number())
 
     ii = 0
-    for k, elm in enumerate(devices):
+    for k, elm in enumerate(circuit.get_shunts()):
 
         i = bus_dict[elm.bus]
 
@@ -318,6 +342,10 @@ def get_shunt_data(
             data.active[k] = elm.active
             data.Y[k] = complex(elm.G, elm.B)
 
+        # reactive power sharing data
+        if data.active[ii]:
+            bus_data.b_fixed[i] += data.Y[ii].imag
+
         data.C_bus_elm[i, k] = 1
         ii += 1
 
@@ -382,6 +410,14 @@ def get_shunt_data(
                                         use_stored_guess=use_stored_guess,
                                         logger=logger)
 
+        # reactive power sharing data
+        if data.active[ii]:
+            if data.controllable[ii]:
+                bus_data.q_shared_total[i] += data.Y[ii].imag
+                data.q_share[ii] = data.Y[ii].imag
+            else:
+                bus_data.b_fixed[i] += data.Y[ii].imag
+
         data.C_bus_elm[i, ii] = 1
         ii += 1
 
@@ -536,6 +572,14 @@ def get_generator_data(
             data.qmin[k] = elm.Qmin
             data.qmax[k] = elm.Qmax
 
+        # reactive power sharing data
+        if data.active[k]:
+            if data.controllable[k]:
+                bus_data.q_shared_total[i] += data.p[k]
+                data.q_share[k] = data.p[k]
+            else:
+                bus_data.q_fixed[i] += data.get_q_at(k)
+
         data.C_bus_elm[i, k] = 1
 
     return data, gen_index_dict
@@ -696,6 +740,14 @@ def get_battery_data(
             data.qmin[k] = elm.Qmin
             data.qmax[k] = elm.Qmax
 
+        # reactive power sharing data
+        if data.active[k]:
+            if data.controllable[k]:
+                bus_data.q_shared_total[i] += data.p[k]
+                data.q_share[k] = data.p[k]
+            else:
+                bus_data.q_fixed[i] += data.get_q_at(k)
+
         data.C_bus_elm[i, k] = 1
 
     return data

src/GridCalEngine/DataStructures/bus_data.py

@@ -46,6 +46,15 @@ def __init__(self, nbus: int):
         self.areas: IntVec = np.empty(nbus, dtype=int)
         self.substations: IntVec = np.empty(nbus, dtype=int)
 
+        # This is the total value used to compute the q_share in generators, batteries and shunts
+        self.q_shared_total = np.zeros(nbus, dtype=float)
+
+        # This is the fixed amount of Q that is Subtrated from Qbus to compute
+        # the Q of each generator, battery and shunt controlling
+        self.q_fixed = np.zeros(nbus, dtype=float)
+        self.ii_fixed = np.zeros(nbus, dtype=float)  # same concept, but imag current from loads
+        self.b_fixed = np.zeros(nbus, dtype=float)  # same concept, but susceptance from shunts
+
         self.original_idx: IntVec = np.zeros(nbus, dtype=int)
 
     def slice(self, elm_idx: IntVec) -> "BusData":
@@ -76,6 +85,11 @@ def slice(self, elm_idx: IntVec) -> "BusData":
         data.areas = self.areas[elm_idx]
         data.substations = self.substations[elm_idx]
 
+        data.q_shared_total = self.q_shared_total[elm_idx]
+        data.q_fixed = self.q_fixed[elm_idx]
+        data.ii_fixed = self.ii_fixed[elm_idx]
+        data.b_fixed = self.b_fixed[elm_idx]
+
         data.original_idx = elm_idx
 
         return data
@@ -115,6 +129,11 @@ def copy(self) -> "BusData":
         data.areas = self.areas.copy()
         data.substations = self.substations.copy()
 
+        data.q_shared_total = self.q_shared_total.copy()
+        data.q_fixed = self.q_fixed.copy()
+        data.ii_fixed = self.ii_fixed.copy()
+        data.b_fixed = self.b_fixed.copy()
+
         data.original_idx = self.original_idx.copy()
 
         return data

src/GridCalEngine/DataStructures/generator_data.py

@@ -14,6 +14,8 @@
 # You should have received a copy of the GNU Lesser General Public License
 # along with this program; if not, write to the Free Software Foundation,
 # Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+from typing import Tuple
+
 import numpy as np
 import scipy.sparse as sp
 import GridCalEngine.Topology.topology as tp
@@ -46,6 +48,7 @@ def __init__(self, nelm: int, nbus: int):
 
         self.qmin: Vec = np.zeros(nelm, dtype=float)
         self.qmax: Vec = np.zeros(nelm, dtype=float)
+        self.q_share: Vec = np.zeros(nelm, dtype=float)
 
         # reliabilty
         self.mttf: Vec = np.zeros(nelm, dtype=float)
@@ -103,6 +106,7 @@ def slice(self, elm_idx: IntVec, bus_idx: IntVec):
 
         data.qmin = self.qmin[elm_idx]
         data.qmax = self.qmax[elm_idx]
+        data.q_share = self.q_share[elm_idx]
 
         data.mttf = self.mttf[elm_idx]
         data.mttr = self.mttr[elm_idx]
@@ -164,6 +168,7 @@ def copy(self):
 
         data.qmin = self.qmin.copy()
         data.qmax = self.qmax.copy()
+        data.q_share = self.q_share.copy()
 
         data.mttf = self.mttf.copy()
         data.mttr = self.mttr.copy()
@@ -217,6 +222,17 @@ def get_injections(self) -> CxVec:
         Q = pf_sign * self.p * np.sqrt((1.0 - pf2) / (pf2 + 1e-20))
         return self.p + 1.0j * Q
 
+    def get_q_at(self, i) -> float:
+        """
+
+        :param i:
+        :return:
+        """
+        pf2 = np.power(self.pf[i], 2.0)
+        pf_sign = (self.pf[i] + 1e-20) / np.abs(self.pf[i] + 1e-20)
+        Q = pf_sign * self.p[i] * np.sqrt((1.0 - pf2) / (pf2 + 1e-20))
+        return Q
+
     def get_Yshunt(self, seq: int = 1) -> CxVec:
         """
         Obtain the vector of shunt admittances of a given sequence per bus
@@ -322,4 +338,11 @@ def get_non_dispatchable_indices(self) -> IntVec:
         :return:
         """
         x = (~self.dispatchable * self.active).astype(int)
-        return np.where(x == 1)[0]
+        return np.where(x == 1)[0]
+
+    def get_controllable_and_not_controllable_indices(self) -> Tuple[IntVec, IntVec]:
+        """
+        Get the indices of controllable generators
+        :return: idx_controllable, idx_non_controllable
+        """
+        return np.where(self.controllable == 1)[0], np.where(self.controllable == 0)[0]

src/GridCalEngine/DataStructures/shunt_data.py

@@ -14,6 +14,7 @@
 # You should have received a copy of the GNU Lesser General Public License
 # along with this program; if not, write to the Free Software Foundation,
 # Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+from typing import Tuple
 import numpy as np
 import scipy.sparse as sp
 import GridCalEngine.Topology.topology as tp
@@ -46,6 +47,7 @@ def __init__(self, nelm: int, nbus: int):
 
         self.qmin: Vec = np.zeros(nelm, dtype=float)
         self.qmax: Vec = np.zeros(nelm, dtype=float)
+        self.q_share: Vec = np.zeros(nelm, dtype=float)
 
         self.cost: Vec = np.zeros(nelm, dtype=float)
 
@@ -88,6 +90,7 @@ def slice(self, elm_idx: IntVec, bus_idx: IntVec) -> "ShuntData":
 
         data.qmax = self.qmax[elm_idx]
         data.qmin = self.qmin[elm_idx]
+        data.q_share = self.q_share[elm_idx]
 
         data.cost = self.cost[elm_idx]
 
@@ -119,6 +122,7 @@ def copy(self) -> "ShuntData":
 
         data.qmax = self.qmax.copy()
         data.qmin = self.qmin.copy()
+        data.q_share = self.q_share.copy()
 
         data.cost = self.cost.copy()
 
@@ -162,6 +166,13 @@ def get_injections_per_bus(self) -> CxVec:
         """
         return self.C_bus_elm * (self.Y * self.active)
 
+    def get_fix_injections_per_bus(self) -> CxVec:
+        """
+        Get fixed Injections per bus
+        :return:
+        """
+        return self.C_bus_elm * (self.Y * self.active * (1 - self.controllable))
+
     def get_qmax_per_bus(self) -> Vec:
         """
         Get generator Qmax per bus
@@ -185,3 +196,10 @@ def get_bus_indices(self) -> IntVec:
         :return: array with the bus indices
         """
         return tp.get_csr_bus_indices(self.C_bus_elm.tocsr())
+
+    def get_controllable_and_not_controllable_indices(self) -> Tuple[IntVec, IntVec]:
+        """
+        Get the indices of controllable generators
+        :return: idx_controllable, idx_non_controllable
+        """
+        return np.where(self.controllable == 1)[0], np.where(self.controllable == 0)[0]

src/GridCalEngine/Simulations/PowerFlow/power_flow_worker.py

@@ -607,33 +607,40 @@ def power_flow_post_process(
     return Sfb, Stb, If, It, Vbranch, loading, losses, Sbus
 
 
-def split_reactive_power_into_devices(nc: NumericalCircuit, Qbus: Vec) -> Tuple[Vec, Vec, Vec]:
+def split_reactive_power_into_devices(nc: NumericalCircuit, Qbus: Vec, results: PowerFlowResults) -> None:
     """
     This function splits the reactive power of the power flow solution (nbus) into reactive power per device that
     is able to control reactive power as an injection (generators, batteries, shunts)
     :param nc: NumericalCircuit
     :param Qbus: Array of nodal reactive power (nbus)
-    :return: Qgen (ngen), Qbatt (nbatt), Qsh (nsh)
+    :param results: PowerFlowResults (values are written to it)
+    :return: Nothing, the results are set in the results object
     """
-    q_max_total = nc.generator_data.C_bus_elm @ (
-            nc.generator_data.qmax * nc.generator_data.active * nc.generator_data.controllable)
 
-    q_max_total += nc.battery_data.C_bus_elm @ (
-            nc.battery_data.qmax * nc.battery_data.active * nc.battery_data.controllable)
+    # generation
+    bus_idx_gen = nc.generator_data.get_bus_indices()
+    gen_q_share = nc.generator_data.q_share / (nc.bus_data.q_shared_total[bus_idx_gen] + 1e-20)
 
-    q_max_total += nc.shunt_data.C_bus_elm @ (
-            nc.shunt_data.qmax * nc.shunt_data.active * nc.shunt_data.controllable)
+    # batteries
+    bus_idx_bat = nc.battery_data.get_bus_indices()
+    batt_q_share = nc.battery_data.q_share / (nc.bus_data.q_shared_total[bus_idx_bat] + 1e-20)
 
-    bus_idx = nc.generator_data.get_bus_indices()
-    Qgen = Qbus[bus_idx] * nc.generator_data.qmax / (q_max_total[bus_idx] + 1e-20)
+    # shunts
+    bus_idx_sh = nc.shunt_data.get_bus_indices()
+    sh_q_share = nc.shunt_data.q_share / (nc.bus_data.q_shared_total[bus_idx_sh] + 1e-20)
 
-    bus_idx = nc.battery_data.get_bus_indices()
-    Qbatt = Qbus[bus_idx] * nc.battery_data.qmax / (q_max_total[bus_idx] + 1e-20)
+    # Fixed injection of reactive power
+    # Zip formul: S0 + np.conj(I0 + Y0 * Vm) * Vm
+    Vm = np.abs(results.voltage)
+    Qfix = nc.bus_data.q_fixed - (nc.bus_data.ii_fixed + nc.bus_data.b_fixed * Vm) * Vm
 
-    bus_idx = nc.shunt_data.get_bus_indices()
-    Qsh = Qbus[bus_idx] * nc.shunt_data.qmax / (q_max_total[bus_idx] + 1e-20)
+    # the remaining Q to share is the total Q computed (Qbus) minus the part that we know is fixed
+    Qvar = Qbus - Qfix
 
-    return Qgen, Qbatt, Qsh
+    # set the results
+    results.gen_q = Qvar[bus_idx_gen] * gen_q_share
+    results.battery_q = Qvar[bus_idx_bat] * batt_q_share
+    results.shunt_q = Qvar[bus_idx_sh] * sh_q_share
 
 
 def multi_island_pf_nc(nc: NumericalCircuit,
@@ -784,7 +791,7 @@ def multi_island_pf_nc(nc: NumericalCircuit,
     results.hvdc_losses = Losses_hvdc * nc.Sbase
 
     # do the reactive power partition and store the values
-    results.gen_q, results.battery_q, results.shunt_q = split_reactive_power_into_devices(nc=nc, Qbus=results.Sbus.imag)
+    split_reactive_power_into_devices(nc=nc, Qbus=results.Sbus.imag, results=results)
 
     return results