Comparison of the Probabilistic Ant Colony Optimization Algorithm and Some Iteration Method in Application for Solving the Inverse Problem on Model With the Caputo Type Fractional Derivative.

This paper presents the algorithms for solving the inverse problems on models with the fractional derivative. The presented algorithm is based on the Real Ant Colony Optimization algorithm. In this paper, the examples of the algorithm application for the inverse heat conduction problem on the model with the fractional derivative of the Caputo type is also presented. Based on those examples, the authors are comparing the proposed algorithm with the iteration method presented in the paper: Zhang, Z. An undetermined coefficient problem for a fractional diffusion equation. Inverse Probl. 2016, 32.

1. Introduction

The inverse problems are an important and commonly encountered class of problems in many branches of technology and mathematics [1,2]. Solving them involves in the appropriate choice of the input parameters of the model (e.g., material coefficients or initial-boundary conditions) in order to model the process in such a way that the suitable output is obtained.

In this paper, we consider the model with the fractional differential equation with the fractional derivative [3,4] of the Caputo type. This type of derivative is being used to describe the phenomenon of the anomalous diffusion. The example of this kind of phenomena is a heat conduction in the porous material. Voller [5] shows that for the materials of anomalous porosity, the models containing the fractional derivative are better for describing the problem that the models with the classical derivative. The experiments conducted in the papers [6,7] show that the fractional derivatives make sense when describing the phenomenon of the heat conduction in such materials as the composites or the porous aluminum. In [6], the authors consider the inverse heat conduction problem in the composite material. For the model presented in this paper the fractional derivative of the Caputo type was used. In [7], the authors are comparing the mathematical models containing various derivatives (classical, fractional of the Riemann–Liouville type and fractional of the Caputo type) basing on the experimental data obtained for the porous aluminum. Moreover, in this case the fractional derivative turned out to be better for describing the phenomenon. More on the application of the fractional calculus in various fields of science can be found in the papers [8,9,10,11,12,13,14].

In this paper, we present the algorithm for solving the inverse problem on the model containing the fractional derivative of the Caputo type. In the presented model, the thermal conductivity coefficient is being reconstructed. For solving the direct problem the finite difference method with the approximation of the Caputo type derivative [15] was used. The algorithm for solving the inverse problem consists in minimizing the functional representing the approximation error. For this purpose, one of the algorithms of the swarm intelligence was used [16,17]. In this paper, the numerical examples used for comparison of the presented algorithm with the algorithm shown in [18] are also being presented.

2. Formulation of the Problem

We consider two methods for solving the inverse problem on the models containing the fractional derivatives. Let us consider the following differential equation with the fractional derivative and the initial-boundary conditions.

In the above Equation there is a fractional derivative with respect to time. In this paper, this derivative is of order and of a Caputo type which is defined as follows,

Equation (1) (with the fractional derivative of order with respect to time) with the initial-boundary conditions (2)–(4) may be used to model the heat conduction phenomenon in the porous material. Then, the respective components of the model take the following names.

—spatial and time variable,

u—function representing the temperature distribution,

c—specific heat,

—density,

—thermal conductivity,

—order of fractional derivative,

g—additional heat source,

f—function representing the initial condition,

—functions representing the boundary condition of the first type.

Generally speaking, the solution of the inverse heat conduction problem consists in such a choice of the model parameters (model input) so the given output temperature distribution is obtained. More information regarding the heat processes modeling in the porous materials, inverse heat conduction problems, and use of the fractional calculus for modeling these kind of processes may be found in [5,7,19,20,21].

Description of the Considered Inverse Problem

In this paper, we consider two algorithms for solving the inverse problem. The description of the first method cas be found in the paper [18]. This method is iterative so in this paper we will call it the iteration method. The second algorithm is the algorithm proposed by the authors of this paper and it is based on the Real Ant Colony Optimization algorithm [16], which is the swarm artificial intelligence algorithm. The details of the algorithm are being described in the Section 3.2. More about application of the swarm intelligence algorithms can be found in [17,22,23].

