Merge branch 'refs/heads/devel'

.env

@@ -1 +1 @@
-PYTHONPATH=src
+PYTHONPATH=./src:./src/tests

src/GridCalEngine/Compilers/circuit_to_data.py

@@ -177,23 +177,12 @@ def set_bus_control_voltage_hvdc(i: int,
     :param logger: Logger
     """
     if bus_data.bus_types[i] != BusMode.Slack_tpe.value:  # if it is not Slack
-        if remote_control and j > -1 and j != i:
-            # remove voltage control
-            # bus_data.bus_types[j] = BusMode.PQV_tpe.value  # remote bus to PQV type
-            # bus_data.set_bus_mode(j, BusMode.PQV_tpe)
-            bus_data.is_p_controlled[j] = True
-            bus_data.is_q_controlled[j] = True
-            bus_data.is_vm_controlled[j] = True
-
-            # bus_data.bus_types[i] = BusMode.P_tpe.value  # local bus to P type
-            # bus_data.set_bus_mode(i, BusMode.P_tpe)
-            bus_data.is_p_controlled[i] = True
-        else:
-            # local voltage control
-            # bus_data.bus_types[i] = BusMode.PV_tpe.value  # set as PV
-            # bus_data.set_bus_mode(i, BusMode.PV_tpe)
-            bus_data.is_p_controlled[i] = True
-            bus_data.is_vm_controlled[i] = True
+        # local voltage control
+        bus_data.bus_types[i] = BusMode.PV_tpe.value  # set as PV
+        # bus_data.set_bus_mode(i, BusMode.PV_tpe)
+        bus_data.is_p_controlled[i] = True
+        bus_data.is_vm_controlled[i] = True
+        bus_data.is_q_controlled[i] = True
 
     if not use_stored_guess:
         if not bus_voltage_used[i]:

src/GridCalEngine/Simulations/PowerFlow/Formulations/pf_generalized_formulation.py

@@ -383,7 +383,7 @@ def calcYbus(Cf, Ct, Yshunt_bus: CxVec,
     return Ybus.tocsc()
 
 
-@njit(cache=True)
+# @njit(cache=True)
 def calcSf(k: IntVec, V: CxVec, F: IntVec, T: IntVec,
            R: Vec, X: Vec, G: Vec, B: Vec, m: Vec, tau: Vec, vtap_f: Vec, vtap_t: Vec):
     """
@@ -528,6 +528,11 @@ def __init__(self, V0: CxVec, S0: CxVec, I0: CxVec, Y0: CxVec,
         self.is_vm_controlled = nc.bus_data.is_vm_controlled.copy()
         self.is_va_controlled = nc.bus_data.is_va_controlled.copy()
 
+        # HVDC LOOP
+        for k in range(self.nc.hvdc_data.nelm):
+            self.nc.bus_data.is_q_controlled[self.nc.hvdc_data.F[k]] = True
+            self.nc.bus_data.is_q_controlled[self.nc.hvdc_data.T[k]] = True
+
         # Controllable Branch Indices
         self.u_cbr_m = np.zeros(0, dtype=int)
         self.u_cbr_tau = np.zeros(0, dtype=int)
@@ -2253,6 +2258,9 @@ def _analyze_hvdc_controls(self) -> None:
         # HVDC LOOP
         for k in range(self.nc.hvdc_data.nelm):
             # hvdc.append(k)
+            # self.nc.bus_data.is_q_controlled[self.nc.hvdc_data.bus_f[k]] = True
+            # self.nc.bus_data.is_q_controlled[self.nc.hvdc_data.bus_t[k]] = True
+            # self.nc.bus_data.bus_types[self.nc.hvdc_data.bus_f[k]] = BusMode.PQV
             if self.nc.hvdc_data.control_mode[k] == HvdcControlType.type_0_free:
                 hvdc_droop_idx.append(k)
 
@@ -2504,6 +2512,18 @@ def compute_f(self, x: Vec) -> Vec:
             Qt_cbr - self.cbr_qt_set
         ]
 
+        # Print index blocks of f
+        print('Lengths: ')
+        print(len(self.i_k_p), len(self.i_k_q), len(loss_vsc), len(loss_hvdc), len(inj_hvdc),
+              len(Pf_cbr), len(Pt_cbr), len(Qf_cbr), len(Qt_cbr))
+
+        print('Pf set: ', self.cbr_pf_set)
+        print('f errors: ')
+        # Get indices of troublesome values
+        for i, ff in enumerate(_f):
+            if abs(ff) > 0.5:
+                print(i, ff)
+
         return _f
 
     def check_error(self, x: Vec) -> Tuple[float, Vec]:

src/tests/test_power_flow.py

@@ -14,6 +14,8 @@
 from GridCalEngine.Compilers.circuit_to_data import compile_numerical_circuit_at
 import GridCalEngine.api as gce
 
+SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
+
 
 def test_ieee_grids():
     """
@@ -541,12 +543,9 @@ def test_hvdc_all_methods() -> None:
     """
     Checks that the HVDC logic is working for all power flow methods
     """
-    fname = os.path.join('data', 'grids', '8_nodes_2_islands_hvdc.gridcal')
-
+    fname = os.path.join(SCRIPT_DIR, 'data', 'grids', '8_nodes_2_islands_hvdc.gridcal')
     grid = gce.open_file(fname)
 
-    grid.hvdc_lines[0].Pset = 10  # this is what we'll check later
-
     for solver_type in [SolverType.NR,
                         SolverType.LM,
                         SolverType.PowellDogLeg,
@@ -585,6 +584,8 @@ def test_hvdc_all_methods() -> None:
 
         assert res.converged
         assert res.Pf_hvdc[0] == 10.0
+        assert np.isclose(abs(res.voltage[6]), 1.01111, atol=1e-4)
+        assert np.isclose(abs(res.voltage[1]), 1.02222, atol=1e-4)
 
     # repeat forcing to use the special formulations
     for solver_type in [SolverType.NR,
@@ -623,7 +624,10 @@ def test_hvdc_all_methods() -> None:
 
         assert res.converged
         assert res.Pf_hvdc[0] == 10.0
+        assert np.isclose(abs(res.voltage[6]), 1.01111, atol=1e-4)
+        assert np.isclose(abs(res.voltage[1]), 1.02222, atol=1e-4)
 
 
 if __name__ == "__main__":
-    test_power_flow_12bus_acdc()
+    # test_power_flow_12bus_acdc()
+    test_hvdc_all_methods()

src/tests/test_raw_cgmes_cross_roundtrip.py

@@ -114,33 +114,37 @@ def run_raw_to_cgmes(import_path: str | list[str],
         logger_nc.print()
 
         # FOR DEBUG
+        adm1 = nc1.get_admittance_matrices()
+        adm2 = nc2.get_admittance_matrices()
+        Sbus1 = nc1.get_power_injections()
+        Sbus2 = nc2.get_power_injections()
         print('Buses')
-        print(nc1.bus_names)
-        print(nc2.bus_names)
+        print(nc1.bus_data.names)
+        print(nc2.bus_data.names)
         print('Loads')
-        print(nc1.load_names)
-        print(nc2.load_names)
+        print(nc1.load_data.names)
+        print(nc2.load_data.names)
         print('Gens')
-        print(nc1.generator_names)
-        print(nc2.generator_names)
+        print(nc1.generator_data.names)
+        print(nc2.generator_data.names)
         print('Sbus1')
-        print(nc1.Sbus)
+        print(Sbus1)
         print('Sbus2')
-        print(nc2.Sbus)
+        print(Sbus2)
         print('S_diff')
-        print(nc2.Sbus - nc1.Sbus)
+        print(Sbus1 - Sbus2)
         print('Y1')
-        print(nc1.Ybus.A)
+        print(adm1.Ybus.A)
         print('Y2')
-        print(nc2.Ybus.A)
-        Y_diff = nc2.Ybus.A - nc1.Ybus.A
+        print(adm2.Ybus.A)
+        Y_diff = adm2.Ybus.A - adm1.Ybus.A
         print('Y_diff', Y_diff)
         # mask = Y_diff != 0
-        mask = ~np.isclose(nc1.Ybus.A, nc2.Ybus.A, atol=1e-4, rtol=0)
+        mask = ~np.isclose(adm1.Ybus.A, adm2.Ybus.A, atol=1e-4, rtol=0)
         # print('mask \n', mask)
-        print("Y1_elements", nc1.Ybus.A[mask])
-        print("Y2_elements", nc2.Ybus.A[mask])
-        print("Y_diff", nc1.Ybus.A[mask] - nc2.Ybus.A[mask])
+        print("Y1_elements", adm1.Ybus.A[mask])
+        print("Y2_elements", adm2.Ybus.A[mask])
+        print("Y_diff", adm1.Ybus.A[mask] - adm2.Ybus.A[mask])
 
     assert np.allclose(np.abs(pf1_res.voltage), np.abs(pf2_res.voltage), atol=1e-5)
 

src/trunk/acdc_pf/generalized_wip/try_large_grid.py

@@ -0,0 +1,134 @@
+# This Source Code Form is subject to the terms of the Mozilla Public
+# License, v. 2.0. If a copy of the MPL was not distributed with this
+# file, You can obtain one at https://mozilla.org/MPL/2.0/.  
+# SPDX-License-Identifier: MPL-2.0
+from typing import Tuple
+import os
+import pandas as pd
+import numpy as np
+from GridCalEngine.basic_structures import Logger
+from GridCalEngine.IO.file_handler import FileOpen
+from GridCalEngine.Simulations.PowerFlow.power_flow_worker import PowerFlowOptions
+from GridCalEngine.Simulations.PowerFlow.power_flow_options import SolverType
+from GridCalEngine.Simulations.PowerFlow.power_flow_results import NumericPowerFlowResults
+import GridCalEngine.api as gce
+from GridCalEngine.Simulations.PowerFlow.Formulations.pf_generalized_formulation import PfGeneralizedFormulation
+from GridCalEngine.Simulations.PowerFlow.Formulations.pf_advanced_formulation import PfAdvancedFormulation
+from GridCalEngine.Simulations.PowerFlow.NumericalMethods.newton_raphson_fx import newton_raphson_fx
+
+SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
+TEST_FOLDER = os.path.join(SCRIPT_DIR, "..", "..", "..", "tests")
+
+
+def solve_generalized(grid: gce.MultiCircuit,
+                      options: PowerFlowOptions,
+                      generalized: bool=True) -> Tuple[PfGeneralizedFormulation, NumericPowerFlowResults]:
+    """
+
+    :param grid:
+    :param options:
+    :param generalized:
+    :return:
+    """
+    nc = gce.compile_numerical_circuit_at(
+        grid,
+        t_idx=None,
+        apply_temperature=False,
+        branch_tolerance_mode=gce.BranchImpedanceMode.Specified,
+        opf_results=None,
+        use_stored_guess=False,
+        bus_dict=None,
+        areas_dict=None,
+        control_taps_modules=options.control_taps_modules,
+        control_taps_phase=options.control_taps_phase,
+        control_remote_voltage=options.control_remote_voltage,
+    )
+
+    islands = nc.split_into_islands(consider_hvdc_as_island_links=True)
+    logger = Logger()
+
+    island = islands[0]
+
+    Vbus = island.bus_data.Vbus
+    S0 = island.get_power_injections_pu()
+    I0 = island.get_current_injections_pu()
+    Y0 = island.get_admittance_injections_pu()
+    Qmax_bus, Qmin_bus = island.get_reactive_power_limits()
+
+    if generalized:
+        problem = PfGeneralizedFormulation(V0=Vbus,
+                                       S0=S0,
+                                       I0=I0,
+                                       Y0=Y0,
+                                       Qmin=Qmin_bus,
+                                       Qmax=Qmax_bus,
+                                       nc=island,
+                                       options=options,
+                                       logger=logger)
+    else:
+        problem = PfAdvancedFormulation(V0=Vbus,
+                                    S0=S0,
+                                    I0=I0,
+                                    Y0=Y0,
+                                    Qmin=Qmin_bus,
+                                    Qmax=Qmax_bus,
+                                    nc=island,
+                                    options=options,
+                                    logger=logger)
+
+    solution = newton_raphson_fx(problem=problem,
+                                 tol=options.tolerance,
+                                 max_iter=options.max_iter,
+                                 trust=options.trust_radius,
+                                 verbose=options.verbose,
+                                 logger=logger)
+
+    logger.print("Logger")
+
+    return problem, solution
+
+
+def try_power_flow() -> None:
+    """
+    Check that a transformer can regulate the voltage at a bus
+    """
+    # fname = os.path.join('/Users', 'josep', 'Documents', 'Grids', '202206011015_vsc.gridcal')
+    fname = os.path.join('/Users', 'josep', 'Documents', 'Grids', '202206011015_original.gridcal')
+    # fname = os.path.join('/Users', 'josep', 'Documents', 'Grids', '202206011015_off_shunts.gridcal')
+    # fname = os.path.join('/Users', 'josep', 'Documents', 'Grids', '202206011015_on_shunts.gridcal')
+    # fname = os.path.join('/Users', 'josep', 'Documents', 'Grids', '202206011015_no_shunts.gridcal')
+    # fname = os.path.join('/Users', 'josep', 'Documents', 'Grids', '202206011015_vsc_v2.gridcal')
+
+    grid = gce.open_file(fname)
+
+    options = PowerFlowOptions(solver_type=gce.SolverType.NR,
+                               verbose=1,
+                               control_q=False,
+                               retry_with_other_methods=False,
+                               initialize_with_existing_solution=False,
+                               control_taps_phase=True,
+                               control_taps_modules=False,
+                               tolerance=1e-5,
+                               max_iter=80)
+
+    problem, solution = solve_generalized(grid=grid, options=options, generalized=True)
+    # problem, solution = solve_generalized(grid=grid, options=options, generalized=False)
+
+    print(solution.elapsed)
+
+    n_act = 0
+    n_inact = 0
+    for bus in grid.buses:
+        if bus.active:
+            n_act += 1
+        else:
+            n_inact += 1
+
+    print(n_act, n_inact)
+    print(solution.norm_f)
+
+    assert solution.converged
+
+
+if __name__ == "__main__":
+    try_power_flow()