Integration of Classical Mathematical Modeling with an Artificial Neural Network for the Problems with Limited Dataset

Abstract

One of the most common problems in science is to investigate a function describing a system. When the estimate is made based on a classical mathematical model (white-box), the function is obtained throughout solving a differential equation. Alternatively, the prediction can be made by an artificial neural network (black-box) based on trends found in past data. Both approaches have their advantages and disadvantages. Mathematical models were seen as more trustworthy as their prediction is based on the laws of physics expressed in the form of mathematical equations. However, the majority of existing mathematical models include different empirical parameters, and both approaches inherit inevitable experimental errors. Simultaneously, the approximation of neural networks can reproduce the solution exceptionally well if fed sufficient data. The difference is that an artificial neural network requires big data to build its accurate approximation, whereas a typical mathematical model needs several data points to estimate an empirical constant. Therefore, the common problem that developers meet is the inaccuracy of mathematical models and artificial neural networks. Another common challenge is the mathematical models’ computational complexity or lack of data for a sufficient precision of the artificial neural networks. Here we analyze a grey-box solution in which an artificial neural network predicts just a part of the mathematical model, and its weights are adjusted based on the mathematical model’s output using the evolutionary approach to avoid overfitting. The performance of the grey-box model is statistically compared to a Dense Neural Network on benchmarking functions. With the use of Shaffer procedure, it was shown that the grey-box approach performs exceptionally well when the overall complexity of a problem is properly distributed with the mathematical model and the Artificial Neural Network. The obtained calculation results indicate that such an approach could increase precision and limit the dataset required for learning. To show the applicability of the presented approach, it was employed in modeling of the electrochemical reaction in the Solid Oxide Fuel Cell’s anode. Implementation of a grey-box model improved the prediction in comparison to the typically used methodology.

1. Introduction

Predicting the output of a system is one of the most frequent tasks in the fields of research in engineering, economics, medicine, etc. It may be a performance estimation of a device, cost analysis, or predictive medicine. Whenever possible, classical mathematical models are developed based on the laws of physics, observations, and, unavoidably, assumptions. In this paper, the phrase “the classical mathematical model” is used as an equivalent of white-box model. Mathematical models often can be inaccurate, incomplete, or very hard to formulate due to gaps in the existing knowledge. In such cases, to be able to predict a system’s output, approximation methods are used. Computational models known as Artificial Neural Networks (ANN) have brought a vast improvement to predictions in many fields. An ANN can achieve excellent performance in function approximation, comparable to accurate mathematical models [1]. For such high accuracy, an ANN needs a lot of data. Obtaining enough samples may be expensive and time-consuming, and, in effect, unprofitable. Data numerousness is mostly dependent on the complexity of the process one wants to predict. As a vivid example of such a problem, we can offer a fuel cell modeling a multiscale and interdisciplinary problem. The fuel cell models are hard to generalize over different types and have very complicated models because of various species transport phenomena and electrochemical reaction kinetics [2,3]. Obtaining the necessary tomographic data for electrodes’ microstructure characterization costs a few months of operator work [4,5]. For such a problem, collecting more than twenty data points per dimension would take years and therefore is infeasible in many cases [6,7,8].

In this paper, we address the problem of data shortage by integrating the classical mathematical model with an artificial neural network as presented in Figure 1. For clarification, we will refer to our implementation of a grey box model as the Interactive Mathematical modeling-Artificial Neural Network (IMANN). The approach presented in this study stems from the fact that every mathematical model’s solution, including the example given above, is represented as a function or a set of functions. The method consists of determining the most uncertain parts of a mathematical model or lacking a theoretical description, and substituting them with a prediction of an ANN, instead of using the mathematical model or the ANN alone. With that being said, the work with the IMANN differs from the practice with a regular ANN. In a conventional procedure, a dataset is divided into training and test data. In addition to these steps, the IMANN requires to divide the problem into two parts. The mathematical model describes one of them, and the artificial neural network approximates the other. The decision regarding this division is a crucial part of working with IMANN, which affects the architecture of the network as well as the data. The underlying research of this paper is to improve the predictive accuracy under a limited dataset, which is understood as less than twenty points per input dimensionality. This research problem is addressed by incorporating an artificial neural network to predict different parts of benchmark functions, which are regarded as a general representation of mathematical models. The ANN is learned based on the mathematical model’s output. A detailed description is presented in Section 3. In practice, it might be understood that the ANN prediction replaces only some equations in a model (or even just a part of an equation) and adapts to the system’s behavior. As an example, let us consider a system of equations in which one of the equations is replaced by the prediction done by an artificial neural network. The obtained approximation values of that function would be forced by the numerical model to fulfill all equations in the system. If the prediction is incorrect there would be a discrepancy between the measured and the expected output during training. As a consequence, the ANN would be forced to improve its weights and biases until the model is in accordance to the experimental data. As we note later in this work, such a replacement can benefit the accuracy and minimal dataset needed for an artificial neural network prediction.

2. Literature Review

The problem with limited datasets is addressed frequently in the literature with many different approaches [9,10,11,12]. One type of method is a data augmenting—the generation of a slightly different sample by modifying existing ones [13]. Baird et al. used a combination of many graphical modifications to improve text recognition [14]. Simard et al. proposed an improvement: the Tangent Prop method [15]. In Tangent Prop, modified images were used to define the tangent vector—an imitation of function derivative, which was included in error estimation [15]. Methods in which such vectors were used are still improved in recent works [16,17]. In the literature, one can find a variety of methods that can modify datasets, improving ANN learning. A remarkably interesting approach, when dealing with two- and three-dimensional images is based on a persistent diagram (PD) technique. The PD changes the representation of the data to extract crucial characteristics [18] and uses as little information as possible to store them [19]. Adcock et al. in [20] presented how persistent homology can improve machine learning. All mentioned techniques are based on the idea to manipulate the dataset on which the ANN is trained.

Another method is to alter the ANN by including knowledge into its structure. Typical models are of black-box or white-box type. White-box models are robust and trustworthy due to the incorporation of physical laws. On the other hand, white-box models need empirical parameters and usually work only for specific conditions. Black box models need big datasets and can return unrealistic results, but they are great for generalization and, after training, their response time is negligible. A natural idea is to combine these two approaches into a grey model, which can result in great accuracy and generality. A successful attempt to add knowledge is made by a Knowledge-Based Artificial Neural Networks (KBANN) [21]. The KBANN starts with some initial logic, which is transformed into the ANN [21]. This ANN is then refined by using the standard backpropagation method [21]. The KBANN utilizes a knowledge which is given by a symbolic representation in the form of logical formulas [21]. In situations when knowledge is given by a functional representation, i.e., containing variables, an interesting approach was presented by Su et al. in [22]. The authors proposed a new type of neural network—an Integrated Neural Network (INN) [22]. In the INN, an ANN is coupled with a mathematical model in such a way that it learns how to bias the model’s output, to improve concurrence with the modeled system [22]. The INN output consists of the sum of the model and ANN output [22]. A similar approach to improving the mathematical model was presented by Wang and Zhang in [23]. They proposed an ANN, in which some of the neurons had their activation functions changed to the empirical functions [23]. The ANN was used to alter the empirical model of the existing device in such a way, that it can be used for a different device. This idea led to a Neuro-Space Mapping (Neuro-SM) in which the functional model is a part of the ANN [24,25]. Neuro-SM can be viewed as a model augmented by the ANN. In the Neuro-SM, the ANN maps the input and output of the model, but it does not interfere with the functional representation itself.

