Research article Topical Sections

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


  • The integration of distributed energy resources (DERs) and, therefore, power electronic devices into distribution networks leads to harmonic distortion injection. However, studying harmonic distortion solely through deterministic approaches presents challenges due to the inherent random behavior of DERs. This study introduced a strategy that leverages PowerFactory's harmonic load flow tool. By combining it with Python co-simulation, probabilistic load flows can be developed. These load flows utilize current sources to represent harmonic distortion emitters with predefined harmonic spectra. The proposed strategy was implemented on a real network, where two different capacities of DERs were integrated at various locations within the network. The distributions for the total harmonic distortion of voltage (THDv) and the total harmonic distortion of current (THDi) were obtained 24 hours a day in nodes and lines of the network. The procedure allowed considering the uncertainty associated to the DERs integration in distribution networks in the study of harmonic distortion, which, speaking from a simulation approach, is scarce in the literature.

    Citation: Cristian Cadena-Zarate, Juan Caballero-Peña, German Osma-Pinto. Simulation-based probabilistic-harmonic load flow for the study of DERs integration in a low-voltage distribution network[J]. AIMS Electronics and Electrical Engineering, 2024, 8(1): 53-70. doi: 10.3934/electreng.2024003

    Related Papers:

    [1] Aniello Buonocore, Luigia Caputo, Enrica Pirozzi, Maria Francesca Carfora . Gauss-diffusion processes for modeling the dynamics of a couple of interacting neurons. Mathematical Biosciences and Engineering, 2014, 11(2): 189-201. doi: 10.3934/mbe.2014.11.189
    [2] Aniello Buonocore, Luigia Caputo, Enrica Pirozzi, Maria Francesca Carfora . A simple algorithm to generate firing times for leaky integrate-and-fire neuronal model. Mathematical Biosciences and Engineering, 2014, 11(1): 1-10. doi: 10.3934/mbe.2014.11.1
    [3] Giuseppe D'Onofrio, Enrica Pirozzi . Successive spike times predicted by a stochastic neuronal model with a variable input signal. Mathematical Biosciences and Engineering, 2016, 13(3): 495-507. doi: 10.3934/mbe.2016003
    [4] Virginia Giorno, Amelia G. Nobile . Exact solutions and asymptotic behaviors for the reflected Wiener, Ornstein-Uhlenbeck and Feller diffusion processes. Mathematical Biosciences and Engineering, 2023, 20(8): 13602-13637. doi: 10.3934/mbe.2023607
    [5] Virginia Giorno, Serena Spina . On the return process with refractoriness for a non-homogeneous Ornstein-Uhlenbeck neuronal model. Mathematical Biosciences and Engineering, 2014, 11(2): 285-302. doi: 10.3934/mbe.2014.11.285
    [6] Shinsuke Koyama, Ryota Kobayashi . Fluctuation scaling in neural spike trains. Mathematical Biosciences and Engineering, 2016, 13(3): 537-550. doi: 10.3934/mbe.2016006
    [7] Giacomo Ascione, Enrica Pirozzi . On a stochastic neuronal model integrating correlated inputs. Mathematical Biosciences and Engineering, 2019, 16(5): 5206-5225. doi: 10.3934/mbe.2019260
    [8] Asma Alshehri, John Ford, Rachel Leander . The impact of maturation time distributions on the structure and growth of cellular populations. Mathematical Biosciences and Engineering, 2020, 17(2): 1855-1888. doi: 10.3934/mbe.2020098
    [9] Jakub Cupera . Diffusion approximation of neuronal models revisited. Mathematical Biosciences and Engineering, 2014, 11(1): 11-25. doi: 10.3934/mbe.2014.11.11
    [10] Massimiliano Tamborrino . Approximation of the first passage time density of a Wiener process to an exponentially decaying boundary by two-piecewise linear threshold. Application to neuronal spiking activity. Mathematical Biosciences and Engineering, 2016, 13(3): 613-629. doi: 10.3934/mbe.2016011
  • The integration of distributed energy resources (DERs) and, therefore, power electronic devices into distribution networks leads to harmonic distortion injection. However, studying harmonic distortion solely through deterministic approaches presents challenges due to the inherent random behavior of DERs. This study introduced a strategy that leverages PowerFactory's harmonic load flow tool. By combining it with Python co-simulation, probabilistic load flows can be developed. These load flows utilize current sources to represent harmonic distortion emitters with predefined harmonic spectra. The proposed strategy was implemented on a real network, where two different capacities of DERs were integrated at various locations within the network. The distributions for the total harmonic distortion of voltage (THDv) and the total harmonic distortion of current (THDi) were obtained 24 hours a day in nodes and lines of the network. The procedure allowed considering the uncertainty associated to the DERs integration in distribution networks in the study of harmonic distortion, which, speaking from a simulation approach, is scarce in the literature.



    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.

    Let us consider a stochastic process X={(X1(t),X2(t)),t0} that is solution of the following stochastic differential system

    {dX1(t)={αX1(t)+β[X2(t)X1(t)]}dtdX2(t)={αX2(t)+β[X1(t)X2(t)]+μ}dt+σdBt (2.1)

    with X(0)=0 and where B is a one-dimensional standard Brownian motion. Here, α>0, β, μ and σ>0 are constants.

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

    dX(t)=[AX(t)+M(t)]dt+GdB(t), (2.2)

    where

    A=(αβββαβ)M(t)=M=(0μ) and G=(000σ).

    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

    {X1(t)=μ2(1eαtα1e(α+2β)tα+2β)+σ2t0(eα(ts)e(α+2β)(ts))dB(s)X2(t)=μ2(1eαtα+1e(α+2β)tα+2β)+σ2t0(eα(ts)+e(α+2β)(ts))dB(s). (2.3)

    It is a Gaussian vector with mean

    m(t)=E(X(t))=[μ2(1eαtα1e(α+2β)tα+2β)μ2(1eαtα+1e(α+2β)tα+2β)] (2.4)

    and variance-covariance matrix Q(ts) where,

    Q(t)=[Q(11)Q(12)Q(12)Q(22)](t) (2.5)

    and

    Q(11)(t)=12(1α2α+β+1α+2βe2αt(1α2e2βta+β+e4βtα+2β))Q(12)(t)=121e2αtα1e2(α+2β)tα+2βQ(22)(t)=12(1α+2α+β+1α+2βe2αt(1α+2e2βtα+β+e4βtα+2β)).

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

    Figure 1.  Three sample paths of the two components of the process X. The parameters of the two compartment model are α=0.33, β=0.2, μ=0 and σ=1.

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

    T=inf{t>0:X1(t)>S(t)} (2.6)

    where S(t) is a continuous function with S(0)X1(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 π:0=t0<t1<<tN=T in N subintervals of constant length h=TN, we can express the position of the process at time tk+1 in terms of the position of the process at time tk

    X(tk+1)=12[eαh+e(α+2β)heαhe(α+2β)heαhe(α+2β)heαh+e(α+2β)h]X(tk)+μ2[1eαhα1e(α+2β)hα+2β1eαhα+1e(α+2β)hα+2β]+σ2Ik (2.7)

    where the term

    Ik=[tk+1tk(eα(tt+1s)e(α+2β)(tk+1s))dB(s)tk+1tk(eα(tt+1s)e(α+2β)(tk+1s))dB(s)], (2.8)

    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

    m(1)(t|(X1(θ),X2(θ)),θ)=E(X1(t)|X(θ)=(X1(θ),X2(θ)))=μ2(2βαeα(tθ)2βeα(tθ)+αe(α+2β)(tθ)α(α+2β))X1(θ)eα(tθ)2(1+e2β(tθ))X2(θ)eα(tθ)2(1e2β(tθ)) (2.9)

    and the conditioned variance of the first component

    Q(11)(t|(X1(θ),X2(θ)),θ)=Var(X1(t)|X(θ)=(X1(θ),X2(θ)))=σ2e2α(tθ)[2αe2β(tθ)(α+2β)αe4β(tθ)(α+β)+2β2e2α(tθ)α23αβ2β2]8α(α+β)(α+2β) (2.10)

    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

    Y1(t)=X1(t)S(t)+Σ. (2.11)

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

    {dX1(t)={αX1(t)+β[X2(t)X1(t)]}dtdX2(t)={αX2(t)+β[X1(t)X2(t)]+μ}dt+σdBtX(0)=x0S(t) (2.12)

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

    {dY1(t)={αY1(t)+β[X2(t)Y1(t)]+μ1(t)}dtdX2(t)={αX2(t)+β[Y1(t)X2(t)]+μ2(t)}dt+σdBtX(0)=x0S(0)+ΣΣ. (2.13)

    Here, the term

    M(t)=[μ1(t)μ2(t)]=[(α+β)(S(t)Σ)S(t)μ+β(S(t)Σ)] (2.14)

    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 β. 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.

    The inverse FPT problem consists in searching the unknown boundary S(t) given that the FPT density fT(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={(X1(t),X2(t)),t0}, solution of the stochastic differential equation (2.2). The proposed method is based on the numerical approximation of the following Volterra integral equation [16]

    1Erf(S(t)m(1)(t)2Q(11)(t))=t0fT(θ)EZ(θ)[1Erf(S(t)m(1)(t|(S(θ),X2(θ)),θ)2Q(11)(t|(S(θ),X2(θ)),θ))]dθ (3.1)

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

    P(Z(t)<z)=P(X2(T)<z|T=t,X(t0)=y). (3.2)

    Let us fix a time interval [0,Θ] and a partition π:0=t0<t1<<tN=Θ in N subintervals of constant length h=ΘN. Using Euler formula for integrals [18], equation (3.1) can be approximated as

    1Erf(S(ti)m(1)(ti)2Q(11)(ti))=hij=1fT(tj)EZ(tj)[1Erf(S(ti)m(1)(ti|(S(tj),X2(tj)),tj)2Q(11)(ti|(S(tj),X2(tj)),tj))] (3.3)

    i=1N.

    Equation (3.3) represents a non linear system of N equations in N unknown S(t1),,S(tN) 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 π. Note that in step i the only unknown quantity is S(ti) and it is estimated using the boundary approximations S(t1),,S(ti1), computed in the previous steps.

    The quantity

    θi,k=EZ(tk)[1Erf(S(ti)m(1)(ti|S(tk|(S(tk),X2(tk)),tk)2Q(11)(ti|(S(tk),X2(tk)),tk))] (3.4)

    is not easily handled because it depends on the unknown time dependent boundary. In general, the computation of θk,k is not trivial but, performing a suitable limit on the considered process we can show that θk,k=2 for each value of k. To compute (3.4) when ki, 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 θi,k for k=1,,i1. The presence of an expectation with respect to Z(tk) determines a difficulty for the estimation of (3.4) through Monte Carlo because we need the value of Z(tk)=X2(tk) at time T=tk. 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 τ[tj1,tj], j=1,,i1 we substitute the exact boundary with

    ˆS(τ)=S(tj)S(tj1)tjtj1τ+tjS(tj1)S(tj)tj1tjtj1 (3.5)

    and we simulate the process up to ti1 or until it reaches the threshold. To compute θi,k, k=1,,i1 we use only the trajectories that crossed the approximated boundary (3.5) in a neighbourhood of tk. Then, in correspondence to each of these sample paths we identify with {Zk,k=1M} 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 θi,j is

    ˜θi,j=1Mk=1Erf(S(ti)m(1)(ti|(S(tj),Zk),tj)2Q(11)(ti|(S(tj),Zk),tj))M.

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

    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.

    The IG random variable T has pdf

    fT(t)=[λ2πt3]1/2exp[λ(tρ)22ρ2t],t0, (4.1)

    where ρ>0 is the mean and λ>0 is the shape parameter. Mean, variance and coefficient of variation are given by

    E(T)=ρVar(T)=ρ3λCV=CV(T)=Var(T)E(T)=ρλ. (4.2)

    Throughout all this paper, when not differently specified, we fix the values of the two compartment model as follows: α=0.33,β=0.2 and σ=1. We choose the shape of the IG distribution by fixing its mean 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 λ. 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.

    Figure 2.  Probability densities (a) and evaluated boundaries (b) in the case of Inverse Gaussian-distributed FPTs with E(T)=4. Different lines correspond to different shapes of the Inverse Gaussian densities, CV=0.5 (red), CV=0.75 (magenta), CV=1 (black), CV=1.5 (blue), CV=2 (green). The parameters of the two compartment model are α=0.33, β=0.2, μ=0 and σ=1.

    When ρ is finite, the IG distribution has light tails but if ρ the IG becomes

    fT(t)=[λ2πt3]1/2exp[λ2t],t0, (4.3)

    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 ν through a constant boundary b with the relation

    λ=b2ν2. (4.4)

    Since E(T)=, it makes no sense to compute the CV but, in order to compare light and heavy tails distributions, we use the same values of λ 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 λ increases, followed by a minimum and by an increasing shape of the boundary. The maximum tends to disappear for large values of λ 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).

    Figure 3.  Probability densities (a) and evaluated boundaries (b) in the case of Inverse Gaussian-distributed FPTs with heavy tails. Different lines correspond to different values of the parameter: λ=16 (red), λ=7.11 (magenta), λ=4 (black). The parameters of the two compartment model are α=0.33, β=0.2, μ=0 and σ=1. Probability mass = [0.37 0.55 0.65] for λ=[16 7.11 4].

    A random variable T is Gamma distributed if its pdf is

    fT(t)=γκΓ(κ)tκ1eγt,t0. (4.5)

    Here, γ>0 is the rate parameter and κ>0 is the shape parameter. Such a random variable is characterized by the following mean, variance and coefficient of variation

    E(T)=κγVar(T)=κγ2CV=CV(T)=1κ. (4.6)

    The shapes of Gamma pdf for different values of the parameters γ and κ 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: α=0.33, β=0.2 and σ=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.

    Figure 4.  Probability densities (a) and evaluated boundaries (b) in the case of Gamma-distributed FPTs with E(T)=4. Different lines correspond to different shapes of the Gamma densities, CV=0.5 (red), CV=0.75 (magenta), CV=1 (black), CV=1.5 (blue), CV=2 (green). The parameters of the two compartment model are α=0.33, β=0.2, μ=0 and σ=1.

    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 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.

    Figure 5.  Comparison between the time varying boundaries corresponding to the Inverse Gaussian (dashed) and the Gamma-distributed (solid) FPTs with E(T)=10 and with varying CV, the probability densities are on the sub-plots. The parameters of the two compartment model are α=0.02, β=0.02, μ=0 and σ=0.4.

    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 CV1, 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 μ2(t) (panels (b), (d) and (f)) is the opposite of that of μ1(t).

    Figure 6.  Comparison between the values of the two components of M(t) (2.14) when the boundary is constant Σ=4 in the case of FPT distributed as IG (a-b), IG with heavy tails (c-d) and Gamma (e-f), respectively. The parameters and the colors of the functions are the same of the examples of Figures 2, 3 and 4. The parameters of the two compartment model are α=0.33, β=0.2, μ=0 and σ=1.

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

    Figure 7.  Inverse-Gaussian (a) and Gamma (c) distributed FPTs with mean FPT equal to 10 and CV=0.5. Corresponding boundaries and mean value of the first component (b, d). Different lines correspond to different values of the mean input: μ=0 (black), μ=0.3 (blue), μ=0.6 (cyan). Thicker lines correspond to the time varying boundaries. The vertical dotted lines give the FPT distribution quantiles. The parameters of the two compartment model are α=0.02, β=0.02, σ=0.4, while μ varies as specified above.
    Figure 8.  Inverse-Gaussian (a) and Gamma (c) distributed FPTs with mean FPT equal to 10 and CV=1. Corresponding boundaries and mean value of the first component (b, d). Different lines correspond to different values of the mean input: μ=0 (black), μ=0.3 (blue), μ=0.6 (cyan). Thicker lines correspond to the time varying boundaries. The vertical dotted lines give the FPT distribution quantiles. The parameters of the two compartment model are α=0.02, β=0.02, σ=0.4, while μ varies as specified above.

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

    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 μ and variability σ, influence the second compartment and a weight β takes into account the interconnection between the parts of the neuron. Moreover, the constant α>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.

    Figure 9.  Schematic representation of the two-compartment model.

    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 α and β will be measured in  ms1, while μ will be measured in mV ms1 and σ in  mV ms1/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 μ and for different values of CV. For low input μ, the curves exhibit a maximum after which the firing threshold starts to decrease. As the input μ increases, the firing threshold from concave and decreasing become convex and increasing. Indeed, a big input μ 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 β 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 β. As β decreases, the boundary becomes almost constant and equal to zero. This is consistent with the fact that, for β=0, the two components of the process gets independent and X1 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.

    Figure 10.  Boundaries corresponding to Inverse Gaussian (a) and Gamma (b) distributed FPTs with E(T)=10 and CV=1. Different lines correspond to different values of β: β=0.01 (green), β=0.02 (red), β=0.1 (magenta), β=0.5 (blue).The parameters of the two compartment model are α=0.02, μ=0, σ=0.4, while β varies as specified above.

    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.

    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.

    The authors declare there is no conflict of interest.



    [1] IRENA (2019) Innovation landscape for a renewable-powered future: Solutions to integrate variable renewables. International Renewable Energy Agency, Abu Dhabi. Available from: https://www.irena.org/publications/2019/Feb/Innovation-landscape-for-a-renewable-powered-future.
    [2] IRENA (2020) Energy Transition. Available from: https://www.irena.org/energytransition.
    [3] Caballero-Peña J, Cadena-Zarate C, Osma-Pinto G (2023) Hourly characterization of the integration of DER in a network from deterministic and probabilistic approaches using Co-simulation PowerFactory-Python. Alex Eng J 63: 283–305. https://doi.org/10.1016/j.aej.2022.08.005. doi: 10.1016/j.aej.2022.08.005
    [4] Zhi H, Zhang M, Li R, Zhao J, Wang J, Li X, et al. (2021) A Power-based Piecewise Probabilistic Harmonic Model of Secondary Residential System. 2021 IEEE 5th Conference on Energy Internet and Energy System Integration (EI2), 1175–1179. https://doi.org/10.1109/EI252483.2021.9713235 doi: 10.1109/EI252483.2021.9713235
    [5] Galvani S, Hagh MT, Sharifian MBB, Mohammadi-Ivatloo B (2019) Multiobjective Predictability-Based Optimal Placement and Parameters Setting of UPFC in Wind Power Included Power Systems. IEEE T Ind Inform 15: 878–888. https://doi.org/10.1109/TII.2018.2818821 doi: 10.1109/TII.2018.2818821
    [6] Xie X, Sun Y (2022) A piecewise probabilistic harmonic power flow approach in unbalanced residential distribution systems. International Journal of Electrical Power and Energy Systems 141: 108114. https://doi.org/10.1016/j.ijepes.2022.108114 doi: 10.1016/j.ijepes.2022.108114
    [7] Nasrfard-Jahromi F, Mohammadi M (2016) Probabilistic harmonic load flow using an improved kernel density estimator. International Journal of Electrical Power and Energy Systems 78: 292–298. https://doi.org/10.1016/j.ijepes.2015.11.076 doi: 10.1016/j.ijepes.2015.11.076
    [8] Mohammadi M. (2015) Probabilistic harmonic load flow using fast point estimate method. IET Generation, Transmission and Distribution 9: 1790–1799. https://doi.org/10.1049/iet-gtd.2014.0669 doi: 10.1049/iet-gtd.2014.0669
    [9] Palomino A, Parvania M (2020) Data-Driven Risk Analysis of Joint Electric Vehicle and Solar Operation in Distribution Networks. IEEE Open Access Journal of Power and Energy 7: 141–150. https://doi.org/10.1109/oajpe.2020.2984696 doi: 10.1109/oajpe.2020.2984696
    [10] Yan J, Zhang J, Liu Y, Lv G, Han S, Alfonzo IEG (2020) EV charging load simulation and forecasting considering traffic jam and weather to support the integration of renewables and EVs. Renew Energ 159: 623–641. https://doi.org/10.1016/j.renene.2020.03.175 doi: 10.1016/j.renene.2020.03.175
    [11] Hemmatpour MH, Rezaeian Koochi MH, Dehghanian P, Dehghanian P (2022) Voltage and energy control in distribution systems in the presence of flexible loads considering coordinated charging of electric vehicles. Energy 239: 121880. https://doi.org/10.1016/j.energy.2021.121880 doi: 10.1016/j.energy.2021.121880
    [12] Alzahrani A, Alharthi H, Khalid M (2020) Minimization of Power Losses through Optimal Battery Placement in a Distributed Network with High Penetration of Photovoltaics. Energies 13: Article Number 140. https://doi.org/10.3390/en13010140 doi: 10.3390/en13010140
    [13] Schneider KP, Mather BA, Pal BC, Ten CW, Shirek GJ, Zhu H, et al. (2018) Analytic Considerations and Design Basis for the IEEE Distribution Test Feeders. IEEE T Power Syst 33: 3181–3188. https://doi.org/10.1109/TPWRS.2017.2760011 doi: 10.1109/TPWRS.2017.2760011
    [14] Galvani S, Rezaeian Marjani S, Morsali J, Ahmadi Jirdehi M (2019) A new approach for probabilistic harmonic load flow in distribution systems based on data clustering. Electr Power Syst Res 176: 105977. https://doi.org/10.1016/j.epsr.2019.105977 doi: 10.1016/j.epsr.2019.105977
    [15] Xie X, Peng F, Zhang Y (2022) A data-driven probabilistic harmonic power flow approach in power distribution systems with PV generations. Appl Energy 321: 119331. https://doi.org/10.1016/j.apenergy.2022.119331 doi: 10.1016/j.apenergy.2022.119331
    [16] Ma Y, Azuatalam D, Power T, Chapman AC, Verbič G (2019) A novel probabilistic framework to study the impact of photovoltaic-battery systems on low-voltage distribution networks. Appl Energy 254: 113669. https://doi.org/10.1016/j.apenergy.2019.113669 doi: 10.1016/j.apenergy.2019.113669
    [17] Lehtonen M (1998) A method for probabilistic harmonic load-flow analysis in power systems. Eur T Electr Power 8: 47–50. https://doi.org/10.1002/etep.4450080108 doi: 10.1002/etep.4450080108
    [18] Russo A, Varilone P, Caramia P (2014) Point estimate schemes for probabilistic harmonic power flow. 2014 16th International Conference on Harmonics and Quality of Power (ICHQP), 19–23. https://doi.org/10.1109/ICHQP.2014.6842890 doi: 10.1109/ICHQP.2014.6842890
    [19] Yu G, Lin T (2016) 2m+1 point estimate method for probabilistic harmonic power flow. 2016 IEEE Power and Energy Society General Meeting (PESGM), 1–5. https://doi.org/10.1109/PESGM.2016.7741657 doi: 10.1109/PESGM.2016.7741657
    [20] Ismael SM, Abdel Aleem SHE, Abdelaziz AY, Zobaa AF (2019) Probabilistic Hosting Capacity Enhancement in Non-Sinusoidal Power Distribution Systems Using a Hybrid PSOGSA Optimization Algorithm. Energies 12: Article Number 1018. https://doi.org/10.3390/en12061018 doi: 10.3390/en12061018
    [21] Gandoman FH, Sharaf AM, Abdel Aleem SHE, Jurado F (2017) Distributed FACTS stabilization scheme for efficient utilization of distributed wind energy systems. Int T Electr Energy 27: e2391. https://doi.org/10.1002/etep.2391 doi: 10.1002/etep.2391
    [22] Martínez-Peñaloza A, Osma-Pinto G (2021) Analysis of the performance of the Norton equivalent model of a photovoltaic system under different operating scenarios. Int Rev Electr Eng 16: 328–343. https://doi.org/10.15866/iree.v16i4.20278 doi: 10.15866/iree.v16i4.20278
    [23] Caro L, Ramos G, Montenegro D, Celeita D (2020) Variable Harmonic Distortion in Electric Vehicle Charging Stations. 2020 IEEE Industry Applications Society Annual Meeting, 1–6. https://doi.org/10.1109/IAS44978.2020.9334798 doi: 10.1109/IAS44978.2020.9334798
    [24] Angelim JH, Affonso CD (2019) Probabilistic Impact Assessment of Electric Vehicles Charging on Low Voltage Distribution Systems. 2019 IEEE PES Innovative Smart Grid Technologies Conference - Latin America (ISGT Latin America), 1–6. https://doi.org/10.1109/ISGT-LA.2019.8895494 doi: 10.1109/ISGT-LA.2019.8895494
    [25] Prusty BR, Jena D (2017) A critical review on probabilistic load flow studies in uncertainty constrained power systems with photovoltaic generation and a new approach. Renew Sust Energy Rev 69: 1286–1302. https://doi.org/10.1016/j.rser.2016.12.044 doi: 10.1016/j.rser.2016.12.044
    [26] Jahic A, Eskander M, Schulz D (2019) Bus depot simulator: Steady-state python and digsilent co-simulation for large-scale electric bus depots. NEIS 2019 - Conference on Sustainable Energy Supply and Energy Storage Systems, 85–91.
  • This article has been cited by:

    1. Alexander Klump, The Inverse First-passage Time Problem as Hydrodynamic Limit of a Particle System, 2023, 25, 1387-5841, 10.1007/s11009-023-10020-7
    2. Alexander Klump, Mladen Savov, Conditions for existence and uniqueness of the inverse first-passage time problem applicable for Lévy processes and diffusions, 2025, 35, 1050-5164, 10.1214/25-AAP2157
  • Reader Comments
  • © 2024 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1524) PDF downloads(135) Cited by(1)

Figures and Tables

Figures(9)  /  Tables(1)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog