A directed tabu search method for solving controlled Volterra integral equations

This paper proposes a pseudo-spectral scheme for obtaining approximate optimal control and state. This computational scheme represents the solution of the optimal control problem (OCP) by an mth degree Lagrange interpolating polynomial, using Legendre nodes. Then, OCP of non-linear Volterra integral equation is transformed into an optimization problem. Directed tabu search (DTS) method is utilized to derive the solutions of the optimal control and state as well as the optimal value of the objective function. In the DTS method, two neighborhood-local search strategies based on Nelder–Mead method and adaptive pattern search are applied. In addition, a tabu list with anti-cycling rules, the so-called tabu regions and semi-tabu regions are used. In addition, diversification and intensification search schemes are employed. DTS is able to converge to the global optimum solutions of a set of numerical examples. Therefore, some good balance between the diversification and intensification ensures a faster and efficient convergence to get quality solutions. Numerical examples are provided which confirm the reliability and efficiency of the proposed method. Moreover, a comparison is made with optimal solutions obtained by the other numerical methods in the literature.


Introduction
The controlled Volterra integral equations (VIE) are very momentous, because these equations arise in modelling many classes of phenomena. Many problems in economics, biology, epidemiology, and memory effects can be modeled as Volterra control problems [34]. Since the analytical methods based on the necessary conditions obtained using the Pontryagin's maximum principle and dynamic programming have less implementation ability, different numerical approaches have been devised to overcome the problems arising from the application of analytical methods [8-12, 34, 35, 40, 46].
In this paper, we have introduced a spectral approach based on collocation method to obtain optimal control of systems governed by VIE. Spectral techniques have been demonstrated to provide effective and flexible methods to solve diverse problems numerically [5-7, 44, 46, 47]. The presented method consists of reducing the OCP to a nonlinear programming (NLP) by first expanding the state and control functions in terms of Lagrange interpolating polynomial with unknown coefficients. These polynomials together with numerical integration and collocation approach are utilized to convert OCP to finite-dimensional programming problem. Since the resulted finite-dimensional programming problem is large scale, there are various optimization methods which can be implemented to solve it. It is tried that to find a method which is useable and effective on solving these problems in large-scale setting. Metaheuristic methods are becoming one of the important tools for producing robust and efficient methods that compute approximate solutions of problems with high quality in reasonable computation time. Therefore, metaheuristic methods are applied for solving real-world optimization problems [18,19,24,30,33,42,43]. Among those metaheuristic algorithms, tabu search (TS) [18] is one of the most efficient memory-based methods. However, the original TS cannot completely exert the optimization problems. The modified versions of the TS have been presented to ameliorate its performance, for example, One of the enhanced versions of the TS is referred to as directed tabu search (DTS) [26]. DTS has shown a notable performance when used to solve non-linear continuous optimization problems. Due to its neighborhood mechanisms, DTS had quick convergence, since it detected promising search directions in the neighborhood of a global minimum. Indeed, the slow convergence of TS was overcome by incorporating local search strategies, such as the Nelder-Mead (NM) method [39] and the adaptive pattern search (APS) [25], into the main algorithm. On the other hand, the DTS has been successfully applied to various real-world optimization problems [26,28,29,31,38,41]. Therefore, we survey application of DTS to solve OCP and use the values of the parameters that are required in the implementation of the DTS method of Hedar and Fukushima [26].
This paper is arranged as follows: OCP is defined in Sect. 2. In Sect. 3, we illustrate that how the pseudo-spectral method using Lagrange interpolation polynomials converts the continuous OCP to a finite-dimensional optimization problem. Directed tabu search algorithm is exerted to solve this non-linear optimization problem which has global convergence. In Sect. 4, basic convergence results are given. Section 5 contains numerical examples that demonstrate the efficiency and accuracy of the proposed framework. Section 6 ends this paper with a brief conclusion.

The formulation of optimal control problem
In the present work, we focus on a numerical approach to solve the following OCP of non-linear VIE: Problem A: Specify the real-valued continuous optimal control u Ã ðtÞ and the corresponding optimal state x Ã ðtÞ, t 2 ½0; 1, that maximize (or minimize) the functional J ðx; uÞ ¼ It is assumed that F , K and y are real valued and continuously differentiable with respect to their arguments, and both x and u belong to Sobolev space W l;1 . The interval [0, 1] can be transformed to an arbitrary interval [a, b] via a proper change of variable. It is also supposed that the optimal control of this problem is unique. Analytical discussions about the existence and uniqueness for optimal control of systems governed by non-linear VIE 2 can be found in [1,27] and references there in.

The method of discretization
Pseudo-spectral approaches are powerful tools that are frequently employed in various fields of numerical analysis. These methods have a higher order of accuracy in the case of smooth solution of any considered problem. In this section, implementation of the pseudo-spectral method on problem A is given. At first, the control and state functions are approximated in terms of Mth degree Lagrange interpolating polynomials of the form where w k ðtÞ for k ¼ 0; . . .; M À 1 are the Lagrange interpolating polynomials defined by in which s i , i ¼ 0; 1; . . .; M À 1 can be Gauss-Legendre (GL) or Gauss-Chebyshev (GC) nodes. GC nodes are the roots of Chebyshev polynomial of the first kind in ½À1; 1 [36]. These nodes in the interval [0, 1] are defined by GL nodes are the zeros of Legendre polynomial. Explicit formulas for GL nodes are not known. However, they can be computed numerically using existing subroutines [14]. The best choice of interpolation points to ensure uniform convergence is GC nodes (see Section 2.2 [36]). It can be verified that The approximation process of problem A includes the discretization of both the cost function and the controlled integral equation constraint.

Controlled integral equation discretization
We discretize Eq. 2 using shifted GC nodes, t p , p ¼ 0; 1; . . .; M À 1, as follows: Kðs; t p ; xðsÞ; uðsÞÞdt: ð7Þ We should notice that GC nodes should be transformed into the interval [0, 1]. By substituting 3 in 7, we obtain The GL quadrature formula is utilized to approximate the integral term in Eq. 8. For this purpose, this change of variable must be made with the following form: Then, by applying GL quadrature for approximating the integral, we derive and s j s are the GL nodes, in the interval ½À1; 1, and w j s are the corresponding weights. The quadrature weight, w j , can be obtained from the following relation: In Eq. 11, L M is Legendre polynomial of degree M; For more details about these polynomials, see Tohidi and Samadi [48]. Finally, the controlled Volterra integral Eq. 2 is reduced to M non-linear algebraic equations given in Eq. 10.

Cost functional discretization
For approximating the cost functional stated in Eq. 1, we utilize the GL quadrature after the proper interval transformation Z 1 0 F ðt; xðtÞ; uðtÞÞ where w 0 j ¼ 1 2 w j and s 0 j ¼ s j þ1 2 ; and s j and w j are GL nodes and weights stated in Eq. 11.

The discretized optimization problem
Finally, problem A is approximated by the following NLP: subject to Þ are the unknown parameters of our discrete problem, and p ¼ 0; 1; . . .; M À 1: As shown above, the resulted discrete problem can be reformulated as the following non-linear programming problem: subject to When the continuous problem A is discretized, the infinitedimensional problem A is reduced to the finite-dimensional non-linear optimization problem A M . Many well-developed NLP techniques can be used to solve this problem [32]. The method has been used to solve the non-linear constrained optimization problem which is based on DTS method.

Directed tabu search
DTS is one of the recent development metaheuristics proposed for combinatorial optimization problems. In the DTS method, three search procedures were used: Exploration, Diversification, and Intensification. In the Exploration Search, a local search procedure is applied to produce trial moves, based on the Nelder-Mead method (NM) and the adaptive pattern search method (APS). Besides, memory elements called Tabu Region (TR), Semi-TR, and a multiranked Tabu List (TL) were presented to obtain anti-cycling rules. Visited Regions List (VRL) was also presented for the Diversification Search to diversify the search to unvisited areas of the solution region. At last, supposing that one of the best points generated by the Exploration and Diversification Searches was close to a global solution, the Intensification Search was utilized again as the final step to rectify the best solutions visited so far. Indeed, the Exploration and Diversification search processes were congregated to prepare the DTS main loop and were repeated until the termination conditions were satisfied. Therewith, the Exploration Search process was comprised as an inner loop within the diversification loop [26].
Here, the DTS method applies pattern search to get better performance. The main operators used in the DTS are: Memory elements Memory elements are used to provide anti-cycling rules and to diversify the search to unvisited areas of the solution space.
• TL Due to the recency and objective function values, solutions are ranked and saved in the TL. • VRL The centers of the visited regions and the frequency of them are saved in the VRL for generating new diverse solutions. • TR and Semi-TR No new trial solution is allowed to be generated in a TR. A Semi-TR is a neighbor region around TR. New trial solutions are generated in a Semi-TRs, so that returning back to a visited TR is avoided.
APS In the APS method, the approximate descent direction method is used to find promising directions. Exploration search In exploration search process, trial solutions are generated based on neighborhood region structures and the APS procedure. Diversification search While the exploration search is taking a long time to improve the current solution, then the VRL is called to generate new diverse solutions. Intensification search In this phase, NM method IS applied starting from some of the best solution found till now.
In the DTS method, the exploration, diversification, and intensification search procedures are applied to develop the search process on a global strategy. Indeed, these search procedures are used to get a vast exploration and an extensive diversification. The DTS has been used in this work is stated in Algorithm 1.

1.
Initialization. Select an initial solution, and set the TL and the VRL to be empty.

