Physics-Based Neural Network Methods for Solving Parameterized Singular Perturbation Problem

Abstract

This work is devoted to the description and comparative study of some methods of mathematical modeling. We consider methods that can be applied for building cyber-physical systems and digital twins. These application areas add to the usual accuracy requirements for a model the need to be adaptable to new data and the small computational complexity allows it to be used in embedded systems. First, we regard the finite element method as one of the “pure” physics-based modeling methods and the general neural network approach as a variant of machine learning modeling with physics-based regularization (or physics-informed neural networks) and their combination. A physics-based network architecture model class has been developed by us on the basis of a modification of classical numerical methods for solving ordinary differential equations. The model problem has a parameter at some values for which the phenomenon of stiffness is observed. We consider a fixed parameter value problem statement and a case when a parameter is one of the input variables. Thus, we obtain a solution for a set of parameter values. The resulting model allows predicting the behavior of an object when its parameters change and identifying its parameters based on observational data.

1. Introduction

The development of technologies that make it possible to obtain data from environments, observed objects, and processes and process them with sufficient speed and efficiency leads to the widespread implementation of cyber-physical systems (CPSs) [1]. As a part of the operating of these systems, the received information (data) is synchronized with a corresponding physical or information object. The result of such synchronization is a digital twin (DT) of an object, i.e., a model of a complex system. The DT continues to synchronize with the modeled object and reflect its dynamic characteristics that change over time. Such models can be used to conduct digital experiments, implement alternative external scenarios, and predict the behavior of an object in the future without affecting or changing a real prototype.

The issue of building DTs and the simultaneous development of artificial intelligence technologies have led to the emergence of so-called “pure” data-driven modeling “black boxes”, in which a model is created on the basis of deep learning using big data. While the possibility of obtaining such data has been growing in recent times, the cost of storing and processing data has increased along with a rise in data volumes. In addition, for some complex industrial and technological processes, the acquisition of a large amount of data itself can be accompanied by high expenditure and, sometimes, a lack of measurement possibilities. In these cases, the data can be utilized to verify a model built by classical modeling methods based on the physics of an object, including numerical methods for solving corresponding problems of mathematical physics.

The intention to preserve the advantages and eliminate the disadvantages of both types of methods is reflected in the new classifications of existing modeling paradigms [2,3,4]. Physical laws provide interpretability and extrapolation properties to data-driven models which, in turn, connect theoretical models to real objects. We do not aim at repeating or improving the overviews of different models. However, the description of various modeling classes seems essential. Ref. [2] reviews the concepts of DT modeling and highlights the value of the DT in the modern world and its application in various sectors of the economy. Among the models combining physics and data-based modeling, [2] points out four categories formed by the mutual integration of the following modeling approaches: physics-based modeling, data-driven modeling, and big data methods. A similar classification can be found in [5]. The article [3] focuses on the CPS application and provides a detailed classification of hybrid models, distinguishing among physics-based machine learning modeling classes of algorithms that differ in the nature of interaction between physics-based and machine learning (ML) approaches to modeling. Firstly, these are the methods that utilize physics theory for the preprocessing of data, which then arrive at an input of a neural network in ML. The second type of model comprises those in which prior knowledge of the physics of an object is reflected in an architecture of a network. Thirdly, one of the most popular hybrid modeling methods is underlined, which consists of using physics equations as an additional regularization term in the loss function of an ML model. Similar classes are described in [6]. In [3], models in which activation functions are physically interpretable also act as a separate class. In [4,7], in addition to those mentioned above, embedding approaches are pointed out, in which data-driven models are used as components of classical physics-based models to speed up, simplify, or replace a numerical solution.

The authors of one of the most cited recent papers on hybrid modeling of complex systems [8] mention that, among others, they have been inspired by the work [9] dedicated to the use of neural networks to solve ordinary and partial differential equations. Let us note that in our work [10], published in 2005, we talk about solving problems of mathematical physics by using both neural networks and real measurements. From the point of view of modern classification, we can regard this approach as a hybrid model with physics-based regularization.

In this article, we consider our previously developed methods of hybrid modeling and their modifications in the context of the modern classification of paradigms for modeling complex dynamic systems. We present the principles of building physics-based network models naturally capable of further adaptation and refinement using measurement data and other additional data.

We have chosen the problem of modeling a system described by a boundary value problem with an unknown parameter as a specific issue of this paper. Problems involving differential equations with unknown parameters frequently arise in the modeling and control of CPS. These parameters can be regarded as dynamic and changing with time according to an unknown law. When describing complex systems, the parameter can characterize the conditions in which the object exists, or a specific part of this object. Thus, the parameter can influence a type of the differential problem under question. It can be considered a hard problem, an ill-posed problem, a problem with a small parameter at the highest derivative, and, accordingly, a problem with a boundary layer, and other types. Let us note that, when solving a problem with an unknown parameter, we do not set to ourselves the task of identifying this parameter, dubbed data-driven discovery of differential equations. As a result of applying methods described below, we obtain a parametric solution that is universal for various parameter values belonging to a predetermined interval that is used in the construction and training of a neural network solution. The resulting pre-trained model can be further adapted by training on the measured data. Based on these measurements, unknown parameters can also be estimated by any means or adjusted using additional training of the network. A similar problem is solved in [11], where simulating a lake’s temperature is considered. The authors of [11] use the loss function with physics-based regularization to pre-train a recurrent neural network model before adapting it by means of measurements of a certain lake. They also suggest using the finite element method data as a training dataset.

It is necessary to make a remark about the accuracy of a solution to the problem. Nowadays, there are two fundamentally different approaches to mathematical modeling. The first one is traditional and involves passing from a simulated object to a model in the form of a differential equation (system of equations) with additional conditions (initial, boundary, etc.). Further, this model is considered the object itself and an approximate solution is constructed based on it. The higher the accuracy of such a solution, the more accurate the model of the original real object appears to be. The second approach, which we adhere to, is that the differential model built on the basis of physical principles is certainly approximate. For real complex technical objects, the accuracy of similar models is often low. Therefore, the solution of differential equation(s) with high accuracy has no practical sense. More promising is the construction of an adaptive model that corresponds to the differential equation and additional conditions with reasonable accuracy and the subsequent adaptation of this model to the data of object observations. This approach is especially relevant in a situation when the properties of an object can change during its operation.

The paper [12] considers hybrid systems, which are a class of problems that cyclically switch between two phases: the first is described by a stiff system of differential algebraic equations of motion (interaction with the ground), the second corresponds to the flight phase and is specified by a non-stiff system. Using a specific problem (hopping single-legged robots) as an example, the authors of [12] compare various ODE (ordinary differential equation) solvers, including those common in the Matlab package. It is noted that although different solvers work well for different types of problems, due to the constant restart of the numerical solver, the accuracy of the solution diminishes and the cost of computation increases. This problem could be solved with a simple universal approach. Ref. [13] considers the need to solve a parameterized problem that also arises when studying a wide range of complex physical phenomena, when, while modeling many different possible implementations of the system, large requirements for computational resources lead to the forced use of the apparatus for reducing parametric models. The authors underline the importance of such parametric models in global problems of design, control, optimization, and quantification of uncertainty. The development of a controller that works effectively for all values of model parameters is regarded as an urgent task. Ref. [14] proposes a solution to a parameterized differential problem by the FEM for an interval parameter based on interval arithmetic. However, as a result they obtain only the corresponding interval boundaries of the solution. Earlier, the authors of the given article studied the construction of a unified solution to one parameterized differential problem, which, for some values of the parameter, is stiff [15]. The unified neural network approach was applied using various additional information of heterogeneous types to improve the quality of an approximate solution. In [16], we successfully applied this approach to solve various types of problems with parameters. A stiff differential problem and a differential algebraic problem with different numbers of solutions were solved. In the catalyst problem, we solved the problem in a region that is wider than the domain of existence of the solution to the differential equation. The authors of [17] also used trainable neural networks with a single hidden layer to obtain a differentiable closed-form solution for a class of parameterized initial-value systems with varying degrees of stiffness. Ref. [18] highlights the growing interest in the development of numerical methods for solving similar problems and suggests a method based on a rational spectral collocation. Adaptive grid methods are considered in [19,20,21,22]. Ref. [23] proposes a special neural network for solving linear singularly perturbed differential equations with initial and boundary values. We note the advantage of low computational costs of the application of neural network models after the training phase. However, when describing systems with uncertainty, the structure of the differential problem itself may often be unknown in advance. Thus, the need arises for building an analytical solution that does not require real-time training and reflects the model under consideration for different parameter values. This approach related to physics-based network architectures is presented in this paper as an extension of the multilayered techniques [24].