In the example the thermal conductivity dependent on time coefficient which occurs in the Equation (1) is being reconstructed. The additional information on the basis of which the coefficient is being searched is the condition where the function a is given and the function b (the values of the function b at the points of the grid) is being calculated from the formula by first solving the direct problem for the exact value of the parameter . The values are treated as the input data for the inverse problem. For the determined thermal conductivity coefficient by solving the direct problem we obtain the values , and therefore we also obtain the values . By comparing the values and we construct the algorithm:

By minimizing this functional (using the algorithm proposed by the authors) we are reconstructing the thermal conductivity coefficient .

3. Methods of Solution

In the process of solving the described inverse problem there is a need to repeatedly solve the direct problem. In this section, we first describe the solution of the direct problem and then describe the Real Ant Colony Optimization algorithm needed for minimizing the functional (7).

3.1. Solution of the Direct Problem

The direct problem described by the Equation (1) and the initial-boundary conditions (2)–(4) were resolved by the use of the finite difference method.

Let be the size of the grid respectively in the domain of space and time. We take the following grid steps; . Then, the points of the grid from the space interval are the numbers and from the time interval the numbers . The values of the functions at the grid points are denoted as follows: , and . The values of the approximation function at the points are denoted as . The Equation (1) may be written as follows, where in terms of the heat conductivity is the thermal diffusivity coefficient and . The fractional derivative of the Caputo type of the order with respect to time is approximated by the following Equation [15], where

Using the Equation (9), the boundary conditions and the differential quotient for the second derivative , after proper organizing the elements we obtain the following discrete differential equations ():

The Equation (10) may be also written in the matrix form where , A is a tridiagonal matrix of the size and following coefficients, and D is a vector , where:

3.2. Solution of the Inverse Problem

In order to minimize the functional (7), the Real Ant Colony Optimization algorithm [16] is being used. The inspiration for the algorithm was the swarming behavior of ants in nature. In the presented algorithm, the role of the solutions are being played by the pheromone spots in the number of L. At the beginning those spots are randomly placed in the considered region. Next, they are being sorted by the quality and certain probabilities are being assigned to them. The better the solution is, the higher the probability of choosing it is. This way we create the set of solutions. In every iteration each of M ants construct one new solution (new pheromone spot) using the probability density function (Gauss function). First the ant chooses with the certain probability the solution that will be transformed. This probability is dependent on the parameter q of the algorithm. If value of q is small then the choice of the best solution is preferred and the higher the value of q is, the more similar to each other the probabilities of choosing any solution are. When the solution is chosen then the ant samples its neighborhood using the Gauss function with the parameters and , where is dependent on the parameter . This parameter is responsible for the disappearance of the pheromone trace. The smaller value of is, the faster the algorithm converges to its best solution at the expense of further exploration of the region. Then, the set of solutions is updated with the new solutions—the new and previous solutions are being sorted by quality and the M of the worst solutions get rejected. Moreover, the algorithm was adapted for parallel computing. Mathematical fundaments and results concerning the convergence of the different versions of the ACO algorithm are included in [24,25,26,27]. Based on the analysis of the papers of various authors, see for example [28,29,30,31,32], as well as on the computation carried out by us (for example [7,17,20,22,31]) we come to the conclusion that using the artificial intelligence algorithms (ACO, RACO, RACO, IRM, etc.) in most cases does not require the application of the explicit regulatory procedures. There is no mathematical proof for that, only the empirical data. For the examples presented in the paper we have known the exact solutions so we could compare the approximate solutions with them. The reconstruction errors of the reconstructed parameters generally do not exceed the input data errors. They are slightly larger in some cases. It follows that the results obtained for the disturbed input data remain stable. In order to describe the algorithm we are introducing the notation

We will present the consecutive steps of the algorithm.

Initialization of the algorithm

Setting the input parameters of the algorithm: .

Generating L solutions playing the role of the pheromone spots. Assigning them to the set (initial archive).

