Nonlinear Systems of Volterra Equations with Piecewise Smooth Kernels: Numerical Solution and Application for Power Systems Operation

: The evolutionary integral dynamical models of storage systems are addressed. Such models are based on systems of weakly regular nonlinear Volterra integral equations with piecewise smooth kernels. These equations can have non-unique solutions that depend on free parameters. The objective of this paper was two-fold. First, the iterative numerical method based on the modiﬁed Newton–Kantorovich iterative process is proposed for a solution of the nonlinear systems of such weakly regular Volterra equations. Second, the proposed numerical method was tested both on synthetic examples and real world problems related to the dynamic analysis of microgrids with energy storage systems.


Introduction
The conventional models of dynamical systems in various fields of physics, engineering, energy, economics and other fields are based on various nonstandard classes of functional equations including differential, integral and integral-differential with various delay types. Such mathematical tools are able to properly describe the process of interest and it is the only way to do it properly in most of the cases. Integral equations with delays model the various important classes of the dynamical processes. One of the most interesting classes of functional equations with delays is the evolutionary Volterra integral equations (VIE) [1,2].
The Volterra equation is of course the classical problem, which has been intensively studied during the last century. However, only a few authors have considered these equations with discontinuous kernels because of the fundamental difficulty of such equations. The integral equations of the second kind with discontinuous kernel were considered by Evans [3]. The uniform approximation theory for integral equations with discontinuous kernels was constructed by Anselone [4]. The sufficient conditions of solution existence of the linear VIE with jump discontinuity on the single curve were derived by Denisov and Lorenzi [5]. The papers [6][7][8] are devoted to a numerical treatment of integral dynamical systems described, using the nonlinear Volterra integral equations with unknown functions in the lower limits of integration. This type of equations are also connected to the optimal control problems in macro-economical Vintage Capital Management (VCM), which was discussed in a series of publications by Hritonenko and Yatsenko [9]. It is to be noted that the Volterra models take into account the memory effects of the system when the past states influence the present. Such evolutionary integral dynamic models are widely used in various applications including economics, ecology and population dynamics studies, see monograph [9], papers [10][11][12] and references therein. We refer readers to monograph [13] as an excellent introduction to the theory of VIE of the first kind with variable limits of integration.
It is to be outlined that all these models are related to the theory of VIE of the 1st kind with discontinuous kernels along the finite number of continuous curves recently introduced by Sidorov [14,15]. The existence of a continuous solution depending on free parameters and sufficient conditions for the existence of a unique continuous solution of the systems VIE of the 1st kind with discontinuous kernels were derived in [16]. The class of Volterra operator equations of the first kind with piecewise continuous kernels is studied in [17]. The problem of the generalized solution (in the Sobolev-Schwartz sense) to the Volterra equations with piecewise continuous kernels was considered in [18]. Convergence of the successive approximation method for construction of the main solutions in the sense of Kantorovich of nonlinear Volterra equations is proved in [19]. For systematic studies of VIE of the 1st kind with discontinuous kernels readers may refer to book [20] and part 1 in [21].
One of the prominent fields of applications of the VIE is the mathematical modeling of power systems [22]. The Volterra model can effectively generalize conventional discrete integral model and take into account both state of health of the storage and the availability of generation/storage of the microgrid as shown in [23]. The applicability of the linear system of VIEs to model the community of microgrids (MGs) with storage system is first demonstrated in [24].
The linear Volterra models are not able to handle the various important features of the contemporary storage systems such as the nonlinear dependence of the storage efficiency on state of charge or the complex processes of storage degradation. The objective of the present research was to report on the iterative numerical method based on the modified Newton-Kantorovich procedure applied to the system of nonlinear Volterra integral equations. The contribution of this paper is an iterative numerical method (equipped with the convergence theorem), which can use both the piecewise constant approximation and polynomial collocation method.
The efficiency of proposed methods are based both on toy examples with synthetic data and a real dataset from the dynamic models of a AC/DC hybrid isolated power system consisting of four power grids with renewable generation units and energy storage systems.
The remainder of this paper is structured as follows. In Section 1, the problem statement is formulated. Section 2 is dedicated to the modified Newton-Kantorovich iterative process' application to the nonlinear systems of VIEs. Section 3 deals with the discretization process of the system of linear integral equation appearing on each step of proposed iterative process in Section 2. The constructed theory and numerical method are illustrated on model examples in Section 4. Storage system analysis in the microgrid, using the system of the Volterra equation, is fulfilled in Section 5. The concluding remarks finalize the paper.

