Simulation-based probabilistic-harmonic load flow for the study of DERs integration in a low-voltage distribution network

1.   Introduction

In many situations arising from applications (i.e. neuroscience, finance, reliability, ...), the quantity of interest is the first time that a random quantity crosses a given fixed level. In a mathematical framework, this corresponds to the first passage time (FPTs) of a stochastic process through an eventually time dependent boundary. However, it can happen that the FPT distribution is known as well as the random process, while one is interested in determining the corresponding time dependent boundary. This is the so-called Inverse FPT problem. This problem has been investigated both from a theoretical ^[1,2] and an empirical point of view ^[3,4,5,6] in the one-dimensional case. In ^[1] the existence and uniqueness of the solution of the Inverse FPT is studied. In ^[2] the problem is interpreted in terms of an optimal stopping problem. A numerical algorithm has been proposed in ^[6] for the Wiener process. In ^[5] the Inverse FPT of an Ornstein Uhlenbeck (OU) process has been studied and it has been applied to a classification method with applications to neuroscience. The same framework is in ^[4], where possible thresholds corresponding to Gamma distributed FPTs for an OU process has been investigated with modelling purposes.

Here, we resort again to the Inverse FPT method, we generalize the algorithms to a two-dimensional OU process and we study the possibility to have Inverse Gaussian (IG) or Gamma distributed FPTs. The choice of these two distributions grounds on their role in neurosciences ^{[7,8,9,10,11,12,13]} and reliability theory ^[14,15].

In Section 2 we introduce the two dimensional Gauss-Markov process of interest, underlying some properties that we will use to deal with the Inverse FPT algorithm. In Section 3 we introduce the Inverse FPT method for a two dimensional process but we postpone to a future work the mathematical discussion about its convergence. In Section 4 we apply the algorithm to two choices of the FPT distribution, determining the thresholds corresponding to IG or Gamma FPTs probability density function (pdf). We underline the differences between the two models and we explain how heavy or light tails influence the boundary behavior. The last section discusses the obtained results in neuroscience contest. The two compartment model of Leaky Integrate and Fire type presented in ^[9] and studied in ^[16] describes the membrane potential evolution of a neuron as a two-dimensional OU process. Hence, in Section 5 we reinterpret the Inverse FPT results in this framework.

2.   The two dimensional Ornstein Uhlenbeck process

Let us consider a stochastic process $X = \{(X_1(t), X_2(t)), t\geq 0\}$ that is solution of the following stochastic differential system

with $X(0) = 0$ and where $B$ is a one-dimensional standard Brownian motion. Here, $\alpha > 0$, $\beta$, $\mu$ and $\sigma > 0$ are constants.

To solve the stochastic differential system (2.1), we rewrite it in matrix form

where

It is an autonomous linear stochastic differential equation, in particular it is a two-dimensional Ornstein Uhlenbeck process, special case of a Gauss-Markov diffusion process ^[17]. The solution of (2.2) is

It is a Gaussian vector with mean

and variance-covariance matrix $Q(t-s)$ where,

and

Trajectories of the process are plotted in Figure 1. The different behavior of the two components is evident: the noisy behavior is prevalent in $X_2$, while the first component $X_1$ is smoother. Indeed, as shown in (2.3), the multiplying function of the random term on the first component reduces the noise effect.

We consider the first passage time of the first component of the process (2.1)

where $S(t)$ is a continuous function with $S(0)\geq X_1(0) = 0$.

Note that it is possible to rewrite (2.3) in iterative form. This version is useful for simulation purposes, in order to generate the trajectories in an exact way. Discretizing the time interval $[0, T]$ with the partition $\pi: 0 = t_0 < t_1 < \dots < t_N = T$ in $N$ subintervals of constant length $h = \frac{T}{N}$, we can express the position of the process at time $t_{k+1}$ in terms of the position of the process at time $t_k$

where the term

known as innovation, is a Gaussian vector with zero mean and variance-covariance matrix $Q(h)$.

In the following we will also need the conditioned mean of the first component

and the conditioned variance of the first component

In some instances can be useful to transfer the time dependency from the boundary shape $S(t)$ to an input $M(t)$. Mathematically it is possible to relate these two situations with a simple space transformation. Indeed, the space transformation

changes our process $X$ given by (2.1), originated in $X(0) = x_0$ in presence of a time dependent boundary $S(t)$

into a two dimensional process characterized by time dependent input $M(t)$ and constant threshold $\Sigma$

Here, the term

can be interpreted as an external input acting with different weights on the two compartments. Note that in (2.14) $S(t)$ should be interpreted as a function of time and not as a boundary of a FPT problem. Indeed, in this case the boundary of the model is constant.

The stochastic process (2.1) can be used to describe a system whose behavior depends by two components that are strictly correlated by the parameter $\beta$. The second component is driven by a random Gaussian noise and its evolution is stopped when the first component reaches a given fixed level. This model, known as two compartment model, has many interesting applications, for example in neuroscience, reliability and finance.

3.   Inverse first passage time method

The inverse FPT problem consists in searching the unknown boundary $S(t)$ given that the FPT density $f_T(t)$ is known. We work under the assumption that the boundary $S(t)$ exists, it is unique and sufficiently regular.

Let us consider a diffusion process $X = \{(X_1(t), X_2(t)), t\geq 0\}$, solution of the stochastic differential equation (2.2). The proposed method is based on the numerical approximation of the following Volterra integral equation ^[16]

where $Z(t)$ is a random variable that represents the position of the second component $X_2$ of the process when the first component $X_1$ hits the boundary at time $t$, i.e.

Let us fix a time interval $[0, \Theta]$ and a partition $\pi: 0 = t_0 < t_1 < \dots < t_N = \Theta$ in $N$ subintervals of constant length $h = \frac{\Theta}{N}$. Using Euler formula for integrals ^[18], equation (3.1) can be approximated as

$\forall i = 1 \dots N$.

Equation (3.3) represents a non linear system of $N$ equations in $N$ unknown $S^*(t_1), \dots, S^*(t_N)$ that can be solved by means of root finding iterative algorithms ^[19]. Its solution gives an approximation $S^*(t)$ of the boundary $S(t)$ in the partition points $\pi$. Note that in step $i$ the only unknown quantity is $S(t_i)$ and it is estimated using the boundary approximations $S^*(t_1), \dots, S^*(t_{i-1})$, computed in the previous steps.

The quantity

is not easily handled because it depends on the unknown time dependent boundary. In general, the computation of $\theta_{k, k}$ is not trivial but, performing a suitable limit on the considered process we can show that $\theta_{k, k} = 2$ for each value of $k$. To compute (3.4) when $k\neq i$, we use a Monte Carlo method: we simulate the process $X$ until the first component exceeds the threshold and we save the corresponding value of $Z$. At step $i$, we need to compute $\theta_{i, k}$ for $k = 1, \dots, i-1$. The presence of an expectation with respect to $Z(t_k)$ determines a difficulty for the estimation of (3.4) through Monte Carlo because we need the value of $Z(t_k) = X_2(t_k)$ at time $T = t_k$. To circumvent this problem we introduce an approximate approach as follows. At step $i$ we approximate the threshold with a piecewise linear curve with knots in the already computed boundary values. Hence, for $\tau\in [t_ {j-1}, t_j]$, $\quad j = 1, \dots, i-1$ we substitute the exact boundary with

and we simulate the process up to $t_{i-1}$ or until it reaches the threshold. To compute $\theta_{i, k}$, $k = 1, \cdots, i-1$ we use only the trajectories that crossed the approximated boundary (3.5) in a neighbourhood of $t_k$. Then, in correspondence to each of these sample paths we identify with $\{Z_k, k = 1 \dots M \}$ the sequence of values of the second component of the process $X$ (when the first component has exceeded the threshold). In this way, the Monte Carlo estimate for $\theta_{i, j}$ is

It is possible to prove that this further approximation does not seriously influence the reliability of the algorithm.

4.   Examples

In this Section we illustrate the use of the Inverse FPT method through two examples. The first situation concerns FPTs with Inverse Gaussian distribution. The second one deals with the Gamma distribution. Lastly, a comparison between boundaries and drift terms arising in the two examples is developed.

4.1.   Inverse Gaussian random variable

The IG random variable $T$ has pdf

where $\rho > 0$ is the mean and $\lambda > 0$ is the shape parameter. Mean, variance and coefficient of variation are given by

Throughout all this paper, when not differently specified, we fix the values of the two compartment model as follows: $\alpha = 0.33, \beta = 0.2$ and $\sigma = 1$. We choose the shape of the IG distribution by fixing its mean $\mathbb{E}[T] = 4$ and different values of $CV$ (see Figure 2, panel (a)). From (4.2) we see that changes of $CV$ imply changes of the shape parameter $\lambda$. Moreover, as $CV$ increases, the density becomes more peaked. The corresponding shapes of the time varying thresholds are illustrated in panels (b) where the shapes of the boundary present a maximum that tends to disappear as $CV$ grows to higher values.