In more recent works, the most successful attempt for the interaction between physics and ANN is coupling ANNs with differential equations. The very first attempt to connect ANN with differential equations was performed by Psichogios and Ungar in [26]. The authors presented a way to improve the white-box model by estimating the empirical part of a differential equation describing a fed-batch bioreactor with an ANN. The work of Psichogios and Ungar was continued by Hagge et al. [27] with use of modern machine learning techniques on the same physical problem as Psichogios. The fed-batch bioreactor problem was further addressed in [28,29,30]. Many more researchers attempted to solve the problem in a similar manner—identify the difficult to model parts in first principles and substitute them with machine learning techniques: Cubillos et al. [31] in solid drying process, Piron et al. in crossflow microfiltration process [32], Oliveira in fed-batch bioreactor and bakers’ yeast production [30], Vieira and Mota in water gas heater system [33], Romijn et al. in energy balance for a glass melting process [34] and Chaffart and Ricardez-Sandoval in thin film growth process [35]. Cen et al. [36] presented a method for incorporating several neural networks into nonlinear dynamic systems for fault estimation problems. A different approach, which used differential equations as a basis for ANN was presented by Lagaris et al. [37]. The mathematical model solution was presented as a sum of a predefined function which fulfilled the boundary conditions and Neural Network prediction. Error was estimated by providing the ANN solution to the differential equation. An interesting approach called Physics-Informed Neural Network (PINN) was presented by Raissi et al. [38]. The authors presented a framework in which ANN was used to solve PDE system [38]. In PINNs, an ANN is used to learn how to solve given PDE system with use of sum of two errors—one consists of difference between imposed boundary conditions and second is a residual on arbitrary set of collocation points. The general framework was presented by Parish et al. in [39]. The idea is currently extensively studied as can be seen in [40,41,42] and used in complex applications [43,44]. An advanced approach in the field of Computational Fluid Dynamics (CFD) was proposed by Ling et al. [45] and later on by Wu et al. in [46]. Chan and Elsheikh [47] proposed ANN architecture for application in multiscale methods. An interesting utilization of ANN was presented by Tripathy and Bilionis in [48]. Authors used Neural Networks for high dimensional uncertainty quantification, which can be used for computation-heavy numerical simulations.

As presented in the literature survey, the grey-box models attract continuously increasing attention. Unlike other artificial intelligence methods, artificial neural networks that incorporate knowledge in their structure receive the technology’s continuous growth. The increasing focus on technology is due to the high utilitarian value of this approach. Including artificial neural networks to make a data-driven prediction for the most nonlinear part of the model significantly reduces the required dataset and improves accuracy. The need for a more computationally complex problem creates a need for developing new methods of integrating models and the network.

In this paper, we present a new approach, in which almost any part of a model, from a constant to an entire equation, can be replaced by ANN’s estimation. The interaction is implemented by shifting a part of the mathematical model to be predicted by the ANN and teaching the ANN with an evolutionary algorithm (EA) using a numerical model’s output errors. Learning with evolutionary algorithms allows flexibility and many broad applications, e.g., in standard backpropagation algorithm, an error is needed for every output, which is not possible when the ANN has to approximate a function for every domain point required for the model. The novelty of our approach lies in the interaction between the numerical model and the ANN, in which any intermediate variable or parameter from the numerical model can be an ANN’s input, and the ANN can be called multiple times during one model prediction, without affecting learning algorithm. We have tested our approach statistically with simple benchmarking functions as well as with a real-life application in Solid Oxide Fuel Cell modeling.

3. Methodology

3.1. Implementation

System boundaries are supplied to the ANN, which predicts the assigned part of the mathematical model. System boundaries alongside with an ANN’s output are provided to the numerical model. The numerical model computes the predicted system output. At the learning phase, the predicted output error is used to calculate the proper weights and biases of the ANN. When there is a necessity (e.g., a function is needed to obtain model prediction), IMANN is implemented in such a way that the numerical model sends a request for a part of the mathematical model to ANN with additional parameters (e.g., state of the model, coordinates etc.) Such an approach allows for an interaction between the numerical model and ANN. It should be noted that this is not limited to one request, i.e., the numerical model can request the ANN’s prediction multiple times while producing one output. The proposed approach to a prediction problem is graphically presented in Figure 2. The general approach to developing the IMANN network for any existing numerical model is performed as follows. Having a numerical model $M\left(\mathit{i}\right)$, where $\mathit{i}$ are the system boundaries and parameters, define a part of the model which will be substituted by an ANN prediction. In place of that part, insert a call to the ANN with input as $\mathit{y}$, where $\mathit{y}$ is a vector of all ANN’s input, which will contain selected variables from $\mathit{i}$ as well as from variables known by the numerical model during computation. Add ANN as an input to the numerical model. The described algorithm will be further addressed as $\mathrm{IMANN}(M,ANN)$. The selection of weights and biases is performed by an evolutionary algorithm, which was suggested in the literature for grey-box models [30] but not treated in much detail. Every individual is a representation of one network, i.e., its weights and biases are in the form of a vector. The fitness function is a measure of the discrepancy between a IMANN’s output and a modeled system.

3.2. Learning Process

To train the IMANN, the Covariance Matrix Adaptation-Evolution Strategy (CMA-ES) algorithm is chosen. The CMA-ES algorithm was introduced in [49]. It is an effective and flexible algorithm, which is ideally suited to IMANN learning, where the problem taken under consideration can have different degrees of dimensionality. Here, the dimensionality of the problem is the number of weights and biases in the IMANN. To train the IMANN, the CMA-ES uses a vector made from the network’s weights and biases. The CMA-ES optimizes the vector values based on the error obtained on the training data. Details about the objective function are discussed in the following section. Open-source implementation of CMA-ES in Python was adopted for evolutionary weights adjustments [50].

3.3. Learning Algorithm

1.

Having an empirical dataset D of a N system’s input ($\mathit{x}$) and output ($\Xi$) paired vectors, define a function which divides the dataset D into training and test data sets—$D_{\mathrm{train}}$ and $D_{\mathrm{test}}$.

2.

Having a numerical model M and neural network $NN$, define a function IMANN(M,$NN$) that combines the M and $NN$ into IMANN.

3.

For a given CMA-ES, optimization parameters $\left(p_{\mathrm{CMA}}\right)$ initialize the CMA-ES optimization algorithm.

4.

Perform one optimization step for optimization of the weights and biases vector $\mathit{w}$ based on a fitness function $E_{\mathrm{train}}$ dependent on the $D_{\mathrm{train}}$ and $\mathit{w}$.

5.

For every solution ${\mathit{w}}_i$ check squared error made on $D_{\mathrm{train}}$ ($E_{\mathrm{train}}$) and $D_{\mathrm{test}}$ ($E_{\mathrm{test}}$), save the lowest, throughout whole optimization, value of $E_{\mathrm{train}}^2+E_{\mathrm{test}}^2$ and the corresponding solution ${\mathit{w}}_i$.

