Note

Last update 26/06/2020

Note

If you build your own component using SuperflexPy, you should also share your implementation with the community and contribute to the Elements list page to make other users aware of your implementation.

Expand SuperflexPy: build customized elements

In this page we will illustrate how to create customized elements using the SuperflexPy framework.

The examples include three elements:

• Linear reservoir
• Half-triangular lag function
• Parameterized splitter

The elements presented here are as simple as possible, since the focus is in explaining how the framework works, rather than providing code to use in a practical application. To understand deeper how to exploit all the functionalities of SuperflexPy, please have a look at the elements that have been already implemented (importing path superflexpy.implementation.elements).

In this page we report, for brevity, only the code, without docstring. The complete code used to generate this page is available at the path doc/source/build_element_code.py.

Linear reservoir

The linear reservoir is one of the simplest reservoir that can be built. The idea is that the output flux is a linear function of the state of the reservoir.

The element is controlled by the following differential equation

\[\frac{\textrm{d}S}{\textrm{d}t}=P-Q\]

with

\[Q=kS\]

The solution of the differential equation can be approximated using a numerical method with the equation that, in the general case, becomes:

\[\frac{S_{t+1} - S_{t}}{\Delta t}=P - Q(S)\]

Several numerical methods exist to approximate the solution of the differential equation and, usually, they differ for the state used to evaluate the fluxes: implicit Euler, for example, uses the state at the end of the time step (\(S_{t+1}\))

\[\frac{S_{t+1} - S_{t}}{\Delta t}=P - kS_{t+1}\]

explicit Euler uses the state at the beginning of the time step (\(S_t\))

\[\frac{S_{t+1} - S_{t}}{\Delta t}=P - kS_{t}\]

and so on for other methods.

Note that, even if for this simple case the differential equation can be solved analytically and the solution of the numerical approximation can be found without iteration, we will use anyway the numerical solver offered by SuperflexPy to illustrate how to proceed in a more general case where such option is not available.

The framework provides the class ODEsElement that has most of the methods required to solve the element. The class implementing the element will inherit from this and implement only a few methods.

import numba as nb
from superflexpy.framework.element import ODEsElement

class LinearReservoir(ODEsElement):

The first method to implement is the class initializer

    def __init__(self, parameters, states, approximation, id):

        ODEsElement.__init__(self,
                             parameters=parameters,
                             states=states,
                             approximation=approximation,
                             id=id)

        self._fluxes_python = [self._fluxes_function_python]  # Used by get fluxes, regardless of the architecture

        if approximation.architecture == 'numba':
            self._fluxes = [self._fluxes_function_numba]
        elif approximation.architecture == 'python':
            self._fluxes = [self._fluxes_function_python]
        else:
            message = '{}The architecture ({}) of the approximation is not correct'.format(self._error_message,
                                                                                           approximation.architecture)
            raise ValueError(message)

The main purpose of the method (lines 9-16) is to deal with the numerical solver used for solving the differential equation. In this case we can accept two architectures: pure python or numba. The option selected will change the function used to calculate the fluxes. Keep in mind that, since some operations the python implementation of the fluxes is still used, this must be always present.

The second method to define is the one that maps the (ordered) list of input fluxes to a dictionary that gives a name to these fluxes.

    def set_input(self, input):

        self.input = {'P': input[0]}

Note that the name (key) of the input flux must be the same used for the correspondent variable in the flux functions.

The third method to implement is the one that runs the model, solving the differential equation and returning the output flux.

    def get_output(self, solve=True):

        if solve:
            self._solver_states = [self._states[self._prefix_states + 'S0']]
            self._solve_differential_equation()

            self.set_states({self._prefix_states + 'S0': self.state_array[-1, 0]})

        fluxes = self._num_app.get_fluxes(fluxes=self._fluxes_python,
                                          S=self.state_array,
                                          S0=self._solver_states,
                                          **self.input,
                                          **{k[len(self._prefix_parameters):]: self._parameters[k] for k in self._parameters},
                                          )

        return [- fluxes[0][1]]

The method takes, as input, the parameter solve: if False, the state array of the reservoir will not be calculated again, potentially producing a different result, and the output will be computed based on the state that is already stored. This is the desired behavior in case of post-run inspection of the element.

Line 4 transforms the states dictionary in an ordered list, line 5 call the built-in solver of the differential equation, line 7 updates the state of the model to the one of the updated value, lines 9-14 call the external numerical solver to get the values of the fluxes (note that, for this operation, the python implementation of the fluxes is used always).