When $\rho$ is finite, the IG distribution has light tails but if $\rho\rightarrow \infty$ the IG becomes

and it catches the heavy tails feature of interest for the analysis of some data ^[7,8,12]. The density (4.3) is known to be the density of the FPT of a Brownian motion with zero drift and diffusion coefficient $\nu$ through a constant boundary $b$ with the relation

Since $\mathbb{E}(T) = \infty$, it makes no sense to compute the $CV$ but, in order to compare light and heavy tails distributions, we use the same values of $\lambda$ in Figure 2 and 3. In Figure 3, different shapes of the pdfs (panel (a)) and of the corresponding boundaries (panel (b)) are shown. The heavy tails of this distribution determine a new shape for the threshold that has a decreasing maximum as $\lambda$ increases, followed by a minimum and by an increasing shape of the boundary. The maximum tends to disappear for large values of $\lambda$ and the values of the boundary are essentially positive. The growth of the boundary, for larger values of $t$, stop to allow the crossing of the samples determining the tail of the distribution. Figure 3 refers to the time interval $[0, 20]$, corresponding to a low probabilistic mass. A check for longer intervals does not change the results from a qualitative viewpoint, while higher probability masses are reached (figure not shown).

4.2.   Gamma random variable

A random variable $T$ is Gamma distributed if its pdf is

Here, $\gamma > 0$ is the rate parameter and $\kappa > 0$ is the shape parameter. Such a random variable is characterized by the following mean, variance and coefficient of variation

The shapes of Gamma pdf for different values of the parameters $\gamma$ and $\kappa$ can be seen in Figure 4, panel (a). The shapes of the Gamma pdf strongly change with the value of $CV$. We recall that $CV = 1$ corresponds to the exponential distribution. Tails of the Gamma distribution are light, decaying to zero as an exponential. The corresponding shapes of the time varying thresholds are shown in panel (b) where the values of the two compartment model are: $\alpha = 0.33$, $\beta = 0.2$ and $\sigma = 1$. Here, as $CV$ increases, the maximum of the boundary disappears and the threshold time varying threshold becomes flat or, eventually when $CV = 2$, increasing. This is the main difference of the boundary behavior with respect to the IG case.

4.3.   Comparison

Often, IG and Gamma distributions appear as output of models of the same phenomenon but, for different choices of diffusion parameters. Hence, it seems useful to compare these distributions in terms of corresponding boundaries, using the Inverse FPT method. Hence, in this subsection we compare the boundaries corresponding to IG and to Gamma distributions when mean and $CV$ are the same.

Figure 5 shows the time varying boundaries corresponding to the IG (dashed) and the Gamma-distributed (solid) interspike intervals (ISIs) with the same mean value and the same $CV$. In this example $\mathbb{E}[T] = 10$ while $CV = [0.5, 1, 1.5]$. Densities and corresponding boundaries become more and more different as we increase the $CV$ value (cf. inbox of Figure 5). The different spreading of probability mass of the two classes of distributions is reflected in different shapes of the corresponding boundaries. Since the IG density has heavier tails, the probability mass should not be consumed for short times. For this reason the boundary increases allowing crossings for large times.

To help a physical interpretation of the results in terms of input of a two compartment model, in Figure 6 we compare the behavior of the two components (2.14) of $M(t)$, when the boundary is transformed into a constant through (2.11). We ask what would be the input to both compartments if the output distribution is fixed and threshold is constant. We illustrate the cases of FPT distributed as IG (a-b), IG with heavy tails (c-d) and Gamma (e-f), respectively. Reinterpreting the time dependent boundary in terms of a modification of the drift allows to interpret our results in terms of increasing or decreasing drift. We note that a positive drift on the first component is always necessary to obtain the prescribed FPT distribution. When $CV\leq 1$, FPTs distributed as Gamma or IG imply similar input. On the contrary, when tails of IG are heavy (panel (c)) or $CV$ is large enough (panels (a) or (e)), the drift term of the first component strongly changes becoming decreasing. Interestingly the behavior of the second component $\mu_2(t)$ (panels (b), (d) and (f)) is the opposite of that of $\mu_1(t)$.

Lastly, we change the comparison criterion and we apply the Inverse FPT method varying the values of the parameter $\mu$ in (2.1). We fix the values of the model as follows: $\alpha = 0.02$, $\beta = 0.02$ and $\sigma = 0.4$. We consider examples of boundaries corresponding IG or Gamma spiking densities for $CV = 0.5$ (Figure 7) or $CV = 1$ (Figure 8).

In the figures, we also compare the boundaries (thicker line) with the mean of the first component $\mathbb{E}[X_1(t)]$. We note that boundary always intersects the function $\mathbb{E}[X_1(t)]$. Interestingly, if we fix the firing FPT and its $CV$, the intersection value is the same for different values of the parameter $\mu$. This fact can be easily understood by noting that a change in $\mu$ determines the same shift both on the mean value of $X_1(t)$ and on the boundary. However, this value changes considering different $CV$s.

5.   Application to neuroscience

We give here an example of application of the Inverse FPT method to neuroscience.

Simplest neuronal models resort to one-dimensional processes to describe the membrane potential evolution. This choice implies a strong simplification of the neuronal structure that is identified by a single point. More complex models introduce bivariate stochastic processes to discriminate the membrane potential dynamics in the dendritic or in the trigger zone ^[9]. Neurophysiological reasons suggest the existence of an interaction between the membrane potential dynamics in the two zones and, when the first component (the trigger one) attains a boundary value, the neuron releases a spike. A reasonable simplification allows to add a noisy term only to the dendritic component. The reset after the spike can include both the components or only the trigger zone. Here, we will consider only the case of total resetting of both components to a resting value that we fix, for simplicity, equal to zero.

In this framework, the stochastic process (2.1) describes the depolarization of the trigger zone and the dendritic one, respectively ^[9]. The model assumes that external inputs, with intensity $\mu$ and variability $\sigma$, influence the second compartment and a weight $\beta$ takes into account the interconnection between the parts of the neuron. Moreover, the constant $\alpha > 0$ accounts for the spontaneous membrane potential decay (cf. Figure 9). Then, the FPT $T$ mimics the ISI of the neuron and the boundary $S(t)$ corresponds to the spiking threshold for the neuron.

Often, IG and Gamma distribution fit neuronal data and FPT may help to interpret the presence of these distributions. Here, we reinterpret Figures 2-8 in the neuronal model framework. Hence, constants $\alpha$ and $\beta$ will be measured in $\text{ ms}^{-1}$, while $\mu$ will be measured in $\text{mV ms}^{-1}$ and $\sigma$ in $\text{ mV ms}^{-1/2}$.

Figures 2, 3, 4 and 5 reinterpret the heaviness of the tail of the ISI distributions in terms of the threshold shapes. As we increase the $CV$ value, IG and Gamma densities and the corresponding boundaries become more and more different. This means that $CV$ plays an important role in the formulation of the model. Moreover, in the case of the Gamma ISI distribution, the slope of the threshold becomes increasing when CV is large enough. A similar increasing behavior of the boundary could be obtained with IG distributed ISIs with heavy tails (cf. Figure 3).

In Figures 7 and 8 we investigate not only the behavior of the time varying firing threshold but also the dynamics of the underlying two compartment neuronal model. The mean membrane potentials and the corresponding boundaries are plotted for different values of the mean input $\mu$ and for different values of $CV$. For low input $\mu$, the curves exhibit a maximum after which the firing threshold starts to decrease. As the input $\mu$ increases, the firing threshold from concave and decreasing become convex and increasing. Indeed, a big input $\mu$ facilitates the spiking. Therefore, to obtain the assigned distribution, the threshold must move away, becoming increasing.

Lastly in Figure 10 we study the role of the parameter $\beta$ in the model, applying the Inverse FPT method to IG (panel (a)) and Gamma (panel (b)) FPT distribution and varying the values of the parameter $\beta$. As $\beta$ decreases, the boundary becomes almost constant and equal to zero. This is consistent with the fact that, for $\beta = 0$, the two components of the process gets independent and $X_1$ is deterministic and equal to zero, since $X(0) = 0$. Then, in order to have a crossing and to get a prescribed distribution, the threshold should approach zero.

6.   Conclusions

The extension of the Inverse FPT method to two-dimensional OU diffusion processes allows to study the shape of the boundaries for a given FPT pdf. We applied the algorithm to FPT distributed as an Inverse Gaussian and as a Gamma random variable. Differences in the boundary shape corresponding to FPTs with heavy or light tails enlighten different features of the corresponding two compartment model.

Lastly, we reinterpret the obtained results in a neuroscience framework. The shape of the boundaries corresponding to different firing distributions may enlighten features of the model eventually recognizing instances of scarce physiological significance such as diverging thresholds.

Acknowledgements

The authors are grateful to Professor L. Sacerdote and Professor L. Kostal for their interesting and useful comments and suggestions. This work is partially supported by INDAM-GNCS.

Conflict of interest

The authors declare there is no conflict of interest.