Nowadays, when all kinds of computational and simulation methods are becoming more and more important and the processing power of computers is increasing, it is worth developing and looking for new applications of this type of tools. In practice, there are many classical numerical methods that are successfully used, developed and adapted to current problems. On the other hand, in recent years, scientists have been developing artificial intelligence methods, which include metaheuristic optimization algorithms. In this paper, we focused on combining a classical numerical method for solving a differential equation and selected heuristic algorithms for solving the inverse problem.

In the scientific literature, many works on the use of fractional derivatives to model many processes occurring in physics and engineering can be found (Bhangale et al., 2023; Brociek et al., 2017, 2019; Sowa et al., 2023; Bohaienko and Gladky, 2023; Koleva et al., 2021). In particular, derivatives of this type are used to model anomalous diffusion. As an example, heat flow in porous materials can be mentioned (Brociek et al., 2017, 2019). In article (Brociek et al., 2019), the authors model the phenomenon of heat flow in a porous medium. For this purpose, several mathematical models were used, including those based on fractional derivatives. The results from numerical experiments were compared with measurement data. Models that used fractional derivatives proved to be more accurate than the model with integer-order derivatives. Sowa et al. (2023) used fractional calculus and presented a voltage regulator extended model. The approach presented in the paper for modelling a voltage regulator has been verified with measurement data, and the non-integer order Caputo derivative proved to be an effective tool. Another application of fractional derivatives in mathematical modelling can be found in the paper (Bohaienko and Gladky, 2023), where a model for predicting the dynamics of moisture transport in fractal-structured soils was presented. The model incorporates the Caputo derivative. The numerical solution was obtained using the Crank-Nicholson finite-difference scheme. More examples and trends in the application of fractional calculus in various scientific fields are available in the publications (Hristov, 2023b; Obembe et al., 2017; Ionescu et al., 2017). Publications related to numerical methods in the field of fractional calculus are also worth mentioning (Ciesielski and Grodzki, 2024; Hou et al., 2023; Arul et al., 2022; Hristov, 2023a; Podlubny, 1999; Stanisławski, 2022).

In many engineering problems, there is a need to solve what’s known as an ‘inverse problem’. In simple terms, these issues involve identifying the input parameters of a model (e.g. boundary conditions or material parameters) based on observations (outputs) of the model. Typically, these problems are challenging because they are ill-posed (Aster et al., 2013c; Kaipio and Somersalo, 2005). Examples of inverse problems in various applications are included in the articles (Montazeri et al., 2022; Brociek et al., 2024; Ashurov et al., 2023; Wang et al., 2023; Ibraheem and Hussein, 2023; Brociek et al., 2023; Magacho et al., 2023). For example, the article (Brociek et al., 2024) focuses on an inverse problem concerning the identification of the aerothermal heating of a reusable launch vehicle based on temperature measurements taken in the thermal protection system (TPS) of this vehicle. The mathematical model of the TPS presented in the paper takes into account the dependence on temperature of the material parameters, as well as the thermal resistances occurring in the contact zones of the layers. To solve the inverse problem, the Levenberg-Marquardt method was applied.

One approach to solving inverse problems is to create a fitness function (or loss function, error function), and then optimize it to find the identified parameters. In this context, metaheuristic optimization algorithms can be very effective. In the paper (Hassan and Tallman, 2021), the authors utilize genetic algorithms, simulated annealing, and particle swarm optimization to solve the piezoresistive inverse problem in self-sensing materials. The considered problem was ill-posed and multi-modal. The results obtained in the study indicate that the genetic algorithm proved to be the most effective. As another example, the article (Al Thobiani et al., 2022) addresses an inverse problem for crack identification in two-dimensional structures. The authors utilize the eXtended Finite Element Method (XFEM) associated with the original Grey Wolf Optimization (GWO) and an improved GWO using Particle Swarm Optimization (PSO) (IGWO). The utilization of heuristic optimization algorithms for inverse problems in models with fractional derivatives can be found in papers (Brociek et al., 2020; Brociek and Słota, 2015). In both of these articles, the Ant Colony Optimization (ACO) algorithm was applied. The first publication involved a comparison with an iterative method. In the second article, heat flux on the boundary was identified. Additionally, papers (Kalita et al., 2023; Alyami et al., 2024; Zhang and Chi, 2023) address metaheuristic optimization algorithms and their applications.