The last method(s) to implement is the one that defines the fluxes: in this case the methods are two, one for the python implementation and one for numba.

    @staticmethod
    def _fluxes_function_python(S, S0, ind, P, k):

        if ind is None:
            return (
                [
                    P,
                    - k * S,
                ],
                0.0,
                S0 + P
            )
        else:
            return (
                [
                    P[ind],
                    - k[ind] * S,
                ],
                0.0,
                S0 + P[ind]
            )

    @staticmethod
    @nb.jit('Tuple((UniTuple(f8, 2), f8, f8))(optional(f8), f8, i4, f8[:], f8[:])',
            nopython=True)
    def _fluxes_function_numba(S, S0, ind, P, k):

        return (
            (
                P[ind],
                - k[ind] * S,
            ),
            0.0,
            S0 + P[ind]
        )

They are both private static methods. Their input consists of the state used to compute the fluxes (S), initial state (S0, used to define the maximum possible state for the reservoir), index to use in the arrays (ind, all inputs are arrays and, when solving for a single time step, the index indicates the time step to look for), input fluxes (P), and parameters (k). The output is a tuple containing three elements:

• tuple with the values of the fluxes calculated according to the state; positive sign for incoming fluxes, negative for outgoing;
• lower bound for the search of the state;
• upper bound for the search of the state;

The implementation for the numba solver differs in two aspects:

• the usage of the numba decorator that defines the types of the input variables (lines 24-25)
• the fact that the method works only for a single time step and not for the vectorized solution (use python method for that)

Half-triangular lag function

The half-triangular lag function is a function that has the shape of a right triangle, growing linearly until \(t_{\textrm{lag}}\) and then zero. The growth rate (\(\alpha\)) is designed such as the total area of the triangle is one.

\[\begin{split}& f_{\textrm{lag}}=\alpha t & \quad \textrm{for }t\leq t_{\textrm{lag}}\\ & f_{\textrm{lag}}=0 & \quad \textrm{for }t>t_{\textrm{lag}}\end{split}\]

SuperflexPy provides the class LagElement that contains most of the functionalities needed to solve a lag function. The class implementing a customized lag function will inherit from it and implement only the necessary methods that are needed to calculate the transformation that needs to be applied to the incoming flux.

import numpy as np

class TriangularLag(LagElement):

The only method to implement is the private method used to calculate the weight array, that is used to distribute the incoming flux.

    def _build_weight(self, lag_time):

        weight = []

        for t in lag_time:
            array_length = np.ceil(t)
            w_i = []

            for i in range(int(array_length)):
                w_i.append(self._calculate_lag_area(i + 1, t)
                           - self._calculate_lag_area(i, t))

            weight.append(np.array(w_i))

        return weight

that makes use of a secondary private static method

    @staticmethod
    def _calculate_lag_area(bin, len):

        if bin <= 0:
            value = 0
        elif bin < len:
            value = (bin / len)**2
        else:
            value = 1

        return value

This method returns the value of the area of a triangle that is proportional to the lag function, with a smaller base bin. The _build_weight method, on the other hand, uses the output of this function to calculate the weight array.

Note that this choice of using a second static method to calculate the weight array, while being convenient, is specific to this particular case; other implementation of the _build_weight method are possible and welcome.

Parameterized splitter

A splitter is an element that takes the flux coming from one element upstream and divides it to feed multiple elements downstream. Usually, the behavior of such an element is controlled by some parameters that define the part of the flux that goes into a specific element.

While SuperflexPy can support infinite fluxes (e.g., for simulating transport processes) and an infinite number of downstream elements, the simplest case that we are presenting here has a single flux that gets split into two downstream elements. In this example, the system needs only one parameter (\(\alpha_{\textrm{split}}\)) to be defined, with the downstream fluxes that are

\[\begin{split}& Q_1^{\textrm{out}} = \alpha_{\textrm{split}} Q^{\textrm{in}} \\ & Q_2^{\textrm{out}} = \left(1-\alpha_{\textrm{split}}\right) Q^{\textrm{in}}\end{split}\]

SuperflexPy provides the class ParameterizedElement that is designed for all the generic elements that are controlled by some parameters but do not have a state. The class implementing the parameterized splitter will inherit from this and implement only some methods.

from superflexpy.framework.element import ParameterizedElement

class ParameterizedSingleFluxSplitter(ParameterizedElement):

The first thing to define are two private attributes defining how many upstream and downstream elements the splitter has; this information is used by the unit when constructing the model structure.

    _num_downstream = 2
    _num_upstream = 1

After that we need to define the function that takes the inputs and the one that calculates the outputs of the splitter.

    def set_input(self, input):

        self.input = {'Q_in': input[0]}

    def get_output(self, solve=True):

        split_par = self._parameters[self._prefix_parameters + 'split-par']

        return [
            self.input['Q_in'] * split_par,
            self.input['Q_in'] * (1 - split_par)
        ]

The two methods have the same structure of the one implemented as part of the linear reservoir example. Note that, in this case, the argument solve of the get_output method is not used but it is still required to maintain the interface consistent.