6.

Repeat steps 4 and 5 until one of the optimization’s stopping criteria is satisfied.

The above steps in the form of pseudocode are presented in Algorithm 1. Using the squared values of training and test errors is a way to avoid overfitting.

Algorithm 1 IMANN training process

Require: Empirical dataset $D=\left\{({\mathit{x}}_i,{\Xi }_i)|0\le i\le N\right\}$, CMA-ES optimization parameters $p_{\mathrm{CMA}}$, numerical model M

1: $f_{\min }=\infty$

2: $\left(D_{\mathrm{train}},D_{\mathrm{test}}\right)=\mathrm{DivideData}\left(D\right)$ // Divide data to training and test sets

3: CMAoptimizer = InitializeCMA$\left(p_{\mathrm{CMA}}\right)$ // Initialize optimizer

4: do

5:  NextGeneration(CMAoptimizer)

6:  $P=\mathrm{GetPopulation}\left(\mathrm{CMAoptimizer}\right)$ // Get current population

7:  for $c_i\in P$ do

8:    $I=\mathrm{Network}\left(c_i\right)$ // Initialize network from $c_i$ candidate’s weights vector ${\mathit{w}}_i$

9:    $IM=\mathrm{IMANN}(M,I)$ // Initialize IMANN as an integration of the numerical model and the ANN

10:    $E_{\mathrm{train}}=Error(IM,D_{\mathrm{train}})$

11:    setCandidateFitness($c_i$, $E_{\mathrm{train}}$)

12:    $E_{\mathrm{test}}=Error(IM,D_{\mathrm{test}})$

13:    $E=E_{\mathrm{train}}^2+E_{\mathrm{test}}^2$

14:    if$f_{\min }>E$  then

15:    $f_{\min }\leftarrow E$

16:    save(${\mathit{w}}_i$)

17:    end if

18:  end for

19: while !stopCondition(CMAoptimizer)

The difference between the PINN and the IMANN training scheme is the usage of models output. In PINN’s learning, the physical model is a part of the neural network and, therefore, is treated as an activation function used to backpropagate the error. Additionally, fulfillment of the corresponding PDE is a part of the optimization. In the IMANN, the part of the solution proposed by the ANN is provided to the numerical model. In terms of fulfilling the PDE, IMANN relies on the numerical model to provide a physically correct prediction.

3.4. This Study

Choosing an EA as a primary optimization technique in the grey-box model comes from two factors. First, introduced earlier in this section, the ability to adjust weights based on a full function predicted by an ANN without the knowledge of the error in each point. Second, the BP algorithm uses derivatives to improve the prediction. Treating the model as an activation function in the last neuron can be poorly conditioned or impossible to solve in specific applications. As a simple example, estimating a constant a in the equation $sin\left(ax\right)$ is used. We use two datasets consisting of 20 and 100 points. In both, every fifth point is assigned to the test set (overall 20%). Estimation is performed with the EA and BP algorithms for 500 attempts. The error estimation procedure is described at the beginning of Section 5. Histograms are presented in Figure 3 and Figure 4. For both datasets, the distribution of errors made by the EA’s approach is smaller. The EA and BP approaches struggle with the ambiguous nature of the function - for the finite amount of uniformly distributed data points generated using a trigonometric function, one can find a different function of the same type but with a higher frequency that will have the same function values in provided data points. The BP algorithm additionally has a problem with the last activation function derivative, which is equal to $xcos\left(ax\right)$. Such derivative stops the BP algorithm in a local minimum, if the initial weights are not chosen close enough to a good solution. Since the ambiguity of many real problems can lead to unreal models, the parameters proposed by the IMANN should be monitored. In real applications, physical constraints should be added, such as boundaries for the parameters imposed on the ANN output or in the form of additional error terms during learning.

4. Benchmarking Functions

Every system’s output, physical or theoretical, depends on its boundaries (inputs) and characteristics (adaptable parameters). Mathematical models are functions that try to reflect real system behavior. Mathematical models are often based on assumptions and empirical parameters due to the gap in existing knowledge regarding the phenomena. The simplifications in problem formulation results in the discrepancy between the model and the real system outputs. The difference can dissolve only for hypothetical cases where the system output is described by mathematical equations and well-defined.

Every function can be treated as a system, and any arbitrary function can be treated as its mathematical model. The concept of IMANN strives to be reliable and applicable to any system: physical, economic, biological, or social, only when a part of this system can be represented in a mathematical form. To be able to represent a wide range of possible applications, benchmark functions were employed as a representation of the system. The system’s inputs ($\mathit{x}$) and measurable outputs ($\Xi \left(\mathit{x}\right)$) correspond to the inputs and outputs of benchmarking functions. The ANN predicts a part of the system, and the numerical model calculates the rest of it. From the model’s view, the ANN’s prediction is a parameter or one of the model’s equations. Values calculated by the benchmarking functions represent the measurable system outputs. If the subfunction values predicted by the ANN had the same value as calculated from the extracted part of the benchmark function, this would represent a perfect match between the IMANN and the system output. If the real system output differs from the predicted one, the IMANN will be forced to improve the weights and biases.

Functions

To test the IMANN, two benchmarking functions are used. Typical test functions used in optimization problems were used. One polynomial function—Styblinski-Tang function—for a one-dimensional input and Rosenbrock function for two dimensions. The chosen polynomial function is given by formula:

$$f_{\mathrm{P}}\left(x\right)=x^5-16x^3+5x^2,$$

(1)

the two-dimensional Rosenbrock function is given by:

$$f_{\mathrm{R}}(x,y)=100{\left(y-x^2\right)}^2+{\left(1-x\right)}^2.$$

(2)

In the case of polynomial functions, eight formulations of mathematical models are used, four with one subfunction and four with two subfunctions. The difference between the formulations lies only in the form and number of the subfunctions. Model formulations with one subfunction are given by the following equations:

$$f_1\left(x\right)=a\left(x\right)x^5-16x^3+5x^2,$$

(3a)

$$f_2\left(x\right)=\widehat{a}\left(x\right)x^4-16x^3+5x^2,$$

(3b)

$$f_3\left(x\right)=\tilde{a}\left(x\right)x^3-16x^3+5x^2,$$

(3c)

$$f_4\left(x\right)=\overline{a}\left(x\right)-16x^3+5x^2,$$

(3d)

where $a,\widehat{a},\tilde{a}$ and $\overline{a}$ are subfunctions, and $f_i$ is the i-th model function. We have chosen the formulations to distribute the complexity of prediction between ANN and numerical model in different scenarios. Equation (3a) is a case where ANN is responsible for a very simple task, while Equation (3d) results in a complexity similar to the predicting a value of Equation (1). For a perfect match with the modeled benchmarking function, the subfunctions should be functions of x in the form:

$$a\left(x\right)=1$$

(4a)

$$\widehat{a}\left(x\right)=x$$

(4b)

$$\tilde{a}\left(x\right)=x^2$$

(4c)

$$\overline{a}\left(x\right)=x^5.$$

(4d)

The polynomial function’s model formulations with two subfunctions are defined as:

$$f_5\left(x\right)=a\left(x\right)x^5+b\left(x\right)x^3+5x^2,$$

(5a)

$$f_6\left(x\right)=\widehat{a}\left(x\right)x^4+\widehat{b}\left(x\right)x^2+5x^2,$$

(5b)

$$f_7\left(x\right)=\tilde{a}\left(x\right)x^3+\tilde{b}\left(x\right)x+5x^2,$$

(5c)

$$f_8\left(x\right)=\overline{a}\left(x\right)+\overline{b}\left(x\right)+5x^2,$$

(5d)

where $a,\widehat{a},\tilde{a}$, $\overline{a},b,\widehat{b},\tilde{b}$ and $\overline{b}$ are subfunctions. The idea behind choosing different formulations is the same as for Equation (3a–d). All model functions are defined in the domain $\Omega =\left\{x:x\in [-4,4]\right\}$. Ideally, the subfunctions should be functions of x in the form:

$$\begin{array}{cc}a\left(x\right)=1,\hfill & \hfill b\left(x\right)=-16,\end{array}$$

(6a)

$$\begin{array}{cc}\widehat{a}\left(x\right)=x,\hfill & \hfill \widehat{b}\left(x\right)=-16x,\end{array}$$

(6b)

$$\begin{array}{cc}\tilde{a}\left(x\right)=x^2,\hfill & \hfill \tilde{b}\left(x\right)=-16x^2,\end{array}$$

(6c)

$$\begin{array}{cc}\overline{a}\left(x\right)=x^5,\hfill & \hfill \overline{b}\left(x\right)=-16x^3.\end{array}$$

(6d)

In the case of the Rosenbrock function, the model formulation is defined in $\Omega =\left\{(x,y):x\in [-2,2],y\in [-1,3]\right\}$ and is given by:

$$f_9(x,y)=c_1^2(x,y)+c_2^2(x,y),$$

(7)

where $c_1$ and $c_2$ ideally are subfunctions of x and y in the form of:

$$\begin{array}{cc}{c}_1(x,y)=y-x^2,\hfill & \hfill c_2(x,y)=1-x.\end{array}$$

(8)

The ANN in the grey-box model is responsible for predicting the value of all subfunctions mentioned above and provides them into the model. The ANN is learned with the difference between the model output $\Xi (\mathit{x},\mathit{w})$ and the data generated from a benchmarking function $\widehat{\Xi }\left(\mathit{x}\right)$ in n sample points ${\mathit{x}}_i$, from which $m(m<n)$ are training data points, and $n-m$ are test data points. $\mathit{w}$ is the vector of weights and biases of a neural network. The learning process is performed by an evolutionary algorithm, here with the use of a CMA-ES library for Python as explained in Section 3.2. The vector $\mathit{w}$, which fully describes one network, is optimized by an optimization algorithm based on the fitness value:

$$F\left(\mathit{w}\right)=\sum _{i=1}^m{\left(\Xi \left({\mathit{x}}_i,\mathit{w}\right)-\widehat{\Xi }\left({\mathit{x}}_i\right)\right)}^2,$$

(9)

and the best network is chosen based on the whole dataset according to:

$$min\left({\left(\sum _{i=m+1}^n{\left(\Xi \left({\mathit{x}}_i,\mathit{w}\right)-\widehat{\Xi }\left({\mathit{x}}_i\right)\right)}^2\right)}^2+F^2\left(\mathit{w}\right)\right).$$

(10)

The dimensionality of the optimization of the fully connected network for the considered problem can be expressed with the following formula:

$$D=n_{\mathrm{in}}{n}_1+\sum _{i=2}^k{n}_{i-1}{n}_i+n_k{n}_{\mathrm{out}}+\sum _{i=1}^k{n}_i+2n_{\mathrm{out}},$$

(11)

where D is the optimization dimensionality, k is the number of hidden layers, $n_i$ is the number of neurons in the i-th hidden layer and $n_{\mathrm{in}}$ and $n_{\mathrm{out}}$ are the input and output dimensionality, respectively. For instance, the IMANN architecture for the polynomial prediction with one subfunction, the dimensionality of the optimization problem is equal to 47. The computational complexity of the CMA-ES algorithm is at least $\mathcal{O}\left(D^2\right)$. Every candidate solution consists of calling ANN for output, which is negligible, and a model function call. Assuming that a numerical model is a unit operation and runs at $\mathcal{O}\left(1\right)$, the overall computational complexity would be $\mathcal{O}\left(D^2\right)$. If computational complexity of a numerical model is dependent on some parameter ℓ and is equal to $\mathcal{O}\left(f_{\mathrm{NM}}\left(\ell \right)\right)$ with some unit operation, the overall computational complexity would be $\mathcal{O}\left(D^2{f}_{\mathrm{NM}}\left(\ell \right)\right)$.

5. Results

To quantify the difference between the benchmark function’s output and the predicted value, the integral of the squared error is used as an accuracy indicator:

$$R(\mathit{x},\mathit{w})={\int }_{\Omega }{\left(\Xi \left(\mathit{x},\mathit{w}\right)-\widehat{\Xi }\left(\mathit{x}\right)\right)}^2\mathrm{d}\mathit{x}.$$

(12)

The integral in Equation (12) is computed with eighty points per dimension Gauss-Legendre quadrature. The IMANN is compared to the DNN implemented with the use of the TensorFlow library [51]. The DNN learns in a standard way, treating the system as a black box. To neglect the stochastic error in the computations, one hundred attempts were performed for the IMANN and the DNN, for every numerousness of the dataset. The mean value of R from these 100 attempts is represented as $R_{\mathrm{m}}.$

$$D=n_{\mathrm{in}}{n}_1+\sum _{i=2}^k{n}_{i-1}{n}_i+n_k{n}_{\mathrm{out}}+\sum _{i=1}^k{n}_i+2n_{\mathrm{out}},$$

(13)

For both networks, the training to test dataset numerousness ratio was 7:3. The training set was chosen randomly for each attempt from the same dataset. The best network was chosen according to the algorithm described in the Section 3.3, and an R value was evaluated. An R value was not used for choosing the best network, because it cannot be estimated in real applications. R is used because it is a better measure of a network’s performance than the error made on the training and test data sets. The term ANN is used when referring to ANN part of IMANN and DNN for a typical ANN implementation.

5.1. Architecture

The implementation of the IMANN has a feedforward architecture consisting of fully connected layers with the last two layers called the part of the model layer (PM layer) and the model layer. In the input layer, the number of neurons corresponds to the conditions in which the system is located. Further layers are standard hidden layers. The number of layers and the number of neurons contained therein are selected according to the complexity of the problem modeled by the network. The next layer is the mentioned PM layer, where the number of neurons corresponds to the number of replaced parts of the model. The output from that layer multiplied by the weights is directly entered into the model layer. The numerical model also receives all arguments that were supplied to the network’s input, alongside the system parameters. The calculated numerical model result is the final output of the IMANN. Schematically, the IMANN’s architecture is presented in Figure 5. The implementation of the DNN was done by using the TensorFlow 1.12.0 library. The neural network is trained with the backpropagation algorithm using batches of size equal six. The optimizer used is Adam with default settings. For both DNN and IMANN, Bent Identity activation functions are used in all neurons. The IMANN’s architectures for problems (4a–d), (6a–d) and (8) were 1-5-5-1, 1-5-5-2 and 2-5-5-2, respectively. DNN’s architectures for problems (1) and (2) were 1-32-16-16-1 and 2-32-32-16-1, respectively. Architectures, as well as activation functions, were adjusted by trial and error for the best performance.