Computing the values of the objective function for all of the pheromone spots (parallel computing) and sorting the archive from the best solution to the worst.

Iteration process

Assigning the probabilities to the pheromone spots according to the pattern: where weights are related to the solution number l and are expressed by the formula:

The ant randomly chooses the th solution with the probability .

The ant transforms the j-th coordinate () of the th solution by sampling the neighborhood using the probability density function (in this case the Gauss function): where .

Steps 5–6 are being repeated by every ant. We get M new solutions (pheromone spots) by that.

Partition of the new solutions into the groups. Computing the value of the objective function for the new solutions (parallel computing).

Adding the new solutions to the archive , sorting them by the quality and rejecting the M worst solutions.

The steps 3–9 are being repeated I times.

Knowing the values of the parameters and I we can determine the number of the objective function calls during the algorithm execution. This number is . The Figure 1 presents the block diagram of the RealACO algorithm.

Figure 1

The block diagram of the Real Ant Colony Algorithm.

4. Numerical Examples

In this section, we present two numerical examples that we used for comparing the algorithm proposed by the authors (based on the Real Ant Colony Optimization algorithm) with the iteration algorithm presented in the paper [18]. In both examples we are considering the model (1)–(4), and the reconstructed coefficient is thermal conductivity dependent on time. The input data for the inverse problem is a set of values (6).

4.1. Example 1

For this example we take the following data:

The reconstructed coefficient is in the following form:

The exact values of the parameters , , are 1, , 1. For minimizing the functional (7) the RealACO was applied with the parameters , , . The values of the parameters were chosen basing on the experiments that we had carried out earlier. The number of the objective function calls is 606. Table 1 presents the results of the reconstruction.

Table 1

Results of calculations for the RealACO algorithm (grid ) (—reconstructed value of the coefficient , —the relative error of reconstruction of the coefficient , —standard deviation ()).

Noise	a¯_i	δ_a¯_i[%]	σ
0%	0.9999	0.01	4.36×10^−4
	0.4999	0.02	7.81×10^−5
	0.9999	0.01	4.36×10^−5
1%	1.0001	0.01	1.77×10^−4
	0.4998	0.04	1.74×10^−4
	0.9999	0.01	7.21×10^−5
2%	0.9995	0.05	1.96×10^−4
	0.4996	0.08	3.67×10^−4
	0.9997	0.03	1.96×10^−5
5%	1.0001	0.01	3.79×10^−4
	0.4996	0.08	5.82×10^−5
	1.0003	0.03	4.84×10^−5

The error of the obtained results was calculated according to formula where N indicates the number of measurements, is the reconstructed coefficient and is the exact value of the coefficient. The grid used to generate the input data was of the size as well as the grid used in the algorithm for solving the direct problem. Table 2 presents the results in dependence on the size of input data disturbance.

Table 2

Comparison of the algorithms according to the obtained errors of reconstructing the coefficient for the input data from the example 1.

Noise	RealACO	Iteration Method
	Error	ε	Error
0%	4.2781×10^−4	0.01	4.8035×10^−3
		0.001	4.0918×10^−4
		0.0001	3.5545×10^−5
1%	3.3106×10^−3	0.01	9.0230×10^−2
		0.001	8.9619×10^−2
		0.0001	8.9578×10^−2
2%	1.0889×10^−2	0.01	1.8033×10^−1
		0.001	1.8057×10^−1
		0.0001	1.8059×10^−1
5%	8.7965×10^−3	0.01	4.4052×10^−1
		0.001	4.4016×10^−1
		0.0001	4.4012×10^−1

In case of the exact input data the smaller value of the parameter is, the smaller error of the result obtained by the iteration algorithm we get. For , the approximate error is , while for , the error is . In case of the RealACO algorithm the obtained error is approximately and it is slightly bigger than the error obtained by the iteration algorithm. For the input data disturbed by the pseudorandom error decreasing the parameter in the iteration method slightly reduces the errors of the obtained results. In case of the disturbed input data the reconstruction errors of the thermal conductivity coefficient when using the algorithm proposed in the paper were smaller than when using the iteration method. For the data disturbed by the pseudo-random error the RealACO algorithm reconstructed the coefficient with the error while for the iteration method the error was .