In this article, an algorithm for solving the inverse problem for the equation of anomalous diffusion is presented. This equation is a partial differential equation with fractional derivatives. The Caputo derivative was adopted as the derivative with respect to time, while the Riemann-Liouville derivative was utilized for the derivative with respect to space. In the considered inverse problem, the objective was to identify the function appearing in the fractional boundary condition as well as the orders of the derivatives. To achieve this, several metaheuristic optimization algorithms were used and compared.

In this article, we consider a mathematical model of anomalous diffusion with fractional derivatives both in time and space. Models of this type can effectively be used to model mass and heat transport phenomena in porous media (Brociek et al., 2019; Sobhani et al., 2023; Kukla et al., 2022). The model consists of the following fractional-order partial differential equation: (1) C ∂ α u ( x , t ) ∂ t α = λ ( x , t ) R L ∂ β u ( x , t ) ∂ x β + f ( x , t ) , x ∈ ( 0 , L ) , t ∈ ( 0 , T ] . The equation is supplemented with the initial condition: (2) u ( x , 0 ) = φ ( x ) , x ∈ [ 0 , L ] , and boundary conditions: (3) u ( 0 , t ) = 0 , t ∈ ( 0 , T ] , (4) u ( L , t ) + ( λ ( x , t ) R L ∂ β − 1 u ( x , t ) ∂ x β − 1 ) x = L = ψ ( t ) . It is important to note the boundary condition (4) at the right boundary of the considered domain. It takes the form of a Robin boundary condition with a fractional derivative. In the differential equation (1), two different fractional-order derivatives were assumed. The derivative of order α ∈ ( 0 , 1 ) with respect to time is a Caputo-type derivative, defined by the following formula: (5) C ∂ α u ( x , t ) ∂ t α = 1 Γ ( 1 − α ) ∫ 0 t ∂ f ( x , s ) ∂ s ( t − s ) − α d s . In the case of the derivative with respect to space, a fractional-order derivative of order β ∈ ( 1 , 2 ) Riemann-Liouville type was applied: (6) R L ∂ β u ( x , t ) ∂ x β = 1 Γ ( 2 − β ) ∂ 2 ∂ s 2 ∫ 0 x f ( s , t ) ( t − s ) 1 − α d s . In the model, it is also assumed that λ is a continuous, positive function called the diffusion coefficient, and the functions f, φ, and ψ are also continuous functions. To effectively model and conduct simulations in such models, the first step is to solve equations (1)–(4). This task is known as the direct problem (or forward problem). In the next section, a numerical method to solve the considered equation is described.

This article primarily focuses on the inverse problem. The presented approach requires optimization of the fitness function. However, in the optimization process, it is necessary to repeatedly solve equations (1)–(4), that is, the so-called forward problem. To solve it, an implicit finite difference scheme is applied.

Firstly, the considered domain Ω = [ 0 , T ] × [ 0 , L ] is discretized, resulting in a grid S = { ( x i , t k ) : i = 0 , 1 , … , N , k = 0 , 1 , … , K } with steps Δ x = L N , x i = i Δ x, Δ t = T K , t k = k Δ t. Then, for all functions involved in the model, the values at the grid points S are determined. We use the following notation: λ i k = λ ( x i , t k ), f i k = f ( x i , t k ), φ i = φ ( x i ), ψ k = ψ ( t k ). Let u i k = u ( x i , t k ) denote the values of the exact solution at the points ( x i , t k ), and U i k represent the corresponding values obtained from the numerical solution.

To derive the implicit finite difference scheme, we need to apply the approximation of the Riemann-Liouville derivative (6) in the form of the shifted Grünwald formula (Tian et al., 2015; Tadjeran and Meerschaert, 2007): (7) R L ∂ β u ( x , t ) ∂ x β ( x i , t k + 1 ) ≈ 1 ( Δ x ) β ∑ j = 0 i + 1 g β , j U i − j + 1 k + 1 , where g β , j = Γ ( j − β ) Γ ( − β ) Γ ( j + 1 ) . Then, we approximate the Caputo derivative (5) using the following formula (Lin and Xu, 2007): (8) C ∂ α u ( x , t ) ∂ t α ( x i , t k + 1 ) ≈ 1 Γ ( 2 − α ) ( Δ t ) α ∑ j = 1 k ( j 1 − α − ( j − 1 ) 1 − α ) ( U i k − j + 1 − U i k − j ) . The derivative appearing in the boundary condition (4), after considering the zero condition at the left boundary, is approximated as follows: (9) R L ∂ β − 1 u ( x , t ) ∂ x β − 1 ( x N , t k + 1 ) ≈ 1 ( Δ x ) β − 1 ( ∑ j = 1 N g β − 1 , j U N − j + 1 k + 1 + g β − 1 , 0 ( 3 U N k + 1 − 3 U N − 1 k + 1 + U N − 2 k + 1 ) ) .