In terms of Figure 2 and Figure 5, x or ($x,y$) are the “System boundaries”, i.e., inputs, to both the numerical model and the ANN. “Measurable system parameters” can be viewed as any of the constants present in the corresponding problem. The “request for ANN predicted part of mathematical model” is the calling of ANN for a prediction, here without providing any additional inputs. “Part of the mathematical model” is the corresponding subfunction provided by the ANN. The model prediction is the value of $f_i$, and the “measurable system response” is a real value corresponding to the input.

5.2. Statistical Testing

The first paragraph covers the statistical comparison of IMANN and DNN algorithms over the case problems defined in Section 4. One must keep in mind that all IMANN implementations (i.e., IMANN’s designed for every model formulation $f_1$–$f_8$) were treated as if they were different algorithms. The algorithm and its prediction are represented with the same notation—subscript denotes problem (modeled function), and superscript indicates implementation (DNN or IMANN)—and the meaning is clear based on the context. Every algorithm has to minimize the errors of problems defined as a pair of functions and the dataset’s size. Errors were defined according to Equation (12). Means of errors obtained for the problems for every dataset are presented in

Table 1. Mean integral squared error value $R_{\mathrm{m}}$.

Dataset	${\mathit{f}}_1^{\mathbf{IMANN}}$	${\mathit{f}}_2^{\mathbf{IMANN}}$	${\mathit{f}}_3^{\mathbf{IMANN}}$	${\mathit{f}}_4^{\mathbf{IMANN}}$	${\mathit{f}}_5^{\mathbf{IMANN}}$	${\mathit{f}}_6^{\mathbf{IMANN}}$	${\mathit{f}}_7^{\mathbf{IMANN}}$	${\mathit{f}}_8^{\mathbf{IMANN}}$	${\mathit{f}}_{\mathbf{P}}^{\mathbf{DNN}}$
5	$6.189$	$2.317\times {10}^3$	$2.203\times {10}^4$	$4.362\times {10}^5$	$1.732\times {10}^4$	$3.991\times {10}^4$	$3.797\times {10}^4$	$5.850\times {10}^4$	$9.278\times {10}^4$
6	$3.658\times {10}^{-1}$	$1.126\times {10}^3$	$1.073\times {10}^4$	$2.699\times {10}^5$	$1.276\times {10}^3$	$1.521\times {10}^4$	$2.784\times {10}^4$	$3.810\times {10}^4$	$5.478\times {10}^4$
7	$2.476$	$5.610\times {10}^2$	$9.232\times {10}^3$	$2.566\times {10}^5$	$9.626\times {10}^2$	$1.868\times {10}^4$	$2.346\times {10}^4$	$3.196\times {10}^4$	$4.313\times {10}^4$
8	$1.767\times {10}^{-1}$	$4.275\times {10}^2$	$5.689\times {10}^3$	$1.658\times {10}^5$	$1.418\times {10}^2$	$5.983\times {10}^3$	$1.413\times {10}^4$	$2.517\times {10}^4$	$3.773\times {10}^4$
9	$5.255\times {10}^{-2}$	$1.860\times {10}^2$	$3.649\times {10}^3$	$1.007\times {10}^5$	$4.072\times {10}^1$	$3.841\times {10}^3$	$9.483\times {10}^3$	$1.839\times {10}^4$	$2.669\times {10}^4$
10	$5.108\times {10}^{-2}$	$4.905\times {10}^1$	$1.925\times {10}^3$	$7.772\times {10}^4$	$6.396$	$2.075\times {10}^3$	$5.054\times {10}^3$	$1.548\times {10}^4$	$1.935\times {10}^4$
11	$4.381\times {10}^{-2}$	$1.294\times {10}^2$	$2.556\times {10}^3$	$8.912\times {10}^4$	$1.699\times {10}^1$	$1.644\times {10}^3$	$5.097\times {10}^3$	$1.521\times {10}^4$	$1.926\times {10}^4$
12	$1.787\times {10}^{-2}$	$6.870\times {10}^1$	$1.538\times {10}^3$	$4.757\times {10}^4$	$3.637\times {10}^{-1}$	$8.952\times {10}^2$	$3.343\times {10}^3$	$1.197\times {10}^4$	$1.498\times {10}^4$
13	$1.698\times {10}^{-2}$	$2.118\times {10}^1$	$1.152\times {10}^3$	$3.079\times {10}^4$	$1.558\times {10}^{-1}$	$5.827\times {10}^2$	$2.729\times {10}^3$	$8.841\times {10}^3$	$1.339\times {10}^4$
14	$7.411\times {10}^{-3}$	$2.605\times {10}^1$	$1.516\times {10}^3$	$4.183\times {10}^4$	$1.089\times {10}^{-1}$	$2.632\times {10}^2$	$2.152\times {10}^3$	$1.018\times {10}^4$	$1.119\times {10}^4$
15	$9.115\times {10}^{-3}$	$1.723\times {10}^1$	$9.104\times {10}^2$	$3.369\times {10}^4$	$2.149\times {10}^{-1}$	$3.141\times {10}^2$	$2.487\times {10}^3$	$7.593\times {10}^3$	$9.632\times {10}^3$
16	$7.915\times {10}^{-3}$	$5.880$	$9.728\times {10}^2$	$2.366\times {10}^4$	$1.518\times {10}^{-1}$	$6.545\times {10}^1$	$1.177\times {10}^3$	$5.332\times {10}^3$	$8.875\times {10}^3$
17	$5.446\times {10}^{-3}$	$2.460\times {10}^1$	$7.609\times {10}^2$	$2.291\times {10}^4$	$3.772\times {10}^{-2}$	$3.529\times {10}^2$	$1.432\times {10}^3$	$6.060\times {10}^3$	$8.147\times {10}^3$
18	$6.344\times {10}^{-3}$	$2.899$	$7.958\times {10}^2$	$1.574\times {10}^4$	$8.170\times {10}^{-3}$	$2.649\times {10}^1$	$1.133\times {10}^3$	$4.631\times {10}^3$	$7.363\times {10}^3$
19	$5.733\times {10}^{-3}$	$4.288\times {10}^{-1}$	$5.663\times {10}^2$	$1.247\times {10}^4$	$1.355\times {10}^{-2}$	$1.795\times {10}^1$	$7.061\times {10}^2$	$4.825\times {10}^3$	$1.034\times {10}^4$
20	$4.135\times {10}^{-3}$	$2.295$	$5.869\times {10}^2$	$8.989\times {10}^3$	$1.359\times {10}^{-2}$	$9.447$	$3.349\times {10}^2$	$2.445\times {10}^3$	$6.013\times {10}^3$

and graphically in Figure 6 and Figure 7. Statistical testing was performed according to the guidelines indicated in Derrac et al. [52]. The Quade test was used for testing the null hypothesis, asserting the equality of medians between the distribution of errors made by different algorithms on overall problems to test if the algorithms differ significantly from each other [52]. Average Quade rankings are presented in