Problem Statement
Let us consider the following nonlinear system of the Volterra integral equations of the first kind where t ∈ [0, T] and kernels h i (t, s, x 1 (s), x 2 (s), · · · , x n (s)) are jump discontinuous on the curves α j (t), j = 1, 2, . . . , n − 1, and can be defined as Here Let us recall the existence and uniqueness theorem for the scalar case (i = 1 in (1)) from [25].
In the present paper, we propose an iterative method for solving such nonlinear systems. Our method employs the linearization of the integral vector-operator. To solve the systems of linear equations with discontinuous kernels arising in the iterative process, we suggest two approaches: a generalization of the direct discretization algorithm proposed in [25] and the polynomial-collocation scheme. The first approach is based on a piecewise constant approximation of exact solutions and has the first order of accuracy. The second approach uses the polynomial approximation of exact solutions and allows us to derive more accurate solutions. The proposed Volterra models and numerical methods are applied to perform the dynamic analysis of microgrids with energy storage following [24].

Description of the Iterative Approximate Method
Let P i denotes the left hand side of i-th equation of the system (1): and X denotes the desired vector of functions, and F denotes the vector of source functions f i We also introduce a matrix integral operator whose components are nonlinear integral operators Let us rewrite system (1) in the operator form as follows To construct an iterative numerical method for solving Equation (5), the following modified Newton-Kantorovich scheme [26] is employed where is the vector in the initial approximation function, and the solution of (5) is defined as the limit X * = lim m→∞ X m . Derivative P (X 0 ) of operator (4) is as follows Here the components ∂P ij ∂x j (t) (x 0 j ), i, j = 1, 2, . . . , n, are defined as follows The limit transition under the sign of the integral implies: On iterations with the number m we have an operator equation with respect to the correction which in expanded form is equivalent to the following system of equations The latter system, after a number of simplifications, is brought to mind where System (8) is the system of linear Volterra integral equations with discontinuous kernels K ij (t, s)G ij (s, x 0 j (s)), i = 1, 2, . . . , n, j = 1, 2, . . . , n, which do not change from one iteration to another. A numerical method for solving similar systems of equations involved in iterative process (6) is proposed below.

The Convergence Theorem
Let C [0,T] be a Banach vector-space of continuous vector-functions The following theorem of convergence (based on the general theory proposed in classical monograph [26], see Theorem 6, page 532) for the iterative process (6) takes place.

Theorem 2.
Let the operator P has a continuous second derivative in the sphere Ω 0 ( X − X 0 r) and the following conditions hold: has a unique solution X * in Ω 0 , process (6) converges to X * , and the velocity of convergence is estimated by the inequality

Proof.
In order to prove this theorem, one must show that Equation (8) is uniquely solvable (including the case m = 0), i.e., condition 1 of the theorem holds. The solvability of linear systems of such Volterra integral equations was discussed in monograph [20], chapter 3. Condition 2 follows from condition 1. Indeed, having obtained a bounded solution of a linear system from condition, one can estimate it and obtain a constant η. A further argumentation is performed under conditions provided the uniqueness of solution for such systems (see [20] for details).
Let us also verify the boundedness of the second derivative P (X 0 ) (X,X) for estimating the constant L in condition 3. The second derivative of the operator P is a bilinear integral operator having the following form with the components defined as follows Therefore, assuming additionally that the functions G ij (s, x(s)) have bounded second derivatives with respect to x in the neighborhood of initial approximation X 0 , we can conclude that the operator P (X 0 ) (X,X) is bounded.

Problem Formulation
In this paragraph, we propose a numerical method for solving systems of linear integral equations arising at each step of the iterative process (6).
Let us consider the following system of linear integral equations of the first kind functions h i (t, s), i = 1, 2, . . . , n, undergo a finite discontinuity on the curves α j (t), j = 1, 2, . . . , n − 1 and Here The new kernels K ij (t, s) and right-hand side f i (t) are continuous and sufficiently smooth functions introduced to simplify the exposition.

