Skip to content

Commit

Permalink
Merge pull request #484 from oemof/fix/newton-convergence-critereons
Browse files Browse the repository at this point in the history 
Remove old newton function and use newton_with_kwargs instead
  • Loading branch information
@fwitte
fwitte authored Jan 29, 2024
2 parents 1b581b6 + 33c076c commit 3fddc8d
Show file tree
Hide file tree
Showing 4 changed files with 67 additions and 157 deletions.
58 changes: 29 additions & 29 deletions docs/modules/ude.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,15 +65,15 @@ equation in a function which returns the residual value of the equation.

.. code-block:: python
>>> def my_ude(self):
... return self.conns[0].m.val_SI - self.conns[1].m.val_SI ** 2
>>> def my_ude(ude):
... return ude.conns[0].m.val_SI - ude.conns[1].m.val_SI ** 2
.. note::

The function must only take one parameter, i.e. the UserDefinedEquation
class instance. The **name of the parameter is arbitrary**. We will use
:code:`self` in this example. It serves to access some important parameters
:code:`ude` in this example. It serves to access some important parameters
of the equation:

- connections required in the equation
Expand All @@ -94,7 +94,7 @@ composition of every connection. The derivatives have to be passed to the
Jacobian. In order to do this, we create a function that updates the values
inside the Jacobian of the :code:`UserDefinedEquation` and returns it:

- :code:`self.jacobian` is a dictionary containing numpy arrays for every
- :code:`ude.jacobian` is a dictionary containing numpy arrays for every
connection required by the :code:`UserDefinedEquation`.
- derivatives to **mass flow** are placed in the first element of the numpy
array (**index 0**)
Expand All @@ -114,13 +114,13 @@ derivatives to mass flow are not zero.

