Skip to content

Commit

Permalink
Centralize the low-level solver code
Browse files Browse the repository at this point in the history 
For use with rlf and rlfs
  • Loading branch information
@alihamdan
alihamdan committed Jan 15, 2025
1 parent 5dfe47f commit e730d8b
Show file tree
Hide file tree
Showing 2 changed files with 120 additions and 85 deletions.
124 changes: 107 additions & 17 deletions roseau/load_flow/_solvers.py
View file
Original file line number Diff line number Diff line change
@@ -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__)

Expand All @@ -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
Expand All @@ -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:
Expand Down Expand Up @@ -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"""
Expand All @@ -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
Expand All @@ -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"
Expand All @@ -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.
"""
Expand Down Expand Up @@ -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.
"""
Expand Down Expand Up @@ -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
81 changes: 13 additions & 68 deletions roseau/load_flow/network.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
#
Expand Down

0 comments on commit e730d8b

Please sign in to comment.