We have planned the structure of this manuscript as follows. Section 2 discusses the methods of constructing the models of systems described in the form of a classical non-parameterized boundary value problem. Through expansion in basis functions, the basic principles of the finite element method (FEM) [25] as one of the “pure” physics-based modeling methods, and the general neural network approach [10] as a variant of ML modeling with physics-based regularization, are presented. In Section 2.3, the ML model is proposed to be improved by using prior knowledge of the physics of an object in the form of a training data set. The values of an FEM solution in the nodes are utilized instead of measurement results for additional training of the neural network model. A less popular modeling class, named physics-based network architectures, is suggested in Section 2.4 as a special multilayered model based on classical numerical methods. Section 3 is devoted to methods of constructing models in cases where the parameters included in the corresponding equations may not be known. In particular, guidance for using the methods described in Section 2.3 and Section 2.4 for parameterized problem setting is discussed. We want to emphasize that this certain problem statement is our major interest in this work. Approaches to solving the problem with a fixed parameter are presented to facilitate an understanding of methods for solving parameterized problems. Section 4 features the parameterized singular perturbation problem and presents the results of applying the previously described methods to the model construction. Section 5 contains conclusions and considerations for future research.

2. Materials and Methods. Ordinary Problems

To clarify the methods, let us formulate as an example a boundary value problem

$$\left\{\begin{array}{c}Au=g,u=u\left(\mathbf{x}\right),\mathbf{x}\in \Omega \subset {\mathbb{R}}^p;\hfill \\ {Bu|}_{\Gamma }=f;\hfill \end{array}\right.$$

(1)

where A is a differential operator, i.e., an algebraic expression containing derivatives of an unknown function u, B is some operator specifying the boundary conditions, and $\Gamma$ is the boundary of some domain.

We seek an approximate solution in the form of an expansion in basis functions $v_i$

$$u_N\left(\mathbf{x}\right)=\sum _{i=1}^N{c}_i{v}_i\left(\mathbf{x}\right).$$

(2)

Searching for a solution in form (2) is usually called the Galerkin method. Its special cases include different variants of the FEM and neural network models, in which equations based on physics are used as an additional term for the regularization of the loss function. In this case, this models differ in the features of the selection of basis functions and parameters $c_i$.

2.1. Methods with Fixed Basis Functions

First, let us distinguish a group of methods in which basis functions $v_i\left(\mathbf{x}\right)$ are fixed and only parameters $c_i$ are adjusted. The extreme learning machine (ELM) algorithm [26] for single hidden layer feedforward neural networks (SFLNs) randomly chooses internal weights and biases of basis functions and then analytically obtains the output weights $c_i$. Ref. [27] proposes the incremental method for SFLNs and also shows its efficiency for SLFNs with piecewise continuous basis functions.

A similar approach is typical for the simplest version of the FEM. There are two conventional approaches to finding the parameters $c_i$. The first, the most common of them, is that the space with the dot product $a\left(\mathbf{x}\right)·b\left(\mathbf{x}\right)$ is regarded as the space in which the solution is sought. The simplest version of such a dot product is the product in the integral form $a\left(\mathbf{x}\right)·b\left(\mathbf{x}\right)={\int }_{\Omega }a\left(\mathbf{x}\right)b\left(\mathbf{x}\right)d\mathbf{x}$. In this way, the parameters $c_i$ are found from the condition of orthogonality of the residual of Equation (1) to all the basis functions

$$(A{u}_N\left(\mathbf{x}\right)-g\left(\mathbf{x}\right))·v_i\left(\mathbf{x}\right)=0.$$

(3)

It is easy to ascertain that they can be found as a solution to the system

$$\sum _{i=1}^N{c}_{i}A{v}_i\left(\mathbf{x}\right)·v_j\left(\mathbf{x}\right)=g\left(\mathbf{x}\right)·v_j\left(\mathbf{x}\right).$$

(4)

If the integral cannot be calculated explicitly, it can be calculated by using cubature formulas or through the mean value

$$a\left(\mathbf{x}\right)·b\left(\mathbf{x}\right)\cong \frac{1}{M}\sum _{j=1}^{M}a\left({\xi }_j\right)b\left({\xi }_j\right).$$

(5)

Basis functions usually satisfy the homogeneous boundary condition $B{v}_i{|}_{\Gamma }=0$. To satisfy the inhomogeneous boundary condition, additional term $v_0\left(\mathbf{x}\right)$ is added to sum (2). This term satisfies the inhomogeneous boundary condition $B{v}_0{|}_{\Gamma }=f$. The sum takes the form $u_N\left(\mathbf{x}\right)=v_0\left(\mathbf{x}\right)+{\sum }_{i=1}^N{c}_i{v}_i\left(\mathbf{x}\right)$.

Accordingly, the substitution $u_N\left(\mathbf{x}\right)$ into (3) yields a linear system of the form

$$\sum _{i=1}^N{c}_{i}A{v}_i\left(\mathbf{x}\right)·v_j\left(\mathbf{x}\right)=(g\left(\mathbf{x}\right)-A{v}_0\left(\mathbf{x}\right))·v_j\left(\mathbf{x}\right),j=1,\dots ,M.$$

(6)

Or, for the discrete form of the loss function,

$$\sum _{i=1}^N{c}_i\sum _{k=1}^{M}A{v}_i\left({\xi }_k\right)·v_j\left({\xi }_k\right)=\sum _{k=1}^M(g\left({\xi }_k\right)-A{v}_0\left({\xi }_k\right))·v_j\left({\xi }_k\right),j=1,\dots ,M.$$

(7)

This approach is referred to as the collocation method.

The above basis functions are usually chosen so that their supports intersect only for a small number of functions, while the matrix of system (4) is sparse, which greatly speeds up the solution procedure.

The second approach is that the quality of the solution to problem (1) is characterized by a loss function (dubbed the error functional or error function) J, the minimization of which leads to the needed parameters $c_i$.

As shown above, we can use the integral form of the loss function, for example, $(A{u}_N\left(\mathbf{x}\right)-g\left(\mathbf{x}\right))·(A{u}_N\left(\mathbf{x}\right)-g\left(\mathbf{x}\right))$.

If it is difficult to calculate the dot product explicitly, then we can use the discrete form

$$J=\sum _{j=1}^M{(A{u}_N\left({\xi }_j\right)-g\left({\xi }_j\right))}^2.$$

(8)

The finite element method is characterized by using basis functions (elements) with a compact support localized in small subdomains of $\Omega$. In this case, the supports of the basis functions intersect and cover the entire domain $\Omega$. With this method, if the supports of the basis functions intersect only with a small number of supports of other basis functions, the matrix of the system turns out to be sparse.

For piece-wise linear elements, direct application of Formula (4) may be impossible due to the absence of derivatives of the required order. In several cases, this problem can be solved by applying the Stokes formula or by the corresponding version of the integration by parts formula.

2.2. Physics-Based Neural Network Approach

When solving stiff problems, the previously considered approach with fixed basis functions does not work well, and these functions have to be adjusted during the solving process. In [28], the ELM learning algorithm with a tunable basis function is studied.

At first glance, the general neural network approach is similar to the previous ones, only expansion (2) is specified by the expression

$$u_N\left(\mathbf{x}\right)=\sum _{i=1}^N{c}_{i}v(\mathbf{x},{\mathbf{a}}_i).$$

(9)

By solving the optimization problem for a certain loss function, both linear parameters $c_i$ and vector parameters (weights and biases) ${\mathbf{a}}_i$ are adjusted. Functions typical of neural networks are used as basis functions. These are usually radial basis functions (RBFs) or perceptron-type functions with sigmoid basis functions.

2.2.1. Multilayer Perceptron

This type of neural network is the most widely used and researched. Let us describe its structure. Let the input be an m-dimensional vector. Linear combinations of coordinates of the input of the form $(w,x)$ are fed to the input of the first layer of neurons. The coefficients of these combinations are named the weights of the first layer. Each neuron acts as a one-dimensional non-linear function dubbed an activation function. Linear combinations of outputs of neurons are fed to the next layer, and linear combinations of outputs of neurons of the last layer form the output of the network. An additional neuron, whose output always equals 1, is often added to all or some of the layers. Its meaning is the same as adding bias to linear regression, i.e., in fact, the subtraction of a constant value. Experience has shown that this greatly improves network training. In addition, an input is often added, whose value always equals 1. Usually, an activation (basis) function is a function that approximates in form to $\mathrm{sign}\left(x\right)$. Its meaning is to split the space into two half-spaces by hyperplanes $(\mathbf{w},\mathbf{x})=0$ and $(\mathbf{w},\mathbf{x})+w_0=0$ as long as a bias is used. In one of the half-spaces, the output of the neuron is equal to +1; in the other it is equal to −1. If we take not only one neuron, but a layer of neurons, the entire space is divided by hyperplanes into subsets, each of which corresponds to its own set of outputs of neurons of this layer. In this way, the network output is obtained as a piecewise constant function. Usually, instead of function $\mathrm{sign}\left(x\right)$, smooth functions are used, which makes it possible to calculate the derivatives and apply gradient methods to train the network. In this case, it is advisable that the calculation of the derivatives of the activation function requires a minimum of additional operations. For such functions, the network does not produce a piecewise constant, but a smoothed dependence. Here are some examples.

The activation function recommended in the works of A.N. Gorban and his students has the form $\phi \left(x\right)=\frac{x}{1+\left|x\right|}$ and ${\phi }^{\prime }\left(x\right)=\frac{1}{{(1+|x\left|\right)}^2}={(1-|\phi \left(x\right)\left|\right)}^2$. One function is more commonly used, namely $\phi \left(x\right)=\mathrm{th}\left(x\right)$. Here, ${\phi }^{\prime }\left(x\right)=\frac{1}{{\mathrm{ch}}^2\left(x\right)}=1-{\phi }^2\left(x\right)$. An asymmetric version of this function is sometimes used $\phi \left(x\right)=\frac{expx}{1+expx}$, where ${\phi }^{\prime }\left(x\right)=\frac{expx}{{(1+expx)}^2}=\phi \left(x\right)(1-\phi \left(x\right))$. Other functions are also applied, for example, $\phi \left(x\right)=arctan\left(x\right)$, ${\phi }^{\prime }\left(x\right)=\frac{1}{1+x^2}$.

Let us describe the functioning of the multilayer perceptron more formally. Denote the vector of inputs of the l-th layer ${\mathbf{y}}_{l-1}$, and the vector of outputs ${\mathbf{x}}_l$. Then the network is described by recurrence relations ${\mathbf{y}}_l={\mathbf{W}}_l{\mathbf{x}}_l$, ${\mathbf{x}}_l={\phi }_l\left({\mathbf{y}}_{l-1}\right)$. Here ${\mathbf{x}}_0=\mathbf{x}$ is a network input and ${\mathbf{y}}_L=\mathbf{f}\left(\mathbf{x}\right)$ is a network output. ${\mathbf{W}}_l$ is a matrix of weights of the l-th layer, and ${\phi }_l$ is the activation function, which acts coordinate-wise.

Application of gradient methods of network training, i.e., minimizing the error, requires calculating the derivatives of the network output with respect to the weights, which are easy to obtain using the formula for differentiating a complex function.

To determine the desired derivatives, calculations are performed using the formula

$$\frac{\partial \mathbf{f}}{\partial w_{ij}^{\left(l\right)}}={\mathbf{W}}_L{\mathbf{Z}}_L{\mathbf{W}}_{L-1}{\mathbf{Z}}_{L-1}\cdots {\mathbf{W}}_{l+1}\frac{\partial {\mathbf{W}}_l}{\partial w_{ij}^{\left(l\right)}}\mathbf{x}.$$

(10)

Matrices ${\mathbf{Z}}_L$ are calculated by counting the network outputs. In accordance with (10), computations begin at the last layer and move to the first. In the literature on neural networks, the calculation of derivatives by this algorithm is often combined with calculations by the gradient descent method, and this hybrid is called the back-propagation method. This algorithm is generalized in the form of automatic differentiation and is widely used in training neural networks, including in [8].

2.2.2. Radial Basis Function Networks

The main difference between this type of network and the perceptron is that each neuron is responsible for local approximation in the vicinity of some point, and not for the difference between the values of the function in half-spaces, as in the case of the perceptron. Hence the conclusion that the perceptron is best used to simulate problems with jumps, such as phase transitions and shock waves. RBF networks are preferred for applications where a sufficiently smooth solution is to be expected.

In this case, the activation function has the form $v_i=\phi ({\alpha }_i\left|\right|\mathbf{x}-{\mathbf{x}}_i\left|\right|)$. Here, ${\mathbf{x}}_i$ is the center weight, ${\alpha }_i$ is a scale parameter for each neuron of a hidden layer, and $\left|\right|\ast \left|\right|$ is some metric on $\mathbb{R}$. Parameters $c_i$, ${\alpha }_i$, and ${\mathbf{x}}_i$ are selected during network. Two fundamentally different approaches can be applied to the adjustment of centers.

We consider fixed centers in the distribution domain of the input. They can be chosen regularly in this area according to certain principles or randomly due to some probabilistic law. Each ${\mathbf{x}}_i$ acts the same as parameters $c_i$ and ${\alpha }_i$ included in the process of loss function optimization.

There are also two options for choosing parameters ${\alpha }_i$. First, ${\alpha }_i=\alpha$, where $\alpha$ is a certain size of the considered domain. This method is utilized, as a rule, if the centers are selected in a regular way. Alternatively, ${\alpha }_i$ is defined by some method of non-linear optimization in conjunction with other parameters. We can apply an intermediate option assuming that ${\alpha }_i=\alpha$ and adjusting $\alpha$ as one of the optimized parameters. There are a few options left for $c_i$. If they are found separately, then a linear system is obtained. If we mutually select all the coefficients, then we can use some non-linear optimization algorithm.

To solve problem (1) it is natural to choose the loss function in the form $J=J_1+\delta J_2$, where $J_1$ corresponds to obeying the equation and $J_2$ to obeying the boundary condition. Usually, due to the impossibility of calculating integrals explicitly, a discrete form is used. Thus, $J_1={\sum }_{j=1}^M{(A{u}_N\left({\xi }_j\right)-g\left({\xi }_j\right))}^2$ and $J_2={\sum }_{k=1}^K{(B{u}_N\left({\xi }_k^{\prime }\right)-f\left({\xi }_k^{\prime }\right))}^2.$

Our experience of applying neural networks [29] allows making some remarks about the choice of points ${\left\{{\xi }_j\right\}}_{j=1}^M$ and ${\left\{{\xi }_k^{\prime }\right\}}_{k=1}^K$. They can be chosen inside $\Omega$ (or in a wider set $\tilde{\Omega }\supset \Omega$ if it is necessary to ensure sufficient smoothness of an approximate solution up to the boundary) and on the boundary in a regular way; for example, uniformly in the case of a bounded domain or according to the normal probability law if it is unbounded.

It is often more appropriate to utilize a set of points selected according to the probabilistic distribution law and re-generated after a certain number of learning epochs (optimization steps) using a constant (or other) probability density. This ensures more stable learning.

We had previously shown (see, for example, [15]) that the optimal number of steps in the learning process between point regeneration is from 3 to 5. In addition, this approach allows us to control the quality of training using standard statistical procedures. In this way, the number of sample points can be significantly reduced without compromising accuracy [29]. At the same time, the learning process becomes more stable and does not get stuck in local minima. In fact, not one function but a sequence of functions is optimized, and each is a discrete approximation to the integral form of the loss function.

If the model must satisfy additional conditions, for example, measurement data $u\left({\xi }_l^{″}\right)=u_l$, a term $\tilde{\delta }{J}_3$ is added to the loss function. Thus, $J=J_1+\delta J_2+\tilde{\delta }{J}_3$, where $J_3={\sum }_{l=1}^L\left(u\left({\xi }_l^{″}\right)-u_l\right)$.

We have dedicated a number of works to the study of modeling problems using physics-based neural networks, see [10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30]. A detailed description of the methods can be found in [29]. In recent years, physics-based ML modeling methods have included different types of neural networks. Ref. [11] considers recurrent neural networks and applies physics-based regularization. Ref. [31] solves the problem of modeling non-linear dynamical systems described by partial differential equations without initially given training sets. They utilize data generated on the basis of equations to train auto-regressive networks.

Note that the FEM can be regarded as a special case of the RBF network method (ELM). Most often, when using the FEM, piece-wise linear approximation is applied. In this case, triangulation is performed; that is, the domain is divided into triangles (the boundary is approximated by a broken line) and a linear function is constructed on each triangle. On the boundary of the triangles, the functions are continually joined. To obtain a representation of the solution in the form of an expansion in basis functions (2), one can consider the behavior of the function along the rays emanating from a certain center $c_i$. In the two-dimensional case, we have $v_i=v_i(\rho ,\psi )$, where $(\rho ,\psi )$ are the polar coordinates of the vector $\mathbf{x}-{\mathbf{c}}_i$, and $\mathbf{x}(\rho ,\psi )$ is the current point. If the position of the boundary of the support of the basis function relative to its center is characterized by some function $\rho =a_i\left(\psi \right)$, we can utilize ${\phi }_i={\left(1-\rho /a_i\left(\psi \right)\right)}_{+}$, where $w_{+}=w(\mathrm{sign}\left(w\right)+1)/2$. In this case, for the polygonal boundary, we obtain a piecewise linear function, for which $a_i\left(\psi \right)$ is calculated from the polar equations of the corresponding straight lines, i.e., at the corresponding intervals of the variation of $\psi$ we obtain $a_i\left(\psi \right)={\rho }_i/cos(\psi -{\psi }_i)$. If we want the basis function on the boundary to have a zero derivative, we can use ${\phi }_i={\left({\left(1-\rho /a_i\left(\psi \right)\right)}_{+}\right)}^2$. For a basis function of the form ${\phi }_i=\left(1+2\rho /a_i\left(\psi \right)\right){\left({\left(1-\rho /a_i\left(\psi \right)\right)}_{+}\right)}^2$, we obtain a smooth vertex. Some more complicated form of a basis function allows obtaining a smooth surface with a polygonal base, but we will not go into that.

In the multidimensional case, we can use $v_i=v_i(\rho ,\mathbf{v})$, where $\rho =\mathbf{x}-{\mathbf{c}}_i$ and $\mathbf{v}$ is a vector on the unit sphere, which can be parameterized by the corresponding coordinates.

2.3. Physics-Based Neural Network Approach with FEM Smoothing

Both physics-based models considered, the FEM and ML with physics-based regularization, have their advantages and disadvantages. Higher degrees of polynomial elements lead to an increase in the cost of computation, although we obtain a solution with higher smoothness. In any case, it is not possible to obtain a solution of infinite smoothness as for neural networks.The FEM is long and more extensively applied. It is implemented in numerous commercial software products. However, the FEM using piecewise linear functions leads to an insufficiently smooth solution. In addition, if the number of elements of the method is small, the accuracy is not high enough; if it is large, the system of linear equations may require large computational resources even with a sparse matrix of system (4).

The FEM is poorly suited to solving problems in which additional data are used to build a model. This is especially the case for tasks in which the model must adapt to the newly received data, taking into account changes in the modeling object. This is largely due to the local nature of the basis functions, which leads, in turn, to a local restructuring of the solution when trying to adapt it by minimizing the loss function $J_3$. The neural networks, by their nature, are capable of being used in the problem of adapting a solution to new data. This is evidenced by the extensive utilizing of neural networks to extract the information hidden in noisy data. However, the direct application of the above neural network approach often leads to a computationally complex problem of non-linear optimization. If only parameters $c_i$ of a neural network solution are adjusted, a system of linear equations with a dense matrix is obtained, which dramatically increases the complexity of its solution. A more efficient approach is to first optimize the loss function $J=J_1+\delta J_2$, and then train the neural network when new data arrive by minimizing $J_3$. However, this method is also very resource-intensive.

Alternatively, we propose the following hybrid approach. First, a solution to problem (1) is constructed using the FEM. Further, this solution is smoothed by means of a neural network. For the approach discussed above, all the basis functions except one are equal to zero at each nodal point, and this type of model is the easiest to build. It consists of optimizing $J_3$ in which the values of the FEM solution at the nodes are used instead of the measurement results. The neural network approximate solution obtained can be adapted to the newly received data in the same way as described above. Let us emphasize that in this hybrid model several physics-based modeling classes are combined.

In [32], a FEM solution is also used as training data in the problem of structural damage identification. However, in that work, the values at nodes are not used directly, but the loss function includes two terms with the corresponding weight factors. The first is responsible for a cross-entropy loss and the second corresponds to normalized damage probabilities.

2.4. Physics-Based Multilayer Network Architecture Method

The expressions of approximate solutions constructed by using the above methods are rather cumbersome. If we write them out analytically, they turn out to be poorly visible. Much more compact approximations are obtained through an approximate series expansion or various kinds of asymptotic methods [33]. However, these methods make it possible to build an acceptable solution only in a sufficiently small neighborhood of the point at which the expansion is constructed or for sufficiently small values of the parameter over which the expansion is carried out. In [24] we propose another approach, which has an essentially broader area of applicability. From the point of view of the modern modeling classification, we can attribute this method to physics-based network architectures. The resulting analytical models are fairly compact.

Let us consider the initial value problem

$$\left\{\begin{array}{c}{y}^{\left(n\right)}\left(x\right)=f(x,y\left(x\right),y^{\prime }\left(x\right),\dots ,y^{(n-1)}\left(x\right)),\hfill \\ y\left(x_0\right)=y_0,\hfill \\ y^{\prime }\left(x_0\right)=y_0^{\prime },\hfill \\ \dots \hfill \\ y^{(n-1)}\left(x_0\right)=y_0^{(n-1)},\hfill \end{array}x\in [x_0,x_0+a].\right.$$

(11)

Using the replacement, one can go from a higher-order system to a first-order multi-dimensional problem

$$\left\{\begin{array}{c}{\mathbf{y}}^{\prime }\left(x\right)=\mathbf{f}(x,\mathbf{y}\left(x\right)),\hfill \\ \mathbf{y}\left(x_0\right)={\mathbf{y}}_0,\hfill \end{array}x\in [x_0,x_0+a].\right.$$

(12)

To solve (12), many numerical methods have been developed [34]. A significant part of them consists of the partition of a given interval by points $x_k$ into intervals of length $h_k,$$k=1,\dots ,n,$ and applying the recurrence formula

$${\mathbf{y}}_{k+1}={\mathbf{y}}_k+F(\mathbf{f},h_k,x_k,{\mathbf{y}}_k,{\mathbf{y}}_{k+1}),$$

(13)

where operator F defines a specific method. If in (13) the function F depends on ${\mathbf{y}}_{k+1}$, it is necessary to accurately or approximately solve (13) with respect to ${\mathbf{y}}_{k+1}$.

Applying formula (13) iteratively n times to an interval with a variable right end $[x_0,x]\subset [x_0,x_0+a]$ we obtain function $\mathbf{y}\left(x\right)$, which can be considered an approximate solution of system (12). Let us underline that, in this case, $h_k=h_k\left(x\right),$${\mathbf{y}}_0={\mathbf{y}}_0\left(x\right),$${\mathbf{y}}_k={\mathbf{y}}_k\left(x\right).$

For constructing multilayer solutions, we can choose one of the classical numerical methods [34], for example:

• the implicit Euler method, where

  $$F(\mathbf{f},h_k,x_k,{\mathbf{y}}_k,{\mathbf{y}}_{k+1})=h_k\mathbf{f}(x_{k+1},{\mathbf{y}}_{k+1});$$

  (14)
• the trapezium method, where

  $$F(\mathbf{f},h_k,x_k,{\mathbf{y}}_k,{\mathbf{y}}_{k+1})=\frac{1}{2}{h}_k(\mathbf{f}(x_k,{\mathbf{y}}_k)+\mathbf{f}(x_{k+1},{\mathbf{y}}_{k+1}));$$

  (15)
• the midpoint method, where

  $$\left\{\begin{array}{cc}{\mathbf{y}}_{k+2}=F(\mathbf{f},h,x_k,{\mathbf{y}}_k,{\mathbf{y}}_{k+1}),\hfill & \\ F(\mathbf{f},h,x_{k+1},{\mathbf{y}}_k,{\mathbf{y}}_{k+1})={\mathbf{y}}_k+2h\mathbf{f}(x_{k+1},{\mathbf{y}}_{k+1}).\hfill \end{array}\right.$$

  (16)

or one of many other methods of this kind.

If function $\mathbf{f}(x,\mathbf{y})$ is a neural network (for example, it was obtained by modeling some related process), the result is an arbitrarily accurate solution to the original problem. In this case, we refer the entire modeling procedure to the hybrid methods mentioned in the Introduction as the embedding approach. Otherwise, it can be approximated by a neural network (for sufficiently wide classes of functions, this can be done with arbitrarily high accuracy) and the equation can be solved with the replacement of the right-hand side by its approximation.

The resulting models can be regarded as physics-based pre-trained multilayer networks, the numerical parameters of which can be further trained using measurement data.

Let us consider the boundary value problem

$$\left\{\begin{array}{c}{\mathbf{y}}^{\prime }\left(x\right)=\mathbf{f}(x,\mathbf{y}\left(x\right)),\hfill \\ \mathbf{u}\left(x_0\right)={\mathbf{u}}_0,\mathbf{v}(x_0+a)={\mathbf{v}}_0,\hfill \end{array}x\in [x_0,x_0+a].\right.$$

(17)

Here, vectors $\mathbf{v}$ and $\mathbf{u}$ are composed of the coordinates of vector $\mathbf{y}$; their total dimension is equal to the dimension of $\mathbf{y}$.

The following approach turned out to be an effective way to solve this problem. We choose arbitrary point $t\in (x_0,x_0+a)$ and build a solution to the system

$$\left\{\begin{array}{c}{\mathbf{y}}^{\prime }\left(x\right)=\mathbf{f}(x,\mathbf{y}\left(x\right)),\hfill \\ \mathbf{y}\left(t\right)=\mu ,\hfill \end{array}x\in [x_0,x_0+a].\right.$$

(18)

where parameter $\mu$ is adjusted based on the boundary conditions

$$\left\{\begin{array}{c}\mathbf{u}(x+0)={\mathbf{u}}_0,\hfill \\ \mathbf{v}(x_0+a)={\mathbf{v}}_0.\hfill \end{array}\right.$$

(19)

Parameter t remains undefined. It can be chosen so that a solution satisfies the original system with the highest accuracy. For the problems of modeling real objects and constructing DT, it is more natural to adjust t based on the best fit to the measurement data.

These multilayer models can easily adapt to new data about an object. For this purpose, we can change t and $\mu$ or only $\mu$ and minimize the loss function ${\sum }_{l=1}^L\left|\right|{\mathbf{y}}_n\left({\xi }_l^{″}\right)-{\mathbf{y}}_l{\left|\right|}^2$. $\left|\right|\ast \left|\right|$ takes into account only those coordinates of ${\mathbf{y}}_l$ for which measurement data are available. We have also studied the adaptive properties of the multilayer model in [35].

In [36], another scheme for imposing physics constraints on the neural network architecture is applied. They introduce a method that converts different types of equations into linear constraints and utilize them to build the output layer of a neural network. For inequality conditions, [36] suggests using positive definite activation functions.

3. Materials and Methods. Parameterized Problems

The formalization of parameterized problems is that some parameter $\mathbf{r}$ appears in the formulation of problem (1) in that $A=A\left(\mathbf{r}\right),$$B=B\left(\mathbf{r}\right)$, $\Gamma =\Gamma \left(\mathbf{r}\right)$, $\Omega =\Omega \left(\mathbf{r}\right)$. It is required to find a solution $u\left(\mathbf{r}\right)$ for which

$$\left\{\begin{array}{c}A\left(\mathbf{r}\right)u=g,u=u\left(\mathbf{x}\right),\mathbf{x}\in \Omega \left(\mathbf{r}\right)\subset {\mathbb{R}}^p;\hfill \\ {B\left(\mathbf{r}\right)u|}_{\Gamma }\left(\mathbf{r}\right)=f;\hfill \end{array}\right.$$

(20)

holds for all $\mathbf{r}$ from a given set.

3.1. FEM and Physics-Based Neural Network Approach

Direct application of the FEM to solving problem (20) leads to replacing (2) by an expansion

$$u_N(\mathbf{x},\mathbf{r})=\sum _{i=1}^N{c}_i\left(\mathbf{r}\right)v_i(\mathbf{x},\mathbf{r}),$$

(21)

and system (4) is replaced by system

$$\sum _{i=1}^N{c}_i\left(\mathbf{r}\right)A\left(\mathbf{r}\right)v_i(\mathbf{x},\mathbf{r})·v_j(\mathbf{x},\mathbf{r})=(g(\mathbf{x},\mathbf{r})-A\left(\mathbf{r}\right)v_0(\mathbf{x},\mathbf{r}))·v_j(\mathbf{x},\mathbf{r}).$$

(22)

As a result, we have a system of linear equations with coefficients depending on the parameter. There are three possibilities among physics-based methods to solve it. The first one is analytical. The second approach uses interval analysis. These techniques work in a reasonable amount of time only when the number of items is small. Thus, acceptable accuracy is only achieved for simple tasks. Thirdly, if the problem is solved for a sufficiently representative set of parameters, this can allow estimating a solution for the entire set of parameters. At the same time, computational complexity increases many times over, and we cannot be sure that the chosen set of parameters is sufficiently representative.

The neural network approach can also be applied to a problem with parameters. In this case, the parameters are included in the number of inputs of the neural network [29]. We have successfully solved a number of similar tasks [15,16]. A modification of the neural network approach to problem (20) consists of finding an approximate solution in the form of the output of an artificial neural network of a given architecture as sum

$$u_N(\mathbf{x},\mathbf{r})=\sum _{i=1}^N{c}_i{v}_i(\mathbf{x},\mathbf{r},{\mathbf{a}}_i),$$

(23)

where weights $c_i,{\mathbf{a}}_i$ are determined in the process of step-by-step network training based on minimizing the loss function of the form $J=J_1+\delta J_2$. Here,

$$J_1=\sum _{j=1}^M{(A\left({\mathbf{r}}_j\right)u_N({\xi }_j,{\mathbf{r}}_j)-g({\xi }_j,{\mathbf{r}}_j))}^2$$

(24)

and

$$J_2=\sum _{k=1}^K{(B\left({\mathbf{r}}_k^{\prime }\right)u_N({\xi }_k^{\prime },{\mathbf{r}}_k^{\prime })-f({\xi }_k^{\prime },{\mathbf{r}}_k^{\prime }))}^2$$

(25)

Regeneration of points is performed as follows: first, ${\mathbf{r}}_j$ and ${\mathbf{r}}_k^{\prime }$ are re-generated from a parameter change domain and then ${\xi }_j\in \Omega =\Omega \left({\mathbf{r}}_j\right)$ and ${\xi }_j^{\prime }\in \Gamma =\Gamma \left({\mathbf{r}}_j\right)$.

If additional data are available, for example, measurements, $J_3={\sum }_{l=1}^L{(u_N({\xi }_l^{″},{\mathbf{r}}_l^{″})-u_l)}^2$ is added to the loss function. Here, $u_l$ is a measurement corresponding to the value of a required solution at point ${\xi }_l^{″}$ and parameter value ${\mathbf{r}}_l^{″}$.

Despite the natural ability of neural network models to adapt to data, the process of constructing them often has a high computational complexity.

3.2. Parametric Physics-Based Neural Network Approach with FEM Smoothing

In order to combine the advantages of the FEM and neural network method for problems with a parameter, we propose the following hybrid method. Suppose we have constructed a finite element solution for a set of parameter values

$$u_N(\mathbf{x},{\mathbf{r}}_j)=\sum _{i=1}^N{c}_{i,j}{v}_i(\mathbf{x},{\mathbf{r}}_j),j=1,\dots ,M.$$

(26)

Further, we construct neural network approximations for the coefficients $c_i$

$$c_i\left(\mathbf{r}\right)=\sum _{k=1}^n{d}_{k,i}{\phi }_k(\mathbf{r},{\mathbf{a}}_{k,i}),$$

(27)

where ${\phi }_k(\mathbf{r},\mathbf{a})$ are neural network basis functions. We can apply the previously discussed RBF or perceptron-type functions. Our experience [10,29] shows that the latter perform significantly better in these tasks. Weights (parameters) of all neural networks are selected by minimizing corresponding loss functions, for example, in the form ${\sum }_{j=1}^M{(c_{i,j}-{\sum }_{k=1}^n{d}_{k,i}{\phi }_k({\mathbf{r}}_j,{\mathbf{a}}_{k,i}))}^2.$ The obtained solution

$$u_N(\mathbf{x},\mathbf{r})=\sum _{i=1}^N\sum _{k=1}^n{d}_{k,i}{\phi }_k(\mathbf{r},{\mathbf{a}}_{k,i})v_i(\mathbf{x},\mathbf{r})$$

(28)

can be utilized not only for interpolation in the range of parameters for which the FEM was used. It is also applicable for extrapolation to the region in which the FEM solution is difficult due to the ill-conditioned matrix.

3.3. Parametric Physics-Based Multilayer Network Architecture Method

As mentioned above, in the multilayer method physics defines network architecture. This technique can be applied to solving parameterized problems virtually without any modifications.

Let us consider the Cauchy problem of the form

$$\left\{\begin{array}{c}{\mathbf{y}}^{\prime }\left(x\right)=\mathbf{f}(x,\mathbf{y}\left(x\right),\mathbf{r});\hfill \\ \mathbf{y}\left(t\left(\mathbf{r}\right)\right)={\mathbf{y}}_0\left(\mathbf{r}\right);\hfill \end{array}\right.$$

(29)

and apply the same recurrence formulas to it as for the problem without a parameter. As a result, we have an approximate solution ${\mathbf{y}}_n(x,\mathbf{r})$. A separate issue is the selection of dependencies $t\left(\mathbf{r}\right)$. In simple problems we can use, for example, a linear dependence on $\mathbf{r}$. In tasks with measurements, we can select this dependence in the form of a neural network. In this case, $t\left(\mathbf{r}\right)={\sum }_{k=1}^n{p}_k{\phi }_k(\mathbf{r},{\mathbf{a}}_k)$ and ${\mathbf{y}}_0\left(\mathbf{r}\right)={\sum }_{k=1}^n{\mathbf{q}}_k{\phi }_k(\mathbf{r},{\mathbf{b}}_k)$.

At the first stage, we again write down the formula for approximate solution ${\mathbf{y}}_n(x,\mathbf{r},t\left(\mathbf{r}\right),{\mathbf{y}}_0\left(\mathbf{r}\right))$. At the second stage, we adjust parameters $p_k,{\mathbf{a}}_k,{\mathbf{q}}_k$ and ${\mathbf{b}}_k$ by minimizing the loss function ${\sum }_{l=1}^L\left|\right|{\mathbf{y}}_n({\xi }_l^{″},{\mathbf{r}}_l,t\left({\mathbf{r}}_l\right),{\mathbf{y}}_0\left({\mathbf{r}}_l\right))-{\mathbf{y}}_l{\left|\right|}^2$.

The study of the ability of a multilayer model to learn upon receipt of additional information about an object has been started in [37]. In [8], an arbitrarily accurate implicit Runge–Kutta time stepping scheme with an unlimited number of stages is utilized as the basis for one of two algorithms for constructing data-driven models based on physical constraints described by general non-linear partial differential equations. In addition, the problem of the data-driven discovery of partial differential equations is solved. However, the authors acknowledge the complexity of the proposed algorithms in comparison with classical numerical methods.

4. Results

We have studied the application of methods described above for solving the benchmark parameterized boundary value problem [33]

$$\left\{\begin{array}{c}{u}^{″}\left(x\right)+a\left(s\right)u^{\prime }\left(x\right)+b\left(s\right)u\left(x\right)=0;\hfill \\ u\left(0\right)=0,u\left(1\right)=1;\hfill \end{array}x\in [0,1].\right.$$

(30)

with coefficients $a\left(s\right)=1/s+s$ and $b\left(s\right)=1/s-s$.

For small values of s, problem (30) is stiff and can be regarded as a parameterized singular perturbation problem with a boundary layer in the vicinity of zero. This problem has an analytical solution that allows comparing the results of different methods clearly.

4.1. Physics-Based Neural Network Approach

Let us apply the neural network technique to solving Problem (30). The task is solved for fixed values of parameter s; a solution also has the form (9). The weights of the neural network are selected by minimizing the loss function

$$\sum _{i=1}^m{(u^{″}\left(x_i\right)+a\left(s\right)u^{\prime }\left(x_i\right)+b\left(s\right)u\left(x_i\right))}^2+{\delta }_1{u}^2\left(0\right)+{\delta }_2{(u\left(1\right)-1)}^2.$$

(31)

Test points $\left\{x_i\right\}$ are randomly chosen on the interval $[0,1]$ with regeneration every 5 steps of the optimization process. RProp is used as a minimization algorithm [38]. It is known [39,40] that neural networks with one hidden layer are universal approximators. With regard to the solution accuracy control, there are two different options. The first one is for setting the threshold value for the optimized loss function and the second one is for setting the limit on the number of network learning epochs. In the second case, among the results obtained, the one for which the loss function has the smallest value is selected. In addition, the residual based on real data about the object can be taken into account.

To solve the problem, we use two types of basis functions, sigmoids and Gaussians, and a different number of hidden layer neurons. For each neural network solution the plot of training loss history is presented. One epoch corresponds to 100 training point regenerations. Between each two regenerations, 5 steps of the RProp algorithm are conducted. The number of training points at each step equals 100. Figure 1 represents a solution to (30) for fixed $s=0.1$ in the case of a Gaussian neural network with 5 neurons in the hidden layer. The correspondent error is shown in Figure 2. The analytical solution has the form

$$\begin{array}{cc}\hfill & {\mathtt{nn}}_5\left(x\right)=0.6934e^{-1.0151{(-0.9445+x)}^2}+0.5966e^{-30.8231{(0.0819+x)}^2}+\hfill \\ \hfill & 0.4256e^{-63.199{(0.1164+x)}^2}+2.1744e^{-1.3991{(0.1815+x)}^2}-46.4433e^{-8.6259{(0.5628+x)}^2}.\hfill \end{array}$$

(32)

The training loss history is shown in Figure 3.

For sigmoid basis functions, the solution, and error plots are shown in Figure 4 and Figure 5. The analytical solution, in this case, looks like

$$\begin{array}{cc}\hfill & {\mathtt{n}}_5\left(x\right)=2.9401\mathtt{tanh}\left(0.1548(-1291.23+x)\right)-0.2423\mathtt{tanh}\left(2.8222(-0.5127+x)\right)\hfill \\ \hfill & -3.4427\mathtt{tanh}\left(0.2677\right(-0.4512+x\left)\right)+0.8904\mathtt{tanh}\left(6.6248\right(-0.0469+x\left)\right)\hfill \\ \hfill & +3.7422\mathtt{tanh}\left(9.4914\right(0.0895+x\left)\right).\hfill \end{array}$$

(33)

Figure 6 illustrates the training loss history for solution (33) as the decimal logarithm with respect to the training epoch.

For a problem in a stiff statement, for example, at $s=0.01$, without additional effort, it is possible to obtain an acceptable solution only using a network with 20 hidden layer neurons and sigmoid basis functions. The corresponding plots of the approximation and its error are shown in Figure 7 and Figure 8.

In this case, the neural network solution has a more cumbersome form. Figure 9 shows the training loss history.

In [41,42], physics-informed deep neural networks with adaptive activation functions are presented. Activators allow increasing the learning rate of neural networks. We have considered a network with one hidden layer and a global adaptive Gaussian activation function. However, the dynamic change in the loss function topology has not been observed.

The approximation of a parametric solution by a neural network requires a long training process. For a small number of layers, a result with a large approximation error can be obtained.

Based on the results for various basis functions in the case of the fixed-parameter problem statement, we have chosen a perceptron network with 100 hidden layer neurons and sigmoid activation functions to build a parametric neural network solution. Test points $\{x_i,s_i\}$ are randomly chosen on the domain $[0,1]\times [0.03,0.1]$. RProp is used as a minimization algorithm [38]. One epoch corresponds to 10 training point regenerations. Between each two regenerations, 5 steps of the RProp algorithm are conducted. The number of training points at each step equals 200.

Figure 10 illustrates the training loss history.

A 3D plot of the approximation error is shown in Figure 11.

4.2. Physics-Based Neural Network Approach with FEM Smoothing

To solve the parameterized problem, we use a neural network approach to smooth the FEM solutions constructed for fixed values of the parameter. The peculiarities of our method require utilizing a single set of FEM basis functions regardless of the parameter value. Thus, we do not use a FEM with adaptive grids (adaptive basis functions).

To start, we have solved system (30) by using the FEM with different values of a number K of equidistant vertexes $x_k$ and various values of a parameter s. For simplicity, instead of $u\left(x\right)$, we have sought an approximate solution to the problem with homogeneous boundary conditions, i.e., an approximation for $u\left(x\right)-x$.

It is the one-dimensional case, so an approximate solution to this system takes the form

$$u_N\left(x\right)=\mathbf{c}\left(s\right)·\mathbf{u},$$

(34)

where vector $\mathbf{u}$ consists of K piecewise linear functions $u_k\left(x\right)$

$$u_k\left(x\right)=\left\{\begin{array}{c}\frac{x-x_{k-1}}{x_k-x_{k-1}},x\in [x_{k-1},x_k];\hfill \\ \frac{x_{k+1}-x}{x_{k+1}-x_k},x\in [x_k,x_{k+1}];\hfill \\ 0,\mathrm{otherwise};\hfill \end{array}\right.$$

(35)

the coefficients $\mathbf{c}\left(s\right)=\left\{c_k\left(s\right)\right\}$ are obtained by solving the corresponding systems of linear equations.

Thus, we have calculated the values of coefficients $c_i\left(s_l\right)$ of FEM solutions (34) for $s_l=0.05,\dots ,0.1$ with step 0.01 and $K=10$ nodes. In the end, we have obtained training sets $\{s_l,c_i\left(s_l\right)\}$ which are utilized for the neural network approximations of FEM coefficients depending on parameter s

$$p_i\left(s\right)=a_{1i}+a_{2i}\mathtt{tanh}\left(a_{3i}(s-a_{4i})\right).$$

(36)

After approximating each of the coefficients, we have a general parametric solution

$$\mathtt{znn}(s,x)=\sum _{i=1}^{K-1}{p}_i\left(s\right){\phi }_i\left(x\right).$$

(37)

This solution, in turn, can be considered a neural network output with piecewise-specified activation functions

$${\phi }_i\left(x\right)=\left\{\begin{array}{c}1-|K(x-\frac{i}{K})|,|x-\frac{i}{K}|<\frac{1}{K};\hfill \\ 0,\mathrm{otherwise}.\hfill \end{array}\right.$$

(38)

The neural network (37) is a hybrid version of an RBF network with preselected centers and widths and a perceptron with one hidden layer.

Figure 12 represents neural network extrapolations of the solution to problem (30) for parameter s value 0.01, which is not included in the training set. Parametric solution (37) obtained in this way shows a fairly good approximation of a solution to problem (30) for small values of parameter s. In order to compare, Figure 12 shows the FEM solution with the same number $K=20$ of non-adaptive basis functions (35), which reflects the stiffness of problem (30) at small values of the parameter. In addition, simple smoothing of this FEM solution, from which the term x is subtracted, by a network with 3 hidden layer neurons

$$zn\left(x,\mathbf{p},\mathbf{q},\mathbf{r}\right)=p_0+\sum _{i=1}^3{p}_i\mathtt{tanh}\left(q_i\left(x-r_i\right)\right),$$

(39)

also improves the quality of the solution (see Figure 12, red and orange lines).

The parameters $\mathbf{p}$, $\mathbf{q}$, $\mathbf{r}$ are selected by optimizing the standard loss function using the RProp method.

4.3. Physics-Based Multilayer Network Architecture Method

To apply an idea described in subsection Physics-Based Multilayer Network Architecture Method, we represent a solution of system (30) as a linear combination of two boundary problem solutions $\phi \left(x\right)=({\phi }_1\left(x\right),{\phi }_2\left(x\right))$ and $\xi \left(x\right)=({\xi }_1\left(x\right),{\xi }_2\left(x\right))$. The first problem is

$$\left\{\begin{array}{c}{y}^{\prime }\left(x\right)=z\left(x\right);\hfill \\ z^{\prime }\left(x\right)=-a\left(s\right)z\left(x\right)-b\left(s\right)y\left(x\right),x\in [0,1];\hfill \\ y\left(x_0\right)=1;\hfill \\ z\left(x_0\right)=0.\hfill \end{array}\right.$$

(40)

and the second one is

$$\left\{\begin{array}{c}{y}^{\prime }\left(x\right)=z\left(x\right);\hfill \\ z^{\prime }\left(x\right)=-a\left(s\right)z\left(x\right)-b\left(s\right)y\left(x\right),x\in [0,1];\hfill \\ y\left(x_0\right)=0;\hfill \\ z\left(x_0\right)=1.\hfill \end{array}\right.$$

(41)

Here, optimal point $x_0$ depends on a parameter s and can be selected by a researcher. Throughout this article, an optimal point is considered a point minimizing the maximum error of the approximate solution on the segment $[0,1]$.

Solutions to Problems (40) and (41) can be obtained by means of the multilayer method. Then a solution of system (30) is

$$u\left(x\right)=c_1{\phi }_1\left(x\right)+c_2{\xi }_1\left(x\right).$$

(42)

Coefficients $c_{1,2}$ are determined by the boundary conditions of system (30)

$$\left\{\begin{array}{c}{c}_1{\phi }_1\left(0\right)+c_2{\xi }_1\left(0\right)=0;\hfill \\ c_1{\phi }_1\left(1\right)+c_2{\xi }_1\left(1\right)=1.\hfill \end{array}\right.$$

(43)

To solve problems (40) and (41), different types of function F (see (13)) have been considered. In the case of using the implicit Euler method (14) as a basis for a multilayer solution, systems (40) and (41) are transformed to

$$\left\{\begin{array}{c}{y}_k^{\prime }\left(x\right)=y_{k-1}+h{z}_k;\hfill \\ z_k=z_{k-1}-h(a\left(s\right)z_k+b\left(s\right)y_k).\hfill \end{array}\right.$$

(44)

This system can be solved with respect to $y_k$ and $z_k$. As a starting point $x_0$, we considered 0 and the optimal point. The resulting multilayer solution to problem (30) for $n=10$ layers and optimal point $x_0=0.151$ has the analytical form

$$y_n(x,0.151,10)=\frac{58.9+303x+850x^2+1940x^4+1680x^5+100x^6+391x^7+90.6x^8+9.44x^9}{{10}^{-9}{x}^{-1}{(8.59+9.9x+x^2)}^{10}}.$$

(45)

Corresponding multilayer solutions based on the implicit Euler method for $x_0=0$ and $x_0=0.151$ at $s=0.1$ are shown in Figure 13.

If function F from (13) corresponds to the trapezoidal method, systems (40) and (41) are reduced to the form

$$\left\{\begin{array}{c}{y}_k=y_{k-1}+\frac{1}{2}h(z_k+z_{k-1});\hfill \\ z_k=z_{k-1}-\frac{1}{2}h(a\left(s\right)z_k+b\left(s\right)y_k+a\left(s\right)z_{k-1}+b\left(s\right)y_{k-1}).\hfill \end{array}\right.$$

(46)

At the same time, for parameter $s=0.1$, a three-layer network gives a solution with an error less than 0.02 even without selecting the optimal point. This fact is illustrated in Figure 14. It also shows how the selection of the starting point has a significant effect on a two-layer network result.

Let us consider a more stiff formulation of the problem when the parameter s value is equal to 0.001. In this case, it is possible to obtain a sufficiently accurate approximate solution by using a two-layer network by selecting the optimal point (see Figure 15), and Formula (13) leads to a simple form of the solution

$$y_{\mathtt{tr}}(x,0.0046,2)=\frac{x(-0.0498+44.7x-22.5x^2+2.79x^3)}{{(0.00238-3.99x-x^2)}^2}.$$

(47)

Increasing the number of layers maintains the effect.

Let us solve parameterized problem (30) by a multilayer method. In this case, parameter s is naturally contained in a solution,

$$y(x,s)=c_1\left(s\right){\phi }_1(s,x)+c_2\left(s\right){\xi }_1(s,x).$$

(48)

Thus, a one-layer solution based on the trapezium method and for $x_0=1.7s$ has the form

$$y_1(x,s)=\frac{x(-1.04+0.969s-1.35s^2+s^4)(0.059x+s^3-0.56s^{2}x-s-1.39){(s^3-0.59s^2-s-0.796)}^{-1}}{s^4+1.18(1-x)s^3-(0.692x+0.346x^2)+s(1.18x-0.208)+s^2(0.346x^2-0.692x-1)}$$

(49)

In this case, the point $x_0$, as mentioned earlier, depends on s but should not be beyond the considered interval.

Despite problem (30) being stiff for small s, we obtain a fairly good solution for all values of s from the interval $[0.001,0.1]$ using the multilayer technique based on the trapezium method with 8 layers. This is confirmed by the 3D plot of the approximation error in Figure 16.

Figure 17 illustrates that choosing the required number of layers allows us to obtain an approximation of any accuracy.

The above results of the multilayer method in solving a parameterized problem singularly perturbed by a small parameter make it possible to consider this method as promising in the modeling and control of complex systems. It is important to note that, if the problem of constructing a preliminary model of a real object is considered, a fairly good approximation for small values of the parameter gives a solution with a single layer. In addition, the known estimates of approximation errors of the basic recurrence relations can be used to define the number of layers to obtain a solution with a preselected accuracy. The construction of a multilayer solution takes little time and can immediately be used for different values of the parameter. Unlike various neural network approaches, a multilayer model does not require training. At the same time, the subsequent adjustment of the parameters of such a solution using the apparatus of neural networks is not excluded in the case of additional data appearing. On this occasion, the general neural network approach should show itself well, which is confirmed by both the results of our previous article [15] and the hybrid neural network models constructed in this work using the “poor” results of the classical FEM.

In the context of new challenges facing the modern control theory and, in particular, the methodology for creating CPS, problems are often ill-posed and poorly formalized, and are set using various heterogeneous data, parameterized differential equations of unclear structure, or without them. Requirements for the speed of processing such data and solving dynamic equations are growing. Thus, the development of versatile low-cost methods presented in this work is an important challenge.

5. Discussion

For modeling CPS, an urgent issue is the transition from differential models to functional ones. At the same time, an important factor is a possibility to adapt these models according to the observation data of the object. Our computational experiments allow drawing preliminary conclusions about what kind of models may be more suitable for further adaptation.

The main advantage of finite element models is that they are well studied. Accordingly, there are a large number of software packages and libraries. At the same time, with the appearance of stiffness in the problem being solved, the number of necessary functions to obtain acceptable accuracy increases sharply. Moreover, difficulties in adapting the solution increase, since it becomes unstable against data errors. This problem is aggravated by the local nature of the basis functions. Neural network models are more capable of adaptation. However, the process of their construction has a high computational complexity that also greatly increases with the appearance of stiffness. This situation allows us to look prospectively at the hybrid approaches we have considered. Even more interesting are multilayer models, which are the most compact for a given accuracy. We have already studied their adaptive properties on one problem [37]. We intend to continue this study in further research.

Another important issue is computational costs. The fastest are multi-layered methods as they do not require training. The next ones are methods with constant basis (activation) functions but, as a result of their application, larger models are obtained. It also takes more time to adapt them further. The training of neural networks can be accelerated, but this is relevant for the stage of adaptation from measurements in real time. In the future, we intend to compare the adaptive properties of multilayer models and neural networks with different acceleration options: structure adaptation, faster methods, and the introduction of hyperparameters.

The issue of the extension of the presented results to nonlinear problems is worth specific discussion. For these problems, the advantages of the FEM are substantially leveled, since the construction of a FEM solution is no longer reduced to solving a linear system. At the same time, the complexity of building a neural network and a multilayer solution changes only slightly.