Table 2. Average rankings of the algorithms (Quade).

Algorithm	Ranking
$f_1^{\mathrm{IMANN}}$	9.000
$f_2^{\mathrm{IMANN}}$	7.331
$f_3^{\mathrm{IMANN}}$	5.581
$f_4^{\mathrm{IMANN}}$	1.000
$f_5^{\mathrm{IMANN}}$	7.669
$f_6^{\mathrm{IMANN}}$	5.294
$f_7^{\mathrm{IMANN}}$	4.125
$f_8^{\mathrm{IMANN}}$	3.000
$f_{\mathrm{P}}^{\mathrm{DNN}}$	2.000

. Quade statistic distributed according to F-distribution with 8 and 120 degrees of freedom is equal to 46.624. p-value computed by Quade Test is 2.6476 × 10⁻³³, which rejects the null hypothesis and, in the result, states that the algorithms’ performance differs significantly. The next step is to test the hypothesis of equality between all pairs of algorithms [52]. For this step, the Shaffer post-hoc procedure is used. z-values, associated p-values and adjusted p-values for $\alpha =0.05$ are presented in

Table 3. Shaffer Table for $\alpha =0.05$.

Algorithms	z	p	Adjusted p
$f_1^{\mathrm{IMANN}}$ vs. $f_4^{\mathrm{IMANN}}$	$8.26236$	$1.42815\times {10}^{-16}$	$5.14133\times {10}^{-15}$
$f_1^{\mathrm{IMANN}}$ vs. $f_{\mathrm{P}}^{\mathrm{DNN}}$	$7.22957$	$4.84530\times {10}^{-13}$	$1.35668\times {10}^{-11}$
$f_4^{\mathrm{IMANN}}$ vs. $f_5^{\mathrm{IMANN}}$	$7.03592$	$1.97951\times {10}^{-12}$	$5.54262\times {10}^{-11}$
$f_2^{\mathrm{IMANN}}$ vs. $f_4^{\mathrm{IMANN}}$	$6.39042$	$1.65428\times {10}^{-10}$	$4.63198\times {10}^{-9}$
$f_1^{\mathrm{IMANN}}$ vs. $f_8^{\mathrm{IMANN}}$	$6.19677$	$5.76324\times {10}^{-10}$	$1.61371\times {10}^{-8}$
$f_5^{\mathrm{IMANN}}$ vs. $f_{\mathrm{P}}^{\mathrm{DNN}}$	$6.00312$	$1.93556\times {10}^{-9}$	$5.41958\times {10}^{-8}$
$f_2^{\mathrm{IMANN}}$ vs. $f_{\mathrm{P}}^{\mathrm{DNN}}$	$5.35763$	$8.43221\times {10}^{-8}$	$2.36102\times {10}^{-6}$
$f_1^{\mathrm{IMANN}}$ vs. $f_7^{\mathrm{IMANN}}$	$5.03488$	$4.78152\times {10}^{-7}$	$1.33883\times {10}^{-5}$
$f_5^{\mathrm{IMANN}}$ vs. $f_8^{\mathrm{IMANN}}$	$4.97033$	$6.68395\times {10}^{-7}$	$1.87151\times {10}^{-5}$
$f_4^{\mathrm{IMANN}}$ vs. $f_6^{\mathrm{IMANN}}$	$4.71213$	$2.45141\times {10}^{-6}$	$5.39311\times {10}^{-5}$
$f_3^{\mathrm{IMANN}}$ vs. $f_4^{\mathrm{IMANN}}$	$4.45393$	$8.43123\times {10}^{-6}$	$1.85487\times {10}^{-4}$
$f_2^{\mathrm{IMANN}}$ vs. $f_8^{\mathrm{IMANN}}$	$4.32483$	$1.52649\times {10}^{-5}$	$3.35827\times {10}^{-4}$
$f_1^{\mathrm{IMANN}}$ vs. $f_3^{\mathrm{IMANN}}$	$3.80843$	$1.39850\times {10}^{-4}$	$3.07670\times {10}^{-3}$
$f_5^{\mathrm{IMANN}}$ vs. $f_7^{\mathrm{IMANN}}$	$3.80843$	$1.39850\times {10}^{-4}$	$3.07670\times {10}^{-3}$
$f_6^{\mathrm{IMANN}}$ vs. $f_{\mathrm{P}}^{\mathrm{DNN}}$	$3.67933$	$2.33844\times {10}^{-4}$	$5.14456\times {10}^{-3}$
$f_1^{\mathrm{IMANN}}$ vs. $f_6^{\mathrm{IMANN}}$	$3.55023$	$3.84888\times {10}^{-4}$	$8.08264\times {10}^{-3}$
$f_3^{\mathrm{IMANN}}$ vs. $f_{\mathrm{P}}^{\mathrm{DNN}}$	$3.42114$	$6.23603\times {10}^{-4}$	$1.12249\times {10}^{-2}$
$f_4^{\mathrm{IMANN}}$ vs. $f_7^{\mathrm{IMANN}}$	$3.22749$	$1.24883\times {10}^{-3}$	$2.24790\times {10}^{-2}$
$f_2^{\mathrm{IMANN}}$ vs. $f_7^{\mathrm{IMANN}}$	$3.16294$	$1.56186\times {10}^{-3}$	$2.81136\times {10}^{-2}$
$f_6^{\mathrm{IMANN}}$ vs. $f_8^{\mathrm{IMANN}}$	$2.64654$	$8.13202\times {10}^{-3}$	$1.30112\times {10}^{-1}$
$f_3^{\mathrm{IMANN}}$ vs. $f_5^{\mathrm{IMANN}}$	$2.58199$	$9.82327\times {10}^{-3}$	$1.57172\times {10}^{-1}$
$f_3^{\mathrm{IMANN}}$ vs. $f_8^{\mathrm{IMANN}}$	$2.38834$	$1.69247\times {10}^{-2}$	$2.53870\times {10}^{-1}$
$f_5^{\mathrm{IMANN}}$ vs. $f_6^{\mathrm{IMANN}}$	$2.32379$	$2.01368\times {10}^{-2}$	$2.61778\times {10}^{-1}$
$f_7^{\mathrm{IMANN}}$ vs. $f_{\mathrm{P}}^{\mathrm{DNN}}$	$2.19469$	$2.81858\times {10}^{-2}$	$3.66415\times {10}^{-1}$
$f_4^{\mathrm{IMANN}}$ vs. $f_8^{\mathrm{IMANN}}$	$2.06559$	$3.88671\times {10}^{-2}$	$4.66405\times {10}^{-1}$
$f_2^{\mathrm{IMANN}}$ vs. $f_3^{\mathrm{IMANN}}$	$1.93649$	$5.28075\times {10}^{-2}$	$5.80883\times {10}^{-1}$
$f_1^{\mathrm{IMANN}}$ vs. $f_2^{\mathrm{IMANN}}$	$1.87194$	$6.12146\times {10}^{-2}$	$6.12146\times {10}^{-1}$
$f_2^{\mathrm{IMANN}}$ vs. $f_6^{\mathrm{IMANN}}$	$1.67829$	$9.32900\times {10}^{-2}$	$8.39610\times {10}^{-1}$
$f_6^{\mathrm{IMANN}}$ vs. $f_7^{\mathrm{IMANN}}$	$1.48464$	$1.37638\times {10}^{-1}$	$1.00000$
$f_1^{\mathrm{IMANN}}$ vs. $f_5^{\mathrm{IMANN}}$	$1.22644$	$2.20031\times {10}^{-1}$	$1.00000$
$f_3^{\mathrm{IMANN}}$ vs. $f_7^{\mathrm{IMANN}}$	$1.22644$	$2.20031\times {10}^{-1}$	$1.00000$
$f_7^{\mathrm{IMANN}}$ vs. $f_8^{\mathrm{IMANN}}$	$1.16190$	$2.45278\times {10}^{-1}$	$1.00000$
$f_4^{\mathrm{IMANN}}$ vs. $f_{\mathrm{P}}^{\mathrm{DNN}}$	$1.03280$	$3.01700\times {10}^{-1}$	$1.00000$
$f_8^{\mathrm{IMANN}}$ vs. $f_{\mathrm{P}}^{\mathrm{DNN}}$	$1.03280$	$3.01700\times {10}^{-1}$	$1.00000$
$f_2^{\mathrm{IMANN}}$ vs. $f_5^{\mathrm{IMANN}}$	$0.64550$	$5.18605\times {10}^{-1}$	$1.00000$
$f_3^{\mathrm{IMANN}}$ vs. $f_6^{\mathrm{IMANN}}$	$0.25820$	$7.96253\times {10}^{-1}$	$1.00000$