Combining all of the above approximations together and after appropriate transformations, we obtain the implicit difference scheme, which can be expressed in matrix form as: (10) A 1 U 1 = Q 0 + F 1 , (11) A k + 1 U k + 1 = ( 1 − b 1 ) Q k + ∑ j = 1 k − 1 ( b j − b j + 1 ) Q k − j + b k Q 0 + F k + 1 , where b j = ( j + 1 ) 1 − α − j 1 − α . And the vectors U k , Q k , F k have the following form: U k = [ U 1 k , U 2 k , … , U N k ] T , Q k = [ U 1 k , U 2 k , … , U N − 1 k , 0 ] T , F k = [ ( Δ t ) α Γ ( 2 − α ) f 1 k , ( Δ t ) α Γ ( 2 − α ) f 2 k , … , F k = [ ( Δ t ) α Γ ( 2 − α ) f N − 1 k , ( Δ x ) β − 1 ψ k ] . The matrix A k , dependent on the time step k, is of size N × N. Its coefficients are determined by the following formula: (12) a i j k = − r λ i k + 1 g β , i − j + 1 , 1 ⩽ j ⩽ i − 1 , 1 ⩽ i ⩽ N − 1 , 1 − r λ i k + 1 g β , 1 , 1 ⩽ j = i ⩽ N − 1 , − r λ i k + 1 g β , 0 , j = i + 1 , 1 ⩽ i ⩽ N − 1 , λ N k + 1 g β − 1 , N − j + 1 , 1 ⩽ j ⩽ N − 3 , i = N , λ N k + 1 g β − 1 , 3 + λ N k + 1 g β − 1 , 0 , j = N − 2 , i = N , λ N k + 1 g β − 1 , 2 − 3 λ N k + 1 g β − 1 , 0 , j = N − 1 , i = N , ( Δ x ) β − 1 + λ N k + 1 g β − 1 , 1 + 3 λ N k + 1 g β − 1 , 0 , j = i = N , 0 , i + 2 ⩽ j ⩽ N , 1 ⩽ i ⩽ N − 2 . In the above formula, the symbol r is defined as r = ( Δ t ) α Γ ( 2 − α ) ( Δ x ) β . Equations (10)–(11) define systems of linear equations. Solving these systems will yield the approximate values of the function u at the grid points S. Article (Xie and Fang, 2020) includes theorems regarding the stability and convergence of the numerical scheme (10)–(12). In the case of the scheme under consideration, the convergence order is O ( ( Δ t ) 2 − α + ( Δ x ) 2 ).

In mathematical modelling, as well as in various computer simulations it is essential to use proper mathematical models. In this paper, time-space fractional diffusion equation (TSFDE) with fractional boundary contition is considered. This model was further described in Section 2 and can be used as an effective tool in modelling the heat conduction in porous materials (Brociek et al., 2019). To solve model described by (1)–(4) it is necessary to know full information about the model input, such as material’s parameters, geometry or model coefficients. In many engineering issues, it is impossible to have the knowledge of all model’s information. It might be because of the lack of measuring equipment, or toughness in choosing the parameters. The usual problem is the choice of such entry model’s parameters – input, that the model’s result – output (e.g. temperature measurements in a chosen control point) takes the proper value. Problems of this sort are named inverse problems and are usually hard to solve (Özişik and Orlande, 2021; Aster et al., 2013c, 2013a, 2013b). To put it simply, the problem is the identifying of the input parameters to fit to the measurement data (part of model’s output) as closely as possible.

In this article, in order to solve the inverse problem, an approach is presented which involves optimizing the following fitness function: (13) F ( a 0 , a 1 , … , a dim ) = ∑ w = 1 W ( u w ( a 0 , a 1 , … , a dim ) − u ‾ w ) 2 . Symbols a 0 , a 1 , … , a dim are marked as unknown input parameters of the model-parameters, which are to be identified. Objective function F is dependent on these parameters. By dim + 1 the size of optimization task is marked, W is a number of data in model’s output (e.g. measurement data), u w ( a 0 , a 1 , … , a dim ) ( w = 1 , 2 , … , W ) denotes output values calculated from the model for fixed input parameters a 0 , a 1 , … , a dim and by u ‾ w the data output (measurement data) is marked, to which model should fit itself. So, function F measures “how close” are the calculated values from a model (for fixed input parameters a 0 , a 1 , … , a dim ) to output values given in a problem (e.g. measurement data). Solving given inverse problem is based on finding the minimum of fitness function relative to unknown parameters a 0 , a 1 , … , a dim . Hence, the use of selected metaheuristic algorithms for finding the minimum of the fitness function is justified. In Section 6, a computational example is presented and algorithms’ efficiency is discussed. Figure 1 schematically presents the data flow in forward and inverse problems.

Heuristic optimization algorithms for searching objective function’s minimum are based on simulating group’s intelligence and communication between the individuals in order to effectively search the space.

They are used for finding points in search space close to the optimal one (global minimum) in terms of fitness function. Very commonly fitness function describes a certain dependence (e.g. approximate solution error), which in an engineering problem should have the lowest possible value (problem of minimization). In contrast to classic mathematical methods, they have small requirements for objective function, which is their biggest advantage. Usually, within heuristic optimization algorithms, two phases of searching can be distinguished: exploration part – searching through possibly the vastest part of space, and exploitation part – looking for good quality solutions in a narrow part of searched space. Algorithms of this sort use different techniques, from simple local searching to advanced evolutionary processes. They use mechanisms preventing the method from getting stuck in a limited area of searched space (falling into a local minimum). These algorithms are independent of a specific problem (they do not depend on the fitness function). Algorithms use knowledge about the problem and/or experience accumulated in a process of searching the domain (agents “communicate” with one another), with quality not decreasing as the iteration increases. This kind of algorithms falls into the category of probabilistic algorithms, meaning that their method of work include some random elements. However, a well-tuned algorithm in a vast majority of cases should be able to provide solutions, which are close to one another.

The biggest problem of heuristic optimization algorithms is tuning them properly, according to the problem to be solved. Metaheuristic optimization algorithms can be divided into four major groups:

These algorithms gained interest because of their effectiveness in various optimization engineering problems. Examples of their usefulness include publications (Kalita et al., 2023; Alyami et al., 2024; Zhang and Chi, 2023; Brociek and Słota, 2015). In this paper, some of metaheuristic optimization algorithms were used and compared, from which, three best (in terms of solving inverse problem for model with fractional derivative) were described in further subsections.

This section is designed to present how the presented method of solving an inverse problem for the model of anomalous diffusion works. As an example, the results provided by a few chosen metaheuristic optimization algorithms were compared.

In a considered inverse problem, fractional boundary condition (4) on the right side is recreated, specifically, function ψ appearing in mentioned condition, as well as orders of fractional derivatives α, β. In the model (1)–(4), the following numeric data was assumed: x ∈ [ 0 , 1 ] , t ∈ [ 0 , 400 ] , λ ( x , t ) = 2 x t , φ ( x ) = 150 x , f ( x , t ) = 5 t 2 / 5 x 2 2 Γ ( 2 5 ) − 10 π x 3 / 2 t 10 Γ ( 9 10 ) − 2 t x ( 4 t x − 15 π t x + 150 ) π . The objective of this example is to test proposed algorithm (treated as a benchmark), hence the data sought – orders of derivatives α, β and function ψ are known and have the following values: α = 0.6 , β = 1.5 , ψ ( t ) = ( t − 10 ) 2 + 2 t ( 8 t − 45 π t + 900 ) 3 π + 50. Upon writing function ψ in a numerical form, the following was obtained: 3.00901 t 2 − 53.1736 t 3 / 2 + 339.514 t − 20 t 1 / 2 + 150. So the sought function ψ and orders of derivatives α, β were identified in a following form: (34) α = a 0 , β = a 1 , ψ ( t ) = a 2 t 2 + a 3 t 3 / 2 + a 4 t + a 5 t 1 / 2 + a 6 , where parameters a 0 , a 1 , a 2 , a 3 , a 4 , a 5 , a 6 are unknown. In a Table 1, the domain of each parameter, as well as it’s exact value are presented.

Data necessary for solving the inverse problem is value of the function u in control point x p = 1 (the right boundary). This data was generated as a result of solving the direct problem for the measurement grid 400 × 400. During the algorithm’s work on solving the inverse problem, grid of the size 100 × 200 was used. Using different sizes of grid is for evading the phenomena known as inverse crime (Kaipio and Somersalo, 2005). To find the minimum of fitness function (13), which describes the error of approximated solution, selected metaheuristic optimization algorithms were used. Tested group of algorithms includes described in the Section 5 Group Teaching Optimization Algorithm (GTOA), Artificial Bee Colony (ABC) and Jellyfish Search (JS). Apart from previously mentioned, following heuristics were used: Equilibrium Optimizer (EO), Grey Wolf Optimizer (GWO), War Strategy Optimization (WSO), Tuna Swarm Optimization (TSO) and Ant Colony Optimization (ACO). In each algorithm, values of parameters such as number of iterations and size of the population were chosen in a way that ensures similar number of calls of fitness function. This is to compare the results obtained from these algorithms.

Obtained results of the search of the minimum of fitness function are presented in a Table 2. Colour blue was used to mark three algorithms (GTOA, ABC and JS), which returned a satisfying result. Colour red was used to mark the rest of the algorithms, which did not handle the task so well. Smallest value of objective function was achieved for GTOA and was approximately 0.007. Also, in this case, errors of identifying parameters a 0 , a 1 , … , a 6 are the smallest. Similarly for ABC and JS, obtained errors and value of objective function are on acceptable level. It is also worth to pay attention to the results from ACO, where value of fitness function is low. However, the errors in identification of parameters are on an unacceptable level, which proves that the algorithm got stuck in a local minimum. For other algorithms,a the values of objective function are on a relatively highly level.

Another important indicator is fitness of data from control point (so called measurement data) to data obtained from model. On one hand, the measure determining this fit is the value of the objective function F, while on the other hand, it is the errors of the reconstructed function u for the identified parameters. Figures 6 (for GTOA and ABC) and 7 (for other algorithms) present errors of distribution of the recreated function u in a control point. For the two best algorithms (GTOA and ABC), these errors are much smaller ( ≈ 10 − 2 ) than for the rest of the algorithms ( ≈ 10 1 ) (see Figs. 6, 7).

Because the exact solution is known, we decided to calculate the errors of reconstruction function u in a full domain (not only in control point). In Fig. 8 errors of reconstruction function u for identified model parameters are presented. For three best algorithms, these errors are on an acceptable level and they share similar characteristics. It is the case of GTOA (maximum error ≈ 1.5), ABC (maximum error ≈ 10) and JS (maximum error ≈ 20). In the rest of algorithms, maximum errors are high and unacceptable.

As part of the conducted computations, a comparison of the error of identification of the function ψ occurring in the boundary condition was also performed. These errors were calculated using the formula: (35) E a b s = 1 t ∗ ∫ 0 t ∗ | ψ ( t ) − ψ approx ( t ) | d t , (36) E r e l = E a b s ( 1 t ∗ ∫ 0 t ∗ | ψ ( t ) | d t ) − 1 100 % , where ψ is the exact solution, and ψ approx is the approximate solution. The corresponding results are presented in Table 3. For the two best algorithms, the percentage errors are 3.21 % for GTOA and 6.82 % for ABC. The JS algorithm ranked third with an error of 16.21 %, while for the remaining algorithms, the errors were large ≈ 35 %.

This article focuses on the method of solving the inverse problem for diffusion equation with fractional derivatives. In the considered task, the orders of derivatives and the function occurring in the boundary condition were identified. In the presented approach, the forward problem was solved using an implicit finite difference scheme, while the inverse problem was solved using heuristic optimization algorithms. The inverse problem turned out to be difficult to solve and required identification of seven parameters.

To solve the inverse problem, the following metaheuristic algorithms were compared: GTOA, ABC, ACO, EO, GWO, JS, TSO, WSO. Satisfactory results were obtained for the Group Teaching Optimization Algorithm (GTOA) and the Artificial Bee Colony (ABC) method. Additionally, the Jellyfish Search (JS) algorithm yielded an acceptable result. The remaining algorithms proved unsuitable for this type of problem.

One of the conclusions from the conducted research and next step in the research is the possibility to build a hybrid algorithm. Firstly, a heuristic algorithm could be used for initial solution localization (exploration part), while deterministic methods such as Nelder-Mead or Hooke-Jeeves could be employed for more focused searches (eploitation part). This is the next step and a plan for future research in the work carried out by the authors of this paper.