Piecewise Constant Approximation
In order to construct a numerical solution of (9) on the segment [0, T] (under the conditions of existence of a single continuous solution), we introduce a grid of nodes (not necessarily uniform) An approximate solution of the system (9) will be sought in the form of a vector of piecewise constant functions with uncertain coefficients so far x j i , i = 1, 2, . . . , n, j = 1, 2, . . . , N. To determine the values x 0 i = x i (0), i = 1, 2, . . . , n, we differentiate both parts of each Equation (9) Using the latter ratio, we come to the necessity of solving a system of equations with respect to values x 0 j , j = 1, 2, . . . , n. Here it is assumed that the components of the Equation (9) are such that the system of linear algebraic Equation (10) has a unique nontrivial solution.
Next, we introduce the designation f k i = f i (t k ), I = 1, . . . , n, k = 1, . . . , N. To determine the coefficients x 1 j , write down each equation of the system (9) at t = t 1 .
Since at this step, the lengths of all integration segments α j (t 1 ) − α j−1 (t 1 ) in (17) do not exceed the maximum grid step h, and the components of the approximate solution take the values x 1 j , j = 1, 2, . . . , n, then applying the quadrature formula of the midpoint rectangles, we have a system of , of which shall be determined by the value of x 1 j , j = 1, 2, . . . , n. Let l K j further denote the number of the grid's segment, which gets the value α j (t k ), that is, the condition t l k j −1 ≤ α j (t k ) ≤ t l k j is satisfied.
Assume now that the values x 2 j , x 3 j , · · · , x l K j −1 j , K = 1, 2, . . . , n are already found. Let us transform system (9) and require the following equalities satisfaction at the point The given approximation leads to the following system For the integrals calculation in (12), the compound midpoint rectangle quadrature constructed on the auxiliary grid of nodes bound at each concrete value N to the discontinuity curve α i (t) of kernel K(t, s) is used. It is not difficult to see that such an approximation method has the first order of accuracy and the following error estimate takes place

Remark 1.
It should be noted that the proposed numerical approach also allows for a more accurate approximation of the solution. In particular, when using piecewise linear approximation, the order of accuracy increases by one. However, such discretized systems may have non-unique solutions (due to the discontinuities α i (t), i = 1, 2, . . . , n, they can be under-determined or have many solutions). Taking into account this fact we suggest a new alternative numerical method below.
Let us rewrite the formal notation (13) as follows

Collocation
For construction of numerical solution of system (14) on [0, T] (in conditions of unique continuous solution existence) let us introduce the mesh (not necessarily uniform in general) As an approximate solution of the system (14), we search in the following form with unknown coefficients A ij , i = 1, 2, . . . , n, j = 0, 2, . . . , m.
To define the coefficients A i0 = x i (0), i = 1, 2, . . . , n, we differentiate Equation (14) with respect to t Using the latter equation we obtain the following system with respect to A j0 , j = 1, 2, . . . , n. Here we assume that components of Equation (14) such as system of linear algebraic Equation (16) enjoys unique nontrivial solution.
To define the coefficients A ij each equation of the system (14) we write in the point t = t k .
Here, for integrals calculation in (20), (21) one must employ the compound midpoint quadrature rule constructed on the auxiliary mesh linked to the lines α i (t) of the kernels h i (t, s), i = 1, 2, . . . , n discontinuities.
Finally, as solution to the system (24), the coefficients A ij are be determined and the polynomial approximation of the components of the approximate solution X m (t) is derived.

System of Linear Integral Equations
In order to illustrate the convergence of proposed algorithm in Section 4.3.2, we consider two system of equations. where The exact solution of the system (28) is x 1 (t) = cos(t), x 2 (t) = sin(t). The numerical solution of the system (28) is given in Table 1, where m is order of approximating polynomial,  Let us consider the following system of equations where f i (t), i = 1, 2, 3, are selected, the exact solution is as follows In Table 2, the numerical solution of system (29) is presented. Therefore, we suggested the efficient method for solution of systems (13) for relatively small T. The polynomial spline-approximation can be used to find the solution larger interval.

Nonlinear Equations
Let us first consider the following scalar VIE: where Here x * (t) = t 2 , t ∈ [0, 1] is the exact solution to Equation (30).
Results concerning solution of this equation, using Newton-Kantorovich method (combined with direct discretization of linear equations described in Section 4.2), are shown in Table 3, where the following notations were used: h is the mesh step, m is iterations number of Newton-Kantorovich method, ε = max

Nonlinear Systems of Equations
Let us consider the following system of nonlinear integral equations where t ∈ [0, 1], and functions f i (t), i = 1, 2, are selected such as is solution of the system. Let us select the initial approximation far enough from the exact one Results are given in Table 4, where N it is iterations number, m is order of approximating polynomials. Let us now consider the same system, but not the exponential exact solution where t ∈ [0, 1] and functions f i (t), i = 1, 2, are selected such as the following functions are exact solutions of the example under consideration Results are given in Tables 5 and 6.  Here, to improve the convergence of the method, the initial approximation was chosen sufficiently close to the exact one: x 0 1 (t) = 0.9 cos(t), x 0 2 (t) = 0.9 sin(t).

