diff --git a/roseau/load_flow/_solvers.py b/roseau/load_flow/_solvers.py index 73f77ed2..6f695a90 100644 --- a/roseau/load_flow/_solvers.py +++ b/roseau/load_flow/_solvers.py @@ -1,16 +1,22 @@ import logging +import time import warnings from abc import ABC, abstractmethod -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Generic, TypeVar import numpy as np -from typing_extensions import Self from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode from roseau.load_flow.license import activate_license, get_license from roseau.load_flow.typing import JsonDict, Solver from roseau.load_flow.utils import find_stack_level -from roseau.load_flow_engine.cy_engine import CyAbstractSolver, CyBackwardForward, CyNewton, CyNewtonGoldstein +from roseau.load_flow_engine.cy_engine import ( + CyAbstractNewton, + CyAbstractSolver, + CyBackwardForward, + CyNewton, + CyNewtonGoldstein, +) logger = logging.getLogger(__name__) @@ -25,7 +31,10 @@ from roseau.load_flow.network import ElectricalNetwork -class AbstractSolver(ABC): +_CySolverT = TypeVar("_CySolverT", bound=CyAbstractSolver) + + +class AbstractSolver(ABC, Generic[_CySolverT]): """This is an abstract class for all the solvers.""" name: str | None = None @@ -38,10 +47,10 @@ def __init__(self, network: "ElectricalNetwork") -> None: The electrical network for which the load flow needs to be solved. """ self.network = network - self._cy_solver: CyAbstractSolver | None = None + self._cy_solver: _CySolverT | None = None @classmethod - def from_dict(cls, data: JsonDict, network: "ElectricalNetwork") -> Self: + def from_dict(cls, data: JsonDict, network: "ElectricalNetwork") -> "AbstractSolver": """AbstractSolver constructor from dict. Args: @@ -83,7 +92,38 @@ def solve_load_flow(self, max_iterations: int, tolerance: float) -> tuple[int, f lic = get_license() if lic is None: activate_license(key=None) - return self._cy_solver.solve_load_flow(max_iterations=max_iterations, tolerance=tolerance) + + start = time.perf_counter() + try: + iterations, residual = self._cy_solver.solve_load_flow(max_iterations=max_iterations, tolerance=tolerance) + except RuntimeError as e: + code, msg = e.args[0].split(" ", 1) + code = int(code) + if code == 0: + msg = f"The license cannot be validated. The detailed error message is {msg!r}" + exception_code = RoseauLoadFlowExceptionCode.LICENSE_ERROR + else: + msg, exception_code = self._parse_solver_error(code, msg) + logger.error(msg) + raise RoseauLoadFlowException(msg=msg, code=exception_code) from e + end = time.perf_counter() + + if iterations == max_iterations: + msg = ( + f"The load flow did not converge after {iterations} iterations. The norm of the " + f"residuals is {residual:5n}" + ) + logger.error(msg=msg) + raise RoseauLoadFlowException( + msg, RoseauLoadFlowExceptionCode.NO_LOAD_FLOW_CONVERGENCE, iterations, residual + ) + logger.debug(f"The load flow converged after {iterations} iterations and {end - start:.3n} s.") + return iterations, residual + + @abstractmethod + def _parse_solver_error(self, code: int, msg: str) -> tuple[str, RoseauLoadFlowExceptionCode]: + """Parse the solver's error code and message into a user-friendly message and an exception code.""" + raise NotImplementedError def reset_inputs(self) -> None: """Reset the input vector (which is used for the first step of the newton algorithm) to its initial value""" @@ -108,7 +148,10 @@ def params(self) -> JsonDict: return {} -class AbstractNewton(AbstractSolver, ABC): +_CyNewtonT = TypeVar("_CyNewtonT", bound=CyAbstractNewton) + + +class AbstractNewton(AbstractSolver[_CyNewtonT], ABC): """This is an abstract class for all the Newton-Raphson solvers.""" DEFAULT_TAPE_OPTIMIZATION: bool = True @@ -127,23 +170,66 @@ def __init__(self, network: "ElectricalNetwork", optimize_tape: bool = DEFAULT_T super().__init__(network=network) self.optimize_tape = optimize_tape + def _parse_solver_error(self, code: int, msg: str) -> tuple[str, RoseauLoadFlowExceptionCode]: + assert code == 1, f"Unexpected error code {code} for a Newton solver." + zero_elements_index, inf_elements_index = self.analyse_jacobian() + if inf_elements_index: + inf_elements = [self.network._elements[i] for i in inf_elements_index] + printable_elements = ", ".join(f"{type(e).__name__}({e.id!r})" for e in inf_elements) + msg += f"The problem seems to come from the elements [{printable_elements}] that induce infinite values." + power_load = False + flexible_load = False + for inf_element in inf_elements: + load = self.network.loads.get(inf_element.id) + if load is inf_element and load.type == "power": + power_load = True + if load.is_flexible: + flexible_load = True + break + if power_load: + msg += " This might be caused by a bad potential initialization of a power load" + if flexible_load: + msg += ", or by flexible loads with very high alpha or incorrect flexible parameters voltages." + raise RoseauLoadFlowException(msg=msg, code=RoseauLoadFlowExceptionCode.NAN_VALUE) + elif zero_elements_index: + zero_elements = [self.network._elements[i] for i in zero_elements_index] + printable_elements = ", ".join(f"{type(e).__name__}({e.id!r})" for e in zero_elements) + msg += ( + f"The problem seems to come from the elements [{printable_elements}] that have at " + f"least one disconnected phase." + ) + return msg, RoseauLoadFlowExceptionCode.BAD_JACOBIAN + + # Debugging methods + # ----------------- def save_matrix(self, prefix: str) -> None: - """Output files of the jacobian and vector matrix of the first newton step. Those files can be used to launch an - eigen solver benchmark (see https://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html) + """Output files of the jacobian and vector matrix of the first newton step. + + Those files can be used to launch an eigen solver benchmark + (see https://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html) Args: prefix: - The prefix of the name of the files. They will be output as prefix.mtx and prefix_m.mtx to follow Eigen - solver benchmark convention. + The prefix of the name of the files. They will be output as `prefix.mtx` and + `prefix_m.mtx` to follow Eigen solver benchmark convention. """ self._cy_solver.save_matrix(prefix) - def current_jacobian(self) -> np.ndarray: - """Show the jacobian of the current iteration (useful for debugging)""" + def current_jacobian(self) -> np.ndarray[tuple[int, int], np.dtype[np.float64]]: + """Get the jacobian of the current iteration (useful for debugging).""" return self._cy_solver.current_jacobian() + def analyse_jacobian(self) -> tuple[list[int], list[int]]: + """Analyse the jacobian to try to understand why it is singular. -class Newton(AbstractNewton): + Returns: + The vector of elements associated with a column of zeros and the vector of elements + associated with a line which contains a NaN. + """ + return self._cy_solver.analyse_jacobian() + + +class Newton(AbstractNewton[CyNewton]): """The classical Newton-Raphson algorithm.""" name = "newton" @@ -168,7 +254,7 @@ def update_network(self, network: "ElectricalNetwork") -> None: self._cy_solver = CyNewton(network=network._cy_electrical_network, optimize_tape=self.optimize_tape) -class NewtonGoldstein(AbstractNewton): +class NewtonGoldstein(AbstractNewton[CyNewtonGoldstein]): """The Newton-Raphson algorithm with the Goldstein and Price linear search. It has better stability than the classical Newton-Raphson, without losing performance. """ @@ -229,7 +315,7 @@ def params(self) -> JsonDict: return {"m1": self.m1, "m2": self.m2} -class BackwardForward(AbstractSolver): +class BackwardForward(AbstractSolver[CyBackwardForward]): """A backward-forward implementation, less stable than Newton-Raphson based algorithms, but can be more performant in some cases. """ @@ -262,3 +348,7 @@ def solve_load_flow(self, max_iterations: int, tolerance: float) -> tuple[int, f logger.error(msg) raise RoseauLoadFlowException(msg=msg, code=RoseauLoadFlowExceptionCode.NO_BACKWARD_FORWARD) return super().solve_load_flow(max_iterations, tolerance) + + def _parse_solver_error(self, code: int, msg: str) -> tuple[str, RoseauLoadFlowExceptionCode]: + assert code == 2, f"Unexpected error code {code} for a Backward-Forward solver." + return msg, RoseauLoadFlowExceptionCode.NO_BACKWARD_FORWARD diff --git a/roseau/load_flow/network.py b/roseau/load_flow/network.py index 70113846..6d7330c9 100644 --- a/roseau/load_flow/network.py +++ b/roseau/load_flow/network.py @@ -6,14 +6,13 @@ import logging import re import textwrap -import time import warnings from collections.abc import Iterable, Mapping from importlib import resources from itertools import chain from math import nan from pathlib import Path -from typing import TYPE_CHECKING, NoReturn, TypeVar +from typing import TYPE_CHECKING, TypeVar import geopandas as gpd import numpy as np @@ -518,10 +517,9 @@ def solve_load_flow( network (e.g. ``print(net.res_buses``). To get the results for a specific element, use the `res_` properties on the element (e.g. ``print(net.buses["bus1"].res_potentials)``. - You need to activate the license before calling this method. Alternatively you may set the - environment variable ``ROSEAU_LOAD_FLOW_LICENSE_KEY`` to your license key and it will be - picked automatically when calling this method. See the :ref:`license` page for more - information. + You need to activate the license before calling this method. You may set the environment + variable ``ROSEAU_LOAD_FLOW_LICENSE_KEY`` to your license key and it will be picked + automatically when calling this method. See the :ref:`license` page for more information. Args: max_iterations: @@ -537,9 +535,14 @@ def solve_load_flow( solver: The name of the solver to use for the load flow. The options are: - - ``'newton'``: the classical Newton-Raphson solver. - - ``'newton_goldstein'``: the Newton-Raphson solver with the Goldstein and - Price linear search. + + - ``newton``: The classical *Newton-Raphson* method. + - ``newton_goldstein``: The *Newton-Raphson* method with the *Goldstein and Price* + linear search algorithm. It generally has better convergence properties than the + classical Newton-Raphson method. This is the default. + - ``backward_forward``: the *Backward-Forward Sweep* method. It usually executes + faster than the other approaches but may exhibit weaker convergence properties. It + does not support meshed networks or floating neutrals. solver_params: A dictionary of parameters used by the solver. Available parameters depend on the @@ -563,25 +566,8 @@ def solve_load_flow( if not warm_start: self._reset_inputs() - start = time.perf_counter() - try: - iterations, residual = self._solver.solve_load_flow(max_iterations=max_iterations, tolerance=tolerance) - except RuntimeError as e: - self._handle_error(e) - - end = time.perf_counter() - - if iterations == max_iterations: - msg = ( - f"The load flow did not converge after {iterations} iterations. The norm of the residuals is " - f"{residual:5n}" - ) - logger.error(msg=msg) - raise RoseauLoadFlowException( - msg, RoseauLoadFlowExceptionCode.NO_LOAD_FLOW_CONVERGENCE, iterations, residual - ) + iterations, residual = self._solver.solve_load_flow(max_iterations=max_iterations, tolerance=tolerance) - logger.debug(f"The load flow converged after {iterations} iterations and {end - start:.3n} s.") self._no_results = False # Lazily update the results of the elements @@ -594,47 +580,6 @@ def solve_load_flow( return iterations, residual - def _handle_error(self, e: RuntimeError) -> NoReturn: - msg = e.args[0] - if msg.startswith("0 "): - msg = f"The license cannot be validated. The detailed error message is {msg[2:]!r}" - logger.error(msg) - raise RoseauLoadFlowException(msg=msg, code=RoseauLoadFlowExceptionCode.LICENSE_ERROR) from e - elif msg.startswith("1 "): - msg = msg[2:] - zero_elements_index, inf_elements_index = self._solver._cy_solver.analyse_jacobian() - if inf_elements_index: - inf_elements = [self._elements[i] for i in inf_elements_index] - printable_elements = ", ".join(f"{type(e).__name__}({e.id!r})" for e in inf_elements) - msg += ( - f"The problem seems to come from the elements [{printable_elements}] that induce infinite values." - ) - power_load = False - flexible_load = False - for inf_element in inf_elements: - if isinstance(inf_element, PowerLoad): - power_load = True - if inf_element.is_flexible: - flexible_load = True - if power_load: - msg += " This might be caused by a bad potential initialization of a power load" - if flexible_load: - msg += ", or by flexible loads with very high alpha or incorrect flexible parameters voltages." - raise RoseauLoadFlowException(msg=msg, code=RoseauLoadFlowExceptionCode.NAN_VALUE) - elif zero_elements_index: - zero_elements = [self._elements[i] for i in zero_elements_index] - printable_elements = ", ".join(f"{type(e).__name__}({e.id!r})" for e in zero_elements) - msg += ( - f"The problem seems to come from the elements [{printable_elements}] that have at least one " - f"disconnected phase. " - ) - logger.error(msg) - raise RoseauLoadFlowException(msg=msg, code=RoseauLoadFlowExceptionCode.BAD_JACOBIAN) from e - else: - assert msg.startswith("2 ") - msg = msg[2:] - raise RoseauLoadFlowException(msg=msg, code=RoseauLoadFlowExceptionCode.NO_BACKWARD_FORWARD) from e - # # Properties to access the load flow results as dataframes #