3D Nonparametric Neural Identification

This paper presents the state identification study of 3D partial differential equations (PDEs) using the differential neural networks (DNNs) approximation. There are so many physical situations in applied mathematics and engineering that can be described by PDEs; these models possess the disadvantage of having many sources of uncertainties around their mathematical representation. Moreover, to find the exact solutions of those uncertain PDEs is not a trivial task especially if the PDE is described in two or more dimensions. Given the continuous nature and the temporal evolution of these systems, differential neural networks are an attractive option as nonparametric identifiers capable of estimating a 3D distributed model. The adaptive laws for weights ensure the “practical stability” of the DNN trajectories to the parabolic three-dimensional (3D) PDE states. To verify the qualitative behavior of the suggested methodology, here a nonparametric modeling problem for a distributed parameter plant is analyzed.


Introduction
1.1.3D Partial Differential Equations.Partial differential equations (PDEs) are of vast importance in applied mathematics, physics, and engineering since so many real physical situations can be modelled by them.The dynamic description of natural phenomenons are usually described by a set of differential equations using mathematical modeling rules [1].Almost every system described in PDE has already appeared in the one-and two-dimensional situations.Appending a third dimension ascends the dimensional ladder to its ultimate rung, in physical space at least.For instance, linear second-order 3D partial differential equations appear in many problems modeling equilibrium configurations of solid bodies, the three-dimensional wave equation governing vibrations of solids, liquids, gases, and electromagnetic waves, and the three-dimensional heat equation modeling basic spatial diffusion processes.These equations define a state representing rectangular coordinates on R 3 .There are some basic underlying solution techniques to solve 3D PDEs: separation of variables and Green's functions or fundamental solutions.Unfortunately, the most powerful of the planar tools, conformal mapping, does not carry over to higher dimensions.In this way, many numerical techniques solving such PDE, for example, the finite difference method (FDM) and the finite element method (FEM), have been developed (see [2,3]).The principal disadvantage of these methods is that they require the complete mathematical knowledge of the system to define a mesh (domain discretization), where the functions are approximated locally.The construction of a mesh in two or more dimensions is a nontrivial task.Usually, in practice, only low-order approximations are employed resulting in a continuous approximation of the function across the mesh but not in its partial derivatives.The approximation discontinuities of the derivative can adversely affect the stability of the solution.However, all those methods are well defined if the PDE structure is perfectly known.Actually, the most of suitable numerical solutions could be achieved only if the PDE is linear.Nevertheless, there are not so many methods to solve or approximate the PDE solution when its structure (even in a linear case) is uncertain.This paper suggests a different numerical solution for uncertain systems (given by a 3D PDE) based on the Neural Network approach [4].

Application of Neural Networks to Model PDEs.
Recent results show that neural networks techniques seem to be very effective to identify a wide class or systems when we have no complete model information, or even when the plant is considered as a gray box.It is well known that radial basis function neural networks (RBFNNs) and MultiLayer Perceptrons (MLPs) are considered as a powerful tool to approximate nonlinear uncertain functions (see [5]): any continuous function defined on a compact set can be approximated to arbitrary accuracy by such class of neural networks [6].Since the solutions of interesting PDEs are uniformly continuous and the viable sets that arise in common problems are often compact, neural networks seem like ideal candidates for approximating viability problems (see [7,8]).Neural networks may provide exact approximation for PDE solutions; however, numerical constraints avoid this possible exactness because it is almost impossible to simulate NN structures with infinite number of nodes (see [7,9,10]).The Differential Neural Network (DNN) approach avoids many problems related to global extremum search converting the learning process to an adequate feedback design (see [11,12]).Lyapunov's stability theory has been used within the neural networks framework (see [4,11,13]).The contributions given in this paper regard the development of a nonparametric identifier for 3D uncertain systems described by partial differential equations.The method produces an artificial mathematical model in three dimensions that is able to describe the PDEs dynamics.The required numerical algorithm to solve the non-parametric identifier was also developed.