Remark 1.
The choice of a suitable initial guess is one of the main problems of all iterative methods. If we do not have any a priory information concerning the solution behavior, as an initial guess for solution of the system of nonlinear Equation (1), one can use an approximate solution of a low order of accuracy obtained by generalizing the direct method described in [25], Section 4.3. For more discussions concerning approximate local solution readers may refer to [14,17].

Storage System Analysis in Microgrid Using System of Volterra Equations
The proposed approach was tested on real annual climatological and load datasets from the resort village of Goryachinsk, located on the coast of Lake Baikal (Eastern Siberia, Russia). A more detailed description of the real data and 11 years of meteorological observations is shown in [24]. The historical datasets include mean hourly wind speed and time series of direct normal solar irradiance as well as typical electrical load profiles for a selected location. Based on this real dataset, we examined an isolated AC/DC network for one of the four resort villages. This isolated network uses three PVs (2 × 120 MW, 33 MW), two wind turbines (220 MW, 186 MW), two hydrogen storages (2 × 11 MW), two diesel generators (2 × 20 MW) and four batteries (2 × 150 MWh, 300 MWh, 180 MWh). The methodology developed in [22,23] let us consider the problem of storage system balancing problem within the hybrid network operational control approach developed in [24] for hybrid AC/DC power systems with renewable energy sources and four storages. As it was outlined in [22], the Volterra integral models can effectively describe the charge-discharge process of energy storage system. The direct problem is a well known ampere-hour integral model. The latter is considered as an inverse problem in terms of VIE (or system of VIEs in general setup) with respect to the instantaneous storage current i(τ) which is assumed to be positive for charge and negative for discharge. Here µ(τ) is the storage efficiency which can be function of state of charge (SOC) in turn.
For a community of microgrids (MGs) with storage systems, it is useful (see [24] for more details) to employ the following system of Volterra integral equations with jump discontinuous kernels with constrains. Here, we follow notations from [24].
where m j = {(t, τ)| α j−1 (t) < τ ≤ α j (t)}, α 0 (t) = 0, α n (t) = t, j = 1, 2, . . . , n. Here, m stands for number of grids; functions α j (t) show the proportions in which units in storage system are used in each subgrid. Here, n is number of units in storage system for i-th grid; the diagonal elements of the matrix K m×m shows efficiency of storage system of each grid, the remaining elements of the matrix show at the Local Level the coefficients of power flow from storage systems of other grids; f i RES (t) is the generation from renewable energy sources (RES); f i Load is predicted electric load of the community, f i AC/DC (t) is AC/DC power flow at the centralized level; v i max is maximum speed of the charge for i-th storage; E i min and E i max are constraints on the storage levels. The alternating power function (APF) based on x i (τ) is possible to find for each storage using this model.
Here, kernel models the nonlinear efficiency change caused by the night temperature drop (Figure 1). In (33), v i max = 3.5MW is the limitation on the maximum power with which the storage system in i-th grid can be charged and discharged, E i min = 0%, E i max = 100%. f i (t) is the disbalance between generation, losses, power flow and load to be compensated by storage system. The calculated APF and SOC are shown in Figure 2. As can be noted, the calculation results for the Volterra energy storage model, based on nonlinear equations, allow us to obtain a more accurate value of the state of charge for the energy storage and alternating power function in comparison with the linear Volterra model, which was verified based on the conventional linear model [22].

Conclusions
In this article, the efficient numerical methods for solution of the novel class of weakly regular linear and nonlinear systems of Volterra integral equations of the first kind are proposed. The iterative method for the solution of the system of nonlinear Volterra integral equations of the first kind with special jump discontinuous kernels employs the linearization of the integral vector-operator, using the Newton-Kantorovich iterative process. The convergence theorem for the proposed iterative process is formulated. The polynomial collocation method is also proposed to solve such systems. It is to be noted that the CESTAC (Controle et Estimation Stochastique des Arrondis de Calculs) method (see [27]) and the CADNA (Control of Accuracy and Debugging for Numerical Applications) library can be employed to control the rounding error estimation and detects possible numerical instabilities.
The efficiently constructed theory and numerical methods are illustrated in five examples and on the practical problem of microgrids balancing using a storage system. The model can effectively handle the nonlinear change of the storage system efficiency drop. In our future work, we will consider the different scenarios of the dynamic analysis of MGs with biomass-based generation and storage systems when efficiency depends on the desired alternating power functions and state of health degradation.