.. code-block:: python
>>> def my_ude_deriv(self):
... c0 = self.conns[0]
... c1 = self.conns[1]
>>> def my_ude_deriv(ude):
... c0 = ude.conns[0]
... c1 = ude.conns[1]
... if c0.m.is_var:
... self.jacobian[c0.m.J_col] = 1
... ude.jacobian[c0.m.J_col] = 1
... if c1.m.is_var:
... self.jacobian[c1.m.J_col] = -2 * self.conns[1].m.val_SI
... ude.jacobian[c1.m.J_col] = -2 * ude.conns[1].m.val_SI
Now we can create our instance of the :code:`UserDefinedEquation` and add it to
the network. The class requires four mandatory arguments to be passed:
Expand Down Expand Up @@ -185,21 +185,21 @@ respectively to calculate the partial derivatives.
>>> from tespy.tools.fluid_properties import dT_mix_dph
>>> from tespy.tools.fluid_properties import dT_mix_pdh
>>> def my_ude_deriv(self):
... c0 = self.conns[0]
... c1 = self.conns[1]
>>> def my_ude_deriv(ude):
... c0 = ude.conns[0]
... c1 = ude.conns[1]
... if c0.m.is_var:
... self.jacobian[c0.m.J_col] = 1 / self.conns[0].m.val_SI
... ude.jacobian[c0.m.J_col] = 1 / ude.conns[0].m.val_SI
... if c0.p.is_var:
... self.jacobian[c0.p.J_col] = - 2 / self.conns[0].p.val_SI
... ude.jacobian[c0.p.J_col] = - 2 / ude.conns[0].p.val_SI
... T = c1.calc_T()
... if c1.p.is_var:
... self.jacobian[c1.p.J_col] = (
... ude.jacobian[c1.p.J_col] = (
... dT_mix_dph(c1.p.val_SI, c1.h.val_SI, c1.fluid_data, c1.mixing_rule)
... * 0.5 / (T ** 0.5)
... )
... if c1.h.is_var:
... self.jacobian[c1.h.J_col] = (
... ude.jacobian[c1.h.J_col] = (
... dT_mix_pdh(c1.p.val_SI, c1.h.val_SI, c1.fluid_data, c1.mixing_rule)
... * 0.5 / (T ** 0.5)
... )
Expand Down Expand Up @@ -265,27 +265,27 @@ instance must therefore be changed as below.
>>> from tespy.tools.fluid_properties import h_mix_pQ
>>> from tespy.tools.fluid_properties import dh_mix_dpQ
>>> def my_ude(self):
... a = self.params['a']
... b = self.params['b']
... c0 = self.conns[0]
... c1 = self.conns[1]
>>> def my_ude(ude):
... a = ude.params['a']
... b = ude.params['b']
... c0 = ude.conns[0]
... c1 = ude.conns[1]
... return (
... a * (c1.h.val_SI - c0.h.val_SI) -
... (c1.h.val_SI - h_mix_pQ(c0.p.val_SI, b, c0.fluid_data))
... )
>>> def my_ude_deriv(self):
... a = self.params['a']
... b = self.params['b']
... c0 = self.conns[0]
... c1 = self.conns[1]
>>> def my_ude_deriv(ude):
... a = ude.params['a']
... b = ude.params['b']
... c0 = ude.conns[0]
... c1 = ude.conns[1]
... if c0.p.is_var:
... self.jacobian[c0.p.J_col] = dh_mix_dpQ(c0.p.val_SI, b, c0.fluid_data)
... ude.jacobian[c0.p.J_col] = dh_mix_dpQ(c0.p.val_SI, b, c0.fluid_data)
... if c0.h.is_var:
... self.jacobian[c0.h.J_col] = -a
... ude.jacobian[c0.h.J_col] = -a
... if c1.p.is_var:
... self.jacobian[c1.p.J_col] = a - 1
... ude.jacobian[c1.p.J_col] = a - 1
>>> ude = UserDefinedEquation(
... 'my ude', my_ude, my_ude_deriv, [c1, c2], params={'a': 0.5, 'b': 1}
Expand Down
22 changes: 16 additions & 6 deletions src/tespy/components/component.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@
from tespy.tools.helpers import _numeric_deriv
from tespy.tools.helpers import bus_char_derivative
from tespy.tools.helpers import bus_char_evaluation
from tespy.tools.helpers import newton
from tespy.tools.helpers import newton_with_kwargs


class Component:
Expand Down Expand Up @@ -642,11 +642,21 @@ def calc_bus_expr(self, bus):
if b['base'] == 'component':
return abs(comp_val / b['P_ref'])
else:
bus_value = newton(
bus_char_evaluation,
bus_char_derivative,
[comp_val, b['P_ref'], b['char']], 0,
val0=b['P_ref'], valmin=-1e15, valmax=1e15)
kwargs = {
"function": bus_char_evaluation,
"parameter": "bus_value",
"component_value": comp_val,
"reference_value": b["P_ref"],
"char_func": b["char"]
}
bus_value = newton_with_kwargs(
derivative=bus_char_derivative,
target_value=0,
val0=b['P_ref'],
valmin=-1e15,
valmax=1e15,
**kwargs
)
return bus_value / b['P_ref']

def calc_bus_efficiency(self, bus):
Expand Down
120 changes: 9 additions & 111 deletions src/tespy/tools/helpers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -471,7 +471,7 @@ def _numeric_deriv(obj, func, dx, conn=None, **kwargs):
return deriv


def bus_char_evaluation(params, bus_value):
def bus_char_evaluation(component_value, char_func, reference_value, bus_value, **kwargs):
r"""
Calculate the value of a bus.
Expand All @@ -497,128 +497,23 @@ def bus_char_evaluation(params, bus_value):
{f\left(\frac{\dot{E}_\mathrm{bus}}
{\dot{E}_\mathrm{bus,ref}}\right)}
"""
comp_value = params[0]
reference_value = params[1]
char_func = params[2]
return bus_value - comp_value / char_func.evaluate(
bus_value / reference_value)
return bus_value - component_value / char_func.evaluate(
bus_value / reference_value
)


def bus_char_derivative(params, bus_value):
def bus_char_derivative(component_value, char_func, reference_value, bus_value, **kwargs):
"""Calculate derivative for bus char evaluation."""
reference_value = params[1]
char_func = params[2]
d = 1e-3
return (1 - (
1 / char_func.evaluate((bus_value + d) / reference_value) -
1 / char_func.evaluate((bus_value - d) / reference_value)
) / (2 * d))


def newton(func, deriv, params, y, **kwargs):
r"""
Find zero crossings with 1-D newton algorithm.
Parameters
----------
func : function
Function to find zero crossing in,
:math:`0=y-func\left(x,\text{params}\right)`.
deriv : function
First derivative of the function.
params : list
Additional parameters for function, optional.
y : float
Target function value.
val0 : float
Starting value, default: val0=300.
valmin : float
Lower value boundary, default: valmin=70.
valmax : float
Upper value boundary, default: valmax=3000.
max_iter : int
Maximum number of iterations, default: max_iter=10.
tol_rel : float
Maximum relative tolerance :math:`|\frac{y - f(x)}{f(x)}|`, default
value: 1e-6.
tol_abs : float
Maximum absolute tolerance :math:`|y - f(x)|`, default value: 1e-6.
tol_mode : str
Check for relative, absolute or both tolerances:
- :code:`tol_mode='abs'` (default)
- :code:`tol_mode='rel'`
- :code:`tol_mode='both'`
Returns
-------
val : float
x-value of zero crossing.
Note
----
Algorithm
.. math::
x_{i+1} = x_{i} - \frac{f(x_{i})}{\frac{df}{dx}(x_{i})}\\
f(x_{i}) \leq \epsilon
"""
# default valaues
x = kwargs.get('val0', 300)
valmin = kwargs.get('valmin', 70)
valmax = kwargs.get('valmax', 3000)
max_iter = kwargs.get('max_iter', 10)
tol_rel = kwargs.get('tol_rel', ERR )
tol_abs = kwargs.get('tol_abs', ERR )
tol_mode = kwargs.get('tol_mode', 'abs')

# start newton loop
expr = True
i = 0
while expr:
# calculate function residual and new value
res = y - func(params, x)
x += res / deriv(params, x)

# check for value ranges
if x < valmin:
x = valmin
if x > valmax:
x = valmax
i += 1

if i > max_iter:
msg = ('Newton algorithm was not able to find a feasible value '
'for function ' + str(func) + '. Current value with x=' +
str(x) + ' is ' + str(func(params, x)) +
', target value is ' + str(y) + '.')
logger.debug(msg)

break
if tol_mode == 'abs' or y == 0:
expr = abs(res) >= tol_abs
elif tol_mode == 'rel':
expr = abs(res / y) >= tol_rel
else:
expr = abs(res / y) >= tol_rel or abs(res) >= tol_abs

return x


def newton_with_kwargs(
derivative, target_value, val0=300, valmin=70, valmax=3000, max_iter=10,
tol_rel=ERR, tol_abs=ERR, tol_mode="rel", **function_kwargs
tol_rel=ERR, tol_abs=ERR ** 2, tol_mode="rel", **function_kwargs
):

# start newton loop
Expand All @@ -629,6 +524,9 @@ def newton_with_kwargs(
function = function_kwargs["function"]
relax = 1

if abs(target_value) <= 1e-6:
tol_mode = "abs"

while expr:
# calculate function residual and new value
function_kwargs[parameter] = x
Expand Down
24 changes: 13 additions & 11 deletions tests/test_tools/test_helpers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,15 +9,16 @@
SPDX-License-Identifier: MIT
"""
from pytest import approx

from tespy.tools.helpers import newton
from tespy.tools.helpers import newton_with_kwargs


def func(params, x):
def func(x, **kwargs):
return x ** 2 + x - 20


def deriv(params, x):
def deriv(x, **kwargs):
return 2 * x + 1


Expand All @@ -36,24 +37,25 @@ def test_newton_bounds():
The function is x^2 + x - 20, there crossings are -5 and 4.
"""
result = newton(func, deriv, [], 0, valmin=-10, valmax=10, val0=0)
kwargs = {"function": func, "parameter": "x"}
result = newton_with_kwargs(deriv, 0, valmin=-10, valmax=10, val0=0, **kwargs)
msg = ('The newton algorithm should find the zero crossing at 4.0. ' +
str(round(result, 1)) + ' was found instead.')
assert 4.0 == result, msg
assert 4.0 == approx(result), msg

result = newton(func, deriv, [], 0, valmin=-10, valmax=10, val0=-10)
result = newton_with_kwargs(deriv, 0, valmin=-10, valmax=10, val0=-10, **kwargs)
msg = ('The newton algorithm should find the zero crossing at -5.0. ' +
str(round(result, 1)) + ' was found instead.')
assert -5.0 == result, msg
assert -5.0 == approx(result), msg

result = newton(func, deriv, [], 0, valmin=-4, valmax=-2, val0=-3)
result = newton_with_kwargs(deriv, 0, valmin=-4, valmax=-2, val0=-3, **kwargs)
msg = ('The newton algorithm should not be able to find a zero crossing. '
'The value ' + str(round(result, 1)) + ' was found, but the '
'algorithm should have found the lower boundary of -4.0.')
assert -4.0 == result, msg
assert -4.0 == approx(result), msg

result = newton(func, deriv, [], 0, valmin=-20, valmax=-10, val0=-10)
result = newton_with_kwargs(deriv, 0, valmin=-20, valmax=-10, val0=-10, **kwargs)
msg = ('The newton algorithm should not be able to find a zero crossing. '
'The value ' + str(round(result, 1)) + ' was found, but the '
'algorithm should have found the upper boundary of -10.0.')
assert -10.0 == result, msg
assert -10.0 == approx(result), msg

0 comments on commit 3fddc8d

Please sign in to comment.