4.2. Example 2

The input data used in the previous example were taken from the paper [18]. On the contrary, this example was not presented in that paper. In this example we compare both algorithms with the input data generated on a grid of a different density than the grid used in the algorithms reconstructing the coefficient. In such a case, testing the algorithm is a base for avoiding the inverse crimes [33] (page 5). The grid used for generating the input data was of size and the grid used in the reconstruction process was of size . The algorithm described in the paper [18] was implemented by us in the C# and then the computation was carried out. The results are presented in the Table 3.

Table 3

Comparison of the algorithms according to the obtained errors of reconstructing the coefficient for the input data from example 2.

Noise	RealACO	Iteration Method
	Error	ε	Error
0%	0.6230	0.01	1.1821
		0.001	1.1825
		0.0001	1.1826
1%	0.5982	0.01	1.2050
		0.001	1.2054
		0.0001	1.2055
2%	0.6146	0.01	1.2369
		0.001	1.2374
		0.0001	1.2375
5%	0.6566	0.01	1.4903
		0.001	1.4907
		0.0001	1.4908

The computational data are as follows.

The coefficient is being reconstructed in the following form, where the exact values of the parameters , , and are 1, 4, and , respectively. In the RealACO algorithm, the parameters , , and were used. The values of the parameters were chosen basing on the experiments that we had carried out earlier. Table 4 presents the results obtained for the RealACO algorithm.

Table 4

Results obtained by the RealACO algorithm (grid ) ( – reconstructed value of a coefficient ,  – relative error of reconstruction of the coefficient ,  – standard deviation ()).

Noise	a¯_i	δ_a¯_i[%]	σ
0%	0.9889	1.10	0.1077
	3.9999	0.01	0.0124
	1.0875	1.14	0.0507
1%	0.9910	0.90	0.0128
	3.9991	0.03	0.0040
	1.0874	1.15	0.0140
2%	0.9900	0.99	0.0266
	3.9998	0.01	0.0076
	1.0873	1.15	0.0156
5%	0.9881	1.18	0.0277
	4.0002	0.01	0.0182
	1.0869	1.18	0.0150

The result of the reconstruction obtained from the RealACO algorithm are very well, the most reconstruction errors are smaller than the input data errors. Even for the input data disturbed by the error the reconstruction errors are no more that . Moreover, the values of the standard deviation for the obtained results are low.

Table 3 presents the results used for comparison of the algorithms for various disturbances of the input data. In case of the iteration method the change of the parameter insignificantly affects the reconstruction errors of the coefficient . In each considered case the reconstruction errors of the coefficient are smaller for the RealACO algorithm than for the iteration method.

Finally, we present how the solutions in the respective iterations converge to the approximate solution for the Real Ant Colony Optimization algorithm (Figure 2 and Figure 3). From about the twentieth iteration the value of the reconstructed coefficient as well as the value of the functional start changing insignificantly.

Figure 2

Values of the coefficients (a) and (b) for the input data disturbed by the error.

Figure 3

Value of the coefficient (a) and value of the functional J (b) for the input data disturbed by the error.

5. Conclusions

The algorithm for solving the inverse problems on the models described by the differential equation with the Caputo type derivative presented in the paper seems to be an effective tool for solving these type of problems. The presented method comes to minimizing the functional representing the approximation error. For the functional, the classical methods of finding the extrema cannot be used and because of that the Real Ant Colony Optimization algorithm was used in the paper. The reconstruction errors obtained in the presented examples were small and did not exceed the input data disturbances. Additionally, the presented method turned out to be more precise than the method presented in the article [18], especially assuming that the input data was generated on different grid than the grid used for solving the direct problem.