3D Finite Differences Approximation
The problem requires the proposal of a non-parametric identifier based on DNN in three dimensions.The problem here may be treated within the PDEs framework.Therefore, this section introduces the DNN approximation characteristics to reconstruct the trajectories profiles for a family of 3D PDEs.
Consider the set of uncertain second-order PDEs here a noise in the state dynamics.This PDE has a set of initial and boundary conditions given by In (1) where L is positive constant and g 2 = (g, g) is used just to ensure that there exists some δ > 0 such that the state equation γ = g(γ, t) with γ(t 0 ) = γ 0 has a unique solution over [t 0 , t 0 + δ] (see [14]).The norm defined above stands for the Sobolev space defined in [15] as follows.
Definition 1.Let Ω be an open set in R n , and let ν ∈ C m (Ω).Define the norm of ν(γ) as . This is the Sobolev norm in which the integration is performed in the Lebesgue sense.The completion of the space of function ν(γ) ∈ C m (Ω): ν m,p < ∞ with respect to • m,p is the Sobolev space H m,p (Ω).For p = 2, the Sobolev space is a Hilbert space (see [14,15]).
Below we will use the norm (4) for the functions u(•, t) for each fixed t.

Numerical Approximation for Uncertain
Functions.Now consider a function h 0 (•) in H m,2 (Ω).By [16], h 0 (•) can be rewritten as where {Ψ i jk (γ)} are functions constituting a basis in H m,2 (Ω).Last expression is referred to as a vector function series expansion of h 0 (γ, θ * ).Based on this series expansion, an NN takes the following mathematical structure: that can be used to approximate a nonlinear function where Following the Stone-Weierstrass Theorem [17], if is the NN approximation error, then for any arbitrary positive constant ε there are some constants The main idea behind the application of DNN [11] to approximate the 3D PDEs solution is to use a class of finite-difference methods but for uncertain nonlinear functions.So, it is necessary to construct an interior set (commonly called grid or mesh) that divides the subdomain equidistant sections, each one of them (Figure 1) defined as (x i , y j , z k ) in such a way that x 0 = y 0 = z 0 = 0 and Using the mesh description, one can use the next definitions: Analogously, we may consider the other cases (u xx , u y , u yy , u z , u zz , u xy , u xz , u yz ).Using the mesh description and applying the finite-difference representation, one gets and it follows for all cases such that the (Δx, Δy, Δz)-approximation of the nonlinear PDE (1) can be represented as ui,j,k (t 2.2.3D Approximation for Uncertain PDE.By simple adding and subtracting the corresponding terms, one can describe (1) as where u t = u t (x, y, z, t), u = u(x, y, z, t), σ = σ(x, y, z), u x = u x (x, y, z, t), u xx = u xx (x, y, z, t), the same for σ 2 , ϕ i , γ i , ψ i , u y , u yy , u z , u zz , u xy , u yz , u xz , and ) any constant matrices and the set of sigmoidal functions have the corresponding size (σ(x, y, z ) and are known as the neural network set of activation functions.These functions obey the following sector conditions: which are bounded in x, y, and z, that is, Following the methodology of DNN [11] and applying the same representation to (12), we get for each i ∈ (1, . . ., L), j ∈ (1, . . ., M), k ∈ (1, . . ., N) the following robust adaptive non-parametric identifier: In the sum we have that s = 1, 3, φ represents functions σ, ϕ, γ, ψ, and U can be taken as the corresponding u i, j,k , In this equation the term f i, j,k (t), which is usually recognized as the modeling error, satisfies the following identify, and here, it has been omitted the dependence to x i , y j , z k of each sigmoidal function: where , 11, φ s , φ represents functions σ, ϕ, γ, ψ and s = 1, 3, U(t) represents the corresponding (u i, j,k , We will assume that the modeling error terms satisfy the following.
Assumption 2. The modeling error is absolutely bounded in Ω: Assumption 3. The error modeling gradient is bounded as , where s = x, y, z and f i, j,k r (r = 1, 3) are constants.

DNN Identification for Distributed Parameters Systems
3.1.DNN Identifier Structure.Based on the DNN methodology [11], consider the DNN identifier for all i = 0, . . ., L; u −1 (t) = u −2 (t) = 0, where φ represents activation functions σ, ϕ, γ, and ψ, s = 1, 3, U is each one of the states u i, j,k , , and A i, j,k ∈ R n×n is a constant matrix to be selected, u i, j,k (t) is the estimate of u i, j,k (t).Obviously that proposed methodology implies the designing of individual DNN identifier for each point x i , y j , z k .The collection of such identifiers will constitute a DNN net containing N × M connected DNN identifiers working in parallel.Here σ 1 (x i , y j , z k ), ϕ 1 (x i , y j , z k ), ϕ 2 (x i , y j , z k ), ϕ 3 (x i , y j , z k ), γ 1 (x i , y j , z k ), γ 2 (x i , y j , z k ), γ 3 (x i , y j , z k ), ψ 1 (x i , y j , z k ), ψ 2 (x i , y j , z k ), ψ 3 (x i , y j , z k ), and σ 2 (x i , y j , z k ) are the NN activation vectors.This means that the applied DNN-approximation significantly simplifies the specification of •) which now are constant for any x i , y j , z k fixed.

Learning Laws for Identifier's Weights.
For each i = 0, . . ., L, j = 0, . . ., M, k = 0, . . ., N, define the vector-functions defining the error between the trajectories produced by the model and the DNN-identifier as well as their derivatives with respect to x, y, and z for each i, j, k: Let W i, j,k r (t) ∈ R n , r = 1, 11 be time-variant matrices.These matrices satisfy the following nonlinear matrix differential equations: where Ω h (x i , y j , z j ) = S l (x i , y j , z j ) (h = 1, 10, l = 1, 3) represents the corresponding sigmoidal functions σ l (x i , y j , z j ), ϕ l (x i , y j , z j ), γ l (x i , y j , z j ), and ψ l (x i , y j , z j ).Here with positive matrices K r (r = 1, 11) and P i, j,k , S i, j,k 1 , and S i, j,k 2 (i = 0, N; j = 0, M) which are positive definite solutions (P i, j > 0, S i, j 1 > 0, S i, j 2 > 0) and S i, j,k 3 (i = 0, N; j = 0, M, k = 0, L) of the algebraic Riccati equations defined as follows: Ric P i, j,k := P i, j,k A i, j,k + A i, j,k P i, j,k where B can be P, S 1 , S 2 , and S 3 and b = (7, 14, 21, 28).Matrices Q i, j,k B have the form where B can be P, S 1 , S 2 , or S 3 , representing the partial derivative; for S 1 it is with respect to x, for S 2 with respect to y, and for S 3 with respect to z, and Λ −1 l (l = 1, 46), where

The special class of Riccati equation
has a unique positive solution P if and only if the four conditions given in [11] (page 65, chapter 2 Nonlinear System Identification: Differential Learning) are fulfilled (see [11]): (1) matrix A is stable, (2) pair (A, R 1/2 ) is controllable, (3) pair (Q 1/2 , A) is observable, and ( 4) matrices (A, Q, R) should be selected in such a way to satisfy the following inequality: which restricts the largest eigenvalue of R guarantying the existence of a unique positive solution.The main result obtained in this part is in the practical stability framework.

Practical Stability and Stabilization
The following definition and proposition are needed for the main results of the paper.Consider the following ODE nonlinear system: with z t ∈ R n , and v t ∈ R m , and t an external perturbation or uncertainty such that t 2 ≤ + .
Definition 4 (Practical Stability).Assume that a time interval T and a fixed function v * t ∈ R m over T are given.Given ε > 0, the nonlinear system (30) is said to be ε-practically stable over T under the presence of t if there exists a δ > 0 (δ depends on ε and the interval T) such that z t ∈ B[0, ε], for all 0 ≤ t, whenever z t0 ∈ B[0, δ].
Similarly to the Lyapunov stability theory for nonlinear systems, it was applied the aforementioned direct method for the ε-practical stability of nonlinear systems using-practical Lyapunov-like functions under the presence of external perturbations and model uncertainties.Note that these functions have properties differing significantly from the usual Lyapunov functions in classic stability theory.
The subsequent proposition requires the following lemma.

Lemma 5. Let a nonnegative function V (t) satisfying the following differential inequality:
where α > 0 and β ≥ 0. Then with μ = β/α and the function [•] + defined as Proof.The proof of this lemma can be obtained directly by the application of the Gronwall-Bellman Lemma.

Proposition 6. Given a time interval T and a function v(•)
over a continuously differentiable real-valued function V (z, t) satisfying V (0, t) = 0, for all t T, is said to be ε-practical Lyapunov-like function over T under v if there exists a constant α > 0 such that with H a bounded nonnegative nonlinear function with upper bound H + .Moreover, the trajectories of z t belong to the zone Proof.The proof follows directly from Lemma 5.
Definition 7. Given a time interval T and a function v(•) over T, nonlinear system (30) is ε-practically stable, T under v if there exists an ε-practical Lyapunov-like function V (x, t) over T under v.

Identification Problem Formulation
The state identification problem for nonlinear systems (13) analyzed in this study, could be now stated as follows.
Problem.For the nonlinear system, given by the vector PDE (20), to study the quality of the DNN identifier supplied with the adjustment (learning) laws ( 22), estimate the upper bound of the identification error δ given by (with P i, j,k from (24)) and, if it is possible, to reduce to its lowest possible value, selecting free parameters participating into the DNN identifier (A, K r , r = 1, 11).
This implicates that the reduction of the identification error δ means that the differential neural network has converged to the solution of the 3D PDE; this can be observed in the matching of the DNN to the PDE state.

Main Result
The main result of this paper is presented in the following theorem.
Proof.The proof is given in the appendix.

Simulation Results
7.1.Numerical Example.Below, the numerical smulation shows the qualitative illustration for a benchmark system.Therefore, consider the following three-dimensional PDE described as follows: where c 1 = c 2 = c 3 = 0.15.It is assumed that there is access to discrete measurements of the state u(x, y, z, t) along the whole domain, which is feasible in practice.ξ(x, y, z, t) is a noise in the state dynamics.This model will be used just to generate the data to test the 3D identifier based on DNN.Boundary conditions and initial conditions were selected as follows: u(0, 0, 0, t) = rand(1), u x (0, 0, 0, t) = 0, u y (0, 0, 0, t) = 0, u z (0, 0, 0, t) = 0.
The trajectories of the model can be seen in Figure 3 as well as the estimated state, produced by the DNN identifier.The

Tumour Growth
Example.The mathematical model of the brain tumour growth is presented in this section based on the results of [18].Here the diffusion coefficient is considered as a constant.Let consider the following three-dimensional  parabolic equation of the tumour growth described as follows: u t x, y, z, t = −Pu x x, y, z, t − Ru y x, y, z, t − Su z x, y, z, t + Qu xx x, y, z, t + Qu yy x, y, z, t + Qu zz x, y, z, t where u(x, y, z, t) is the growth rate of a brain tumour, Q is the diffusion coefficient, W = (P, R, S) is the drift velocity field, Γ = Γ(u) is the proliferation coefficient, and L = L(u) is the decay coefficient of cells.It is assumed that there is access to discrete measurements of the state u(x, y, z, t) along the whole domain, which is feasible in practice by PET-CT (Positron emission tomography-computed tomography) technology.The domain 0 ≤ x ≤ 1, 0 ≤ y ≤ 1, and 0 ≤ z ≤ 1 This model will be used just to generate the data to test the 3D identifier based on DNN.Boundary conditions and initial conditions were selected as follows: u(0, 0, 0, t) = 200 ± 20μ, u x (0, 0, 0, t) = 0, u y (0, 0, 0, t) = 0, u z (0, 0, 0, t) = 0. (41) The trajectories of the model and the estimate state produced by the DNN identifier can be seen in Figure 7.The dissimilarity between both trajectories depends on the learning period required for adjusting the DNN identifier.The error between trajectories produced by the model and the proposed identifier is close to zero almost for all x, y, z and all t that shows the efficiency of the identification process provided by the suggested DNN algorithm is shown in Figure 8. In

Conclusions
The adaptive DNN method proposed here solves the problem of non-parametric identification of nonlinear systems (with uncertainties) given by 3D uncertain PDE.Practical stability for the identification process is demonstrated based on the Lyapunov-like analysis.The upper bound of the identification error is explicitly constructed.Numerical examples demonstrate the estimation efficiency of the suggested methodology.

Theorem 8 .Figure 2 :Figure 3 : 10 −Figure 4 :
Figure 2: Numerical trajectory produced by the mathematical model described by 3D partial differential equation at time 10 s along the whole domain.

Figure 5 : 2 Figure 6 :
Figure 5: Comparison between 3D PDE trajectories and 3D DNN approximation for coordinate z at time 10 s.

Figure 7 :
Figure 7: Numerical trajectory produced by the mathematical model described by 3D partial differential equation at time 46 days along the whole domain.
Figures 9 and 10 there are shown the trajectories of the PDE and the DNN for the coordinate z and the index error of the PDE and the DNN at 46 days for coordinates x, z, respectively.
Let us consider a vector function g(t) ∈ H to be a piecewise continuous in t.By L ∞ (a, b; H) we denote the set of H-valued functions g such that (g(•), u) is Lebesgue measurable for all u ∈ H and ess max t∈[a,b] g(γ, t) < ∞, γ ∈ R n .Suppose that the nonlinear function g(γ, t) satisfies the Lipschitz condition g(γ, t) −