Theoretical consideration
In numerical schemes, discussion on convergence of an propounded algorithm is essential. The convergence of discretization methods for OCP is a topic of active research [4,13,16,16,17,[20][21][22]. In this section, we investigate the following questions: 1. Does the discretized problem A M admit a feasible solution? 2. Does a sequence of discretized optimal solutions of problems A M converge to the optimal solution of problem A?
In Theorem 1, it has been proven that the feasibility of discretized problems A M is guaranteed if a solution to the continuous-time problem A exists? It should be noted that the convergence analysis of the proposed scheme can be done using similar idea provided in [48]. First, we need the following lemma and definition.
Definition 1 [21] A function Q : ½0; 1 ! R belongs to Sobolev space, W l;p , if its jth weak derivative, Q ðjÞ , lies in L p ½0; 1 for all 0 j l.
Lemma 1 [20] Given any function QðtÞ 2 W l;1 ; t 2 ½0; 1; there is a polynomial P M ðtÞ of degree M or less, such that where C is a constant independent of M and C ¼ kQk W l;1 .

Remark 1
The computational interval can be transformed from [0, 1] to [a, b] via an affine transformation.
Theorem 1 (Existence) [48] Let (x(t), u(t)) be a feasible solution for problem A and x(t) and u(t) belong to W l;1 with l ! 2: Then, there exists a positive integer M 1 , such that, for any M [ M 1 , problem A M has a feasible solution, ðX M ðt p Þ; U M ðt p ÞÞ and this feasible solution also satisfies in the following inequalities: in which t p is Gauss-Legendre node or Gauss-Chebyshev node, and L 1 and L 2 are positive constant independent of M.
We establish in Theorem 1 the existence of feasible solutions for discretized problem A M . The next theorem demonstrates that there exists a sequence of optimal solutions of problems A M converging to an optimal solution of problem A? Let X Ã ¼ ½x MÃ 0 ; x MÃ 1 ; . . .; x MÃ MÀ1 T and U Ã ¼ ½u MÃ 0 ; u MÃ 1 ; . . .; u MÃ MÀ1 T be the optimal solutions of the problem A M . The approximate optimal control and state are The next theorem will be demonstrated the convergence of the sequence fðX M Ã ðt p Þ; U M Ã ðt p ÞÞg 1 M¼M 1 .

Illustrative examples
In this section, two numerical examples are considered to illustrate the effectiveness of the pseudo-spectral method along with DTS scheme. The numerical results obtained from our propounded method is also compared with the results in other works in [35,48]. The following numerical implementation is performed using Mathematica. subject to Borzabadi et al. [13] xðtÞ ¼ yðtÞ þ where yðtÞ ¼ t cosðtÞ À 1 2 t 3 : The optimal value of cost function is J Ã ¼ 0. The optimal control u Ã ðtÞ and corresponding optimal state x Ã ðtÞ are as follows: x Ã ðtÞ ¼ sinðtÞ; The approximate solutions for both state and control functions together with the exact solutions are depicted in Figs. 1 and 2. The blue lines in the figures indicate the exact optimal solutions and the red lines are the approximate solutions. The absolute errors of cost functional value, E J , for example 1, are compared with those obtained from the methods Maleknejad and Almasieh [35] and Tohidi and Samadi [48] in Table 1. where yðtÞ ¼ e Àt 2 þ tð1Àe Àt 2 Þ

2
: This problem has the optimal solutions u Ã ðtÞ ¼ t and x Ã ðtÞ ¼ e Àt 2 : Figures 3 and 4 show the exact and approximate optimal control and state. Table 2 gives the results obtained from our proposed scheme and the methods in Maleknejad and Almasieh [35] and Tohidi and Samadi [48].
At the end, we answer a natural question: are there advantages of the proposed pseudo-spectral method along with DTS method compared with the existing ones? To answer this, we summarize what we have observed from numerical experiments and theoretical results as follows: • As seen in Examples 1 and 2, the proposed method with a few collocation points has satisfactory results with respect to other methods. • The proposed orthogonal collocation method leads to rapid convergence, as the number of collocation points increases. The main reason for this fast convergence is that the DTS explore the global search space using local information about promising search direction.

Conclusion
In this study, an advanced numerical PS method has been proposed for solving optimal control of Volterra integral equation by means of Lagrange polynomials via collocation method. The problem has been reduced to a finitedimensional parametric optimization, and there exist many effective algorithms which can be applied to solve the NLP. To solve this problem, directed tabu search (DTS) method is applied. This algorithm comes from the incorporation of a modified Tabu Search (TS) with two local optimizer methods. The local optimizers are direct search    methods, which aim to find a better solution neighbor to the present best solution by conducting the intensification and diversification along the move direction of the best solution in the area. The results indicate that using DTS could help to get better solutions of problems than available approximate methods. Illustrative examples have shown the validity, applicability, and efficiency of the proposed method. The method is in the case of optimal control of systems governed by VIE which is applicable in the field of practical science and engineering [34].