. The z-value is used to determine the p-value from the normal distribution $\mathcal{N}(0,1)$. Adjusted p-value takes into account the accumulation of errors made by multiple comparisons. If adjusted p-value is less than the $\alpha$ value, the algorithms’ performance differs with the level of significance $\alpha$ equal to $0.05$. In Table 3, significant differences are separated with an additional line. It can be seen in the Table 3 that when the ANN in the IMANN has to predict non-complex function (model formulations $f_1$, $f_2$, $f_3$, $f_5$ and $f_6$), the grey box outperforms the DNN, which is indicated by the low adjusted p-value between the corresponding pairs. For more complex problems (model formulations $f_4$, $f_7$ and $f_8$) the difference is not significant. Another dependency worth noticing in Table 3 is that the number of parameters does not change the performance significantly (hypotheses of equality between pairs $f_1$ and $f_5$, $f_2$ and $f_6$, $f_3$ and $f_7$ and $f_4$ and $f_8$ cannot be rejected with 5% significance level).

Another test was performed to test the differences between standard DNN and IMANN formulations for Rosenbrock with two-dimensional input. Mean errors over one hundred trials obtained for every different dataset numerousness are presented in

Table 4. Mean integral squared error value $R_{\mathrm{m}}$.

Dataset	${\mathit{f}}_9^{\mathbf{IMANN}}$	${\mathit{f}}_{\mathbf{R}}^{\mathbf{DNN}}$
25	$1.327\times {10}^5$	$1.983\times {10}^5$
36	$4.969\times {10}^4$	$5.694\times {10}^4$
49	$2.033\times {10}^4$	$1.666\times {10}^4$
64	$2.403\times {10}^4$	$7.762\times {10}^3$
81	$8.267\times {10}^3$	$4.806\times {10}^3$
100	$3.299\times {10}^3$	$2.504\times {10}^3$
121	$4.432\times {10}^3$	$1.996\times {10}^3$
144	$2.440\times {10}^2$	$1.417\times {10}^3$
169	$3.205\times {10}^1$	$1.128\times {10}^3$
196	$1.426\times {10}^2$	$9.741\times {10}^2$
225	$1.312\times {10}^1$	$9.380\times {10}^2$
256	$1.502$	$6.607\times {10}^2$
289	$1.772\times {10}^{-1}$	$5.625\times {10}^2$
324	$1.059\times {10}^{-1}$	$5.409\times {10}^2$
361	$7.778\times {10}^{-1}$	$4.868\times {10}^2$
400	$3.285\times {10}^{-1}$	$4.457\times {10}^2$

and graphically in Figure 8. To compare these two formulations, Wilcoxon procedure was employed, as proposed by Derrac et al. [52]. Wilcoxon ranks are presented in

Table 5. Ranks of the algorithms (Wilcoxon).

Algorithm	Ranking
$f_9^{\mathrm{IMANN}}$	57
$f_{\mathrm{R}}^{\mathrm{DNN}}$	79

. The lower of the two should be less than or equal to the critical value obtained from the Wilcoxon distribution. The critical value obtained from Wilcoxon distribution for 16 degrees of freedom and the level of significance $\alpha =0.05$ is equal to 29, which is lower that 57; therefore, the null hypothesis of equality of means cannot be rejected. Still, a great increase of IMANN’s performance with enough data points presented in Figure 8, shows promising potential. High mean error value for IMANN, when the data numerousness is lower than 11 points per input dimension, is caused by 2 trials per 100 with very high error value. This problem may be addressed with a different learning approach.

5.3. Examples of Prediction

To illustrate the difference between the prediction of IMANN and DNN, an IMANN implementation for $f_3$ and $f_7$ is chosen. Mean with one standard deviation distance of the predicted parameter $\tilde{a}$ for $f_3$ is presented in Figure 9. One can see that the highest deviation is close to $x=0$, where the term $\tilde{a}{x}^3$ vanishes and thus does not count in error significantly. This observation is especially significant when one wants to interpret the value of the parameter predicted by the IMANN. In a domain where the parameter’s impact on system behavior is insignificant, its value should be treated as obtained with very low accuracy. Same figure type for two parameters $\tilde{a}$, $\tilde{b}$ for formulation $f_6$ is presented in Figure 10. Here the deviation between ideal parameters’ values and predicted by the IMANN is very high. The reason for such a difference is the ambiguity that arises when two parameters are predicted. In the Figure 11 the blue line and the blue area indicates the $f_{\mathrm{P}}$ mean prediction made by the IMANN and one standard deviation area, the green line and the green area indicates the $f_{\mathrm{P}}$ prediction made by the DNN and one standard deviation area while the red line is the $f_{\mathrm{P}}$. As can be seen, even though IMANN-predicted parameters have formed different than expected, the resulting prediction of the function $f_{\mathrm{P}}$ is highly accurate. This implies two important remarks. First, interpreting subfunctions’ values provided by IMANN should be checked for ambiguity, and second, proper constraints put on subfunctions’ values can reduce the number of ambiguous solutions. The constraints can be introduced by adding an additional penalty term in the error function during training. The penalty can correspond to subfunction values in specific points, range, or even shape described with derivative (e.g., monotonicity).

Comparison between IMANN and DNN mean predictions of Rosenbrock function for two different dataset numerousness are presented in Figure 12 for 36 data points (25 training points) and in Figure 13 for 256 data points (179 training points). For a low number of data points, the maximum discrepancy between prediction and original function for IMANN is of one order of magnitude lower than for DNN. For a higher number of points, the advantage of IMANN over DNN is greater, reaching the difference of two orders of magnitude.

6. Discussion

From the idea of the IMANN, one can see that by decreasing the load on the ANN part to zero in the IMANN will result in the IMANN’s prediction performance equal to the model (see Figure 2). In the proposed benchmarking case, the numerical model performance is perfect because we already know the form of the subfunctions. In real applications, finding even such a simple thing as the constant fitting parameters might be a problem.

Generalized computations from Section 5 prove that the IMANN can achieve higher performance than the DNN, as presented in Table 3. Gray-box approach is better than black-box when the complexity of subfunctions is lower than the complexity of a predicted function (see the upper part of Table 3). There was no significant difference when one and two subfunctions were modeled by the IMANN with similar complexity of subfunctions (see the lower part of Table 3. Therefore, the IMANN can improve the mathematical model’s performance by modeling its oversimplified, uncertain or missing parts, that are responsible for the discrepancy between the numerical model’s prediction and the modelled system.. The power of artificial neural network is utilized only for the problem part that cannot be solved accurately any other way. Such an approach reduces the complexity of a problem that is solved by an ANN, resulting in better prediction, especially with a small number of learning data points. The obtained results indicate great potential in the integration of mathematical models and artificial neural networks.

There are two major limitations of the IMANN. Since learning requires to run a numerical algorithm for each training data point each time a set of weights is tested, a reasonable learning time may not be obtained with computationally heavy numerical models. A second limitation arises from the fact that a problem of providing a part of the model by the ANN can become ambiguous, which can result in an unphysical or irrational form of subfunctions without additional control. A proposed way to keep ANN outputs physical is to impose additional constraints by adding penalty terms in $E_{\mathrm{train}}$. As for results provided in this paper, this topic is not studied.

7. Practical Application in Solid Oxide Fuel Cells Modeling

Microscale modeling of Solid Oxide Fuel Cell’s (SOFC) anode comes with a highly difficult problem of modeling a multistep electrochemical reaction between oxygen ions and hydrogen. The reaction occurs in the so-called triple phase boundary (TPB) and is a basis of SOFC’s operation. The equation which describes the reaction has empirical parameters that have to be estimated or fitted into the solution. Simplified equation (without elaborating the preexponential term $i_0$) takes the form described by the Butler-Volmer formula:

$$i_{\mathrm{TPB}}=i_0\left[exp\left(\frac{\alpha nF}{RT}{\eta }_{\mathrm{act}}\right)-exp\left(\frac{-\beta nF}{RT}{\eta }_{\mathrm{act}}\right)\right],$$

(14)

where $i_{\mathrm{TPB}}$ is the volumetric exchange current density ($\mathrm{A}/{\mathrm{m}}^3$), $i_0$ is the exchange current density ($\mathrm{A}/{\mathrm{m}}^3$), F = 96,485.3415 s A/mol is the Faraday constant, T is the temperature ($\mathrm{K}$), $\alpha$ and $\beta$ are the charge transfer coefficients (1) that are obtained by fitting to the experimental data [53], ${\eta }_{\mathrm{act}}$ is activation overpotential ($\mathrm{V}$). Detailed mathematical model of SOFC’s anode phenomena can be found elsewhere (e.g., [54]).

The Butler-Volmer equation comes with three empirical parameters—$i_0$, $\alpha$ and $\beta$, for which the values and the functional dependence from microstructural and system parameters is unknown [55]. Here, an ANN is employed to provide $\alpha$ and $\beta$ values as functions of temperature T (K) and current density I ($\mathrm{A}/{\mathrm{m}}^2$) for numerical model of the anode. The architecture is schematically presented in Figure 14.

Data available in the literature [56] was used as a source of microstructural parameters, typically used charge transfer coefficients’ values and both training and test datasets. Training set consisted of twelve data points, six for $800 °\mathrm{C}$ and six for $1000 °\mathrm{C}$. The test set was made of six data points for temperature $900 °\mathrm{C}$.

The ANN’s architecture contains one hidden layer consisting of three neurons with Bent Identity activation functions. The output layer contains two neurons—one for $\alpha$ and one for $\beta$ charge transfer coefficients—with SoftPlus activation functions, as values of charge transfer coefficients should be positive. The output of activation functions was directly provided to the numerical model. The inputs to the ANN were temperature and current density.

To check if the IMANN architecture provides values that are physically feasible, it is shown what is the proposed functional form of charge transfer coefficients in dependence from current density, Figure 15, and temperature, Figure 16. In Figure 15 values of both $\alpha$ and $\beta$ charge transfer coefficients are presented as functions of current density for three different operating temperatures. It can be seen that the $\alpha$ value proposed by the ANN changes nonlinearly with the current density and is increasing monotonically. In Figure 16 charge transfer coefficients are presented as functions of temperature for three different current densities. The $\alpha$ value changes nonlinearly with the temperature and is decreasing monotonically. The $\beta$ value is practically independent from both temperature and current density. Values of charge transfer coefficients proposed by the ANN are close to the ones frequently used in the literature. Researchers in the literature propose different values for charge transfer coefficents. For $\alpha$ and $\beta$, one can find values in the ranges $[0.7,1]$ and $[0.1,0.5]$, corespondingly [55,57,58].

To see the improvement of the overall model, IMANN architecture is compared to the model with typically used charge transfer coefficient values. In the Figure 17 overpotential of the anode versus current density is presented for three different temperatures taken from the literature [56]. Experimental data points are marked as symbols, the standard model is represented as a dashed line, and the proposed method is represented as a solid line. It can be seen that for all data points—both training and test—hybrid model outperforms the classical model in terms of agreement with experimental measurements.

8. Conclusions

The present paper introduced a framework for integrating a classical mathematical model and an artificial neural network, using an evolutionary algorithm to limit the required dataset. The methodology can be applied to any system only if a part of it can be expressed in the form of mathematical equations. The combination of an artificial neural network and a classical mathematical model is interactive, expressed in the possibility to exchange information between a numerical model and ANN as a numerical model’s request and ANN’s response. The Interactive Mathematical Model-Artificial Neural Network was employed to predict several benchmark functions’ values when given a different number of training data. Problems presented in this work were narrowed to theoretical cases with limited complexity compared to real applications. However, the proposed approach can be applied to simulate any system for which the mathematical model is incomplete and a limited dataset of experimental data is available. This situation is widespread in engineering. A practical application in the Solid Oxide Fuel Cell modeling was successfully presented, in which parameters provided by ANN were physically feasible and model prediction was improved. Current limitations include computational complexity due to the use of the derivative-free method for problems of high dimensionality, which requires to call a potentially numerical-heavy numerical model thousands of times. This problem can be addressed by searching for the most optimal optimization method in this application. Another limitation that was described corresponds to the specific problem, not to the IMANN idea. Therefore, it has to be controlled during application by forcing physical and realistic values and comparing the proposed by the IMANN solution to the available experimental or literature data. In future work, the problem of dimensionality should be addressed, i.e., making the possibility of using ANN architectures with a greater number of neurons and layers. More complicated IMANN architectures, including multiple neural networks assigned to different model parts, convolutional neural networks to reduce the number of optimized coefficients or feedback loops (like recurrent ANNs) to achieve additional information for ANN, like the model’s state, can be investigated.