Reconstruction methods for AHP pairwise matrices: How reliable are they?

Applied Mathematics and Computation 279 (2016) 103–124

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Reconstruction methods for AHP pairwise matrices: How reliable are they? Marcelo Karanik a, Leonardo Wanderer a, Jose Antonio Gomez-Ruiz b,∗, Jose Ignacio Pelaez b a b

Artificial Intelligence Research Group, National Technological University, 3500 Resistencia, Argentina Department of Computer Science and Artificial Intelligence, University of Málaga, 29071 Málaga, Spain

a r t i c l e

i n f o

Keywords: AHP Pairwise matrix reconstruction methods Decision support systems Ranking comparison Genetic Algorithms Firefly Algorithm

a b s t r a c t Habitually, decision-makers are exposed to situations that require a lot of knowledge and expertise. Therefore, they need tools to help them choose the best possible alternatives. Analytic hierarchical process (AHP) is one of those tools and it is widely used in many fields. While the use of AHP is very simple, there is a situation that becomes complex: the consistency of the pairwise matrices. In order to obtain the consistent pairwise matrix from the inconsistent one, reconstruction methods can be used, but they cannot guarantee getting the right matrix according to the judgments of the decision maker. This situation does not allow proper evaluation of methods reliability, i.e. it is not possible to obtain a reliable ranking of alternatives based on an inconsistent matrix. In this work, a new way to evaluate the reliability of matrix reconstruction methods is proposed. This technique uses a novel measure for alternatives ranking comparison (based on element positions and distances), which is introduced in order to compare several matrix reconstruction methods. Finally, in order to demonstrate the extensibility of this new reliability measure, two reconstruction methods based on bio-inspired models (a Genetic Algorithm and the Firefly Algorithm) are presented and evaluated by using the aforementioned reliability measure. © 2016 Elsevier Inc. All rights reserved.

1. Introduction Decision making (DM) process involves many aspects to consider simultaneously and, consequently, an efficient decision maker needs cognitive and technical skills in order to do optimal evaluations. Expert decision makers have more probabilities to choose the best option to solve a problem, but they are not infallible. Nowadays, decision support systems (DSS) [1,2] are used to assist the decision maker and, accordingly, to provide more self-confidence for decision tasks. There are many types of DSS and they are used with personal through managerial and enterprise purposes [3–5], assisting either individual or group decisions [6–8] and implementing diverse techniques [9] on stand-alone and web-based architectures [10–12]. In general, a decision problem involves selecting from several alternatives based on multiple criteria. Multi criteria decision analysis (MCDA) methods share certain components: a finite set of alternatives, at least two criteria and a decision maker. MCDA can be classified as outranking methods [13,14], multi attribute utility theory methods [15,16] and non-classical

∗

Corresponding author. Tel.: +34 951952391; fax: +34 952131397. E-mail address: [email protected] (J.A. Gomez-Ruiz).

http://dx.doi.org/10.1016/j.amc.2016.01.008 0096-3003/© 2016 Elsevier Inc. All rights reserved.

104

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

methods, all of which use many environment representations such as hierarchical structures, fuzzy expressions, linguistic terminology, etc. Analytic hierarchy process (AHP) [15], developed by Thomas Saaty, is one of the most widely used MCDA methods. Its usage is highly extended in several areas, such as: technological investment evaluation [17], analysis of financial parameters [18], strategic planning [19], logistic [20], supplier selection [21–24], inventory classification [25], Internet access technology [26], IT project selection [27] and even for reengineering of the health-care system [28]. One important drawback of AHP is the large amount of information needed to complete the pairwise matrices. For n elements to be compared, the number of preference values is (n2 − n )/2. Clearly, this number grows exponentially for high values of n and, under this circumstance, the decision maker should state inconsistent judgments. Saaty [15] proposes a multiplicative reciprocity among the elements of pairwise matrices in order to compute the consistency ratio (CR) of AHP pairwise matrices (see Section 2) and he establishes that a pairwise matrix is consistent if its CR is under 0.1. There are two issues related to the CR: (a) for computing it, the pairwise matrix must be complete, that is without missing values, and (b) if its value is over 0.1 some matrix values must be modified. In order to assist the decision maker, several methods have been proposed to complete missing values [29–32] and to improve the CR [29,30,33,34]. That is, the decision maker fills some (or all) matrix values and, after that, a reconstruction method to complete matrix or to improve CR (or both at the same time) is used. Obviously, this makes the valuation process easy, but here arises a key question: who really makes the final decision, the decision maker or the reconstruction method? Even more, there are several studies about the relation between the CR and the final priority vector [35,36] and they indicate that some expert judgments are not always appropriately reflected in the final ranking. Even though the reconstruction methods have been compared with each other, the main goal has always been to reduce the CR until reaching acceptable value beyond establishing whether the obtained results are reliable or not. Clearly, this situation is very complex since there is no way to establish the correct configuration from an incomplete matrix or even, from a complete but inconsistent one. In other words, given several pairwise matrix configurations obtained by means of different reconstruction methods, there are no reference elements (beyond the consistency ratio) to determine which is more adequate with respect to the expert’s preferences. This article has two goals: first, to propose an original approach to determine the confidence of AHP matrix reconstruction methods, and second, to implement a new reconstruction method based on two bio-inspired algorithms (a Genetic Algorithm and the Firefly Algorithm) in order to complete AHP pairwise matrices improving their CR simultaneously. Additionally, a new technique to compare alternatives rankings is explained. The usage of this technique can be extended to any ranking-based DSS. This paper has been structured as follows: AHP is briefly described in Section 2; in Section 3 a new way to determine the reliability of reconstruction methods and the proposed ranking comparison technique are presented; in Section 4 three reconstruction methods for improving pairwise matrix consistency and three methods for reconstructing incomplete pairwise matrices are compared using the proposed technique; in Section 5 two new reconstruction methods based on a Genetic Algorithm and Firefly Algorithm are implemented and tested; and finally in Section 6 the conclusions about reliability, usability and extensibility treated in this work are presented. 2. The analytic hierarchy process (AHP) 2.1. AHP generalities The analytic hierarchical process helps the decision maker to choose from alternatives of solution based on the analysis of criteria and alternatives using pairwise comparisons. In AHP, the problem is represented using a hierarchy which is an abstraction of the whole problem. The entities of the same level are not connected with each other but they are fully connected with entities of adjacent levels. These relations determine the hierarchy structure of criteria and alternatives to evaluate. The basic idea of AHP is to arrange goals, attributes and issues in a hierarchical structure to achieve two objectives: to obtain a complete view of relationships inherent to the situation and to provide a comparison mechanism at each level that uses the same order of magnitude. One way to summarize judgments is to take a pair of elements and compare them with a single property without considering any other elements [37]. The pairwise comparisons determine the priorities of a set of elements (criteria or alternatives) and are made by means of a value scale. The Fundamental Scale of Saaty [15] is made up of values 1–9 and their multiplicative reciprocals: equal (1), moderately more (3), strongly more (5), very strongly more (7), extremely more (9) and compromise values (2, 4, 6, and 8) among them. These values are put into the superior triangle portion of the pairwise comparison matrix and the inferior triangle portion is completed with the reciprocal judgment values, in order to obtain a reciprocal matrix (reciprocal axiom [38]). Saaty defines three axioms in addition to the reciprocal axiom [38]: the homogeneity axiom (to place in each level the characteristics that can be compared); the dependence axiom (to establish the connections among levels to create a hierarchy of external dependences among them); the expectation axiom (to give a complete vision of the problem). These four axioms govern AHP and use the notion of paired comparisons as the premise.

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

105

Formally, given the reciprocal n × n pairwise matrix A:

⎡

1 ⎢ .. ⎢ . ⎢ ⎢ ai1 A=⎢ ⎢ ... ⎢ ⎢ . ⎣ .. an1

a12 .. . ai2 .. . .. . an2

··· ···

···

a1 j .. . ai j .. . .. . an j

··· ···

···

⎤

a1n .. ⎥ . ⎥ ⎥ ain ⎥ .. ⎥ ⎥ . ⎥ .. ⎥ ⎦ . 1

(2.1)

where

aii = 1 and a ji = 1/ai j for 1 ≤ i,

j≤n

are values of the Fundamental Scale of Saaty [15] and the eigenvector associated to the largest eigenvalue of matrix A:

w = [w1 , . . . , wi , w j , . . . , wn ]T

(2.2)

is the priority vector w associated to matrix A. In this way, the priority vector from criteria pairwise matrices are obtained by calculating the principal eigenvector (using, for example, the power method [39,40]). Then, the eigenvector values must be normalized before being used in the next level. Thereafter, the pairwise matrices for alternatives level are filled in and the principal eigenvector is calculated and normalized for each matrix. These normalized eigenvectors associated to an alternative pairwise matrix form a column in a new matrix used to determine the final values of priorities for each alternative. This matrix and the normalized eigenvector associated to the criteria pairwise matrix are multiplied, obtaining in this way the final ranking of alternatives. 2.2. AHP matrix consistency In AHP, it is necessary to assure that the decision maker judgments are consistent. There are several mathematical ways to measure the consistency of a pairwise reciprocal matrix. For instance, Saaty [15] proposes a multiplicative reciprocity among the elements of pairwise matrices while other authors define a consistency ratio based on percentiles [33], or analyze the matrix consistency based on the characteristic polynomial of a positive reciprocal matrix [32]. In this paper the consistency ratio based on the multiplicative reciprocity is used. This ratio is frequently used for testing the consistency on pairwise matrices reconstruction [30,41]. Before any calculation, the pairwise matrices must be complete and consistent. Actually, since it is hard to obtain a fully consistent matrix, this restriction has a tolerance degree. Saaty defines a method to determine the consistency within a matrix, using the CR. In a reciprocal positive matrix, the largest eigenvalue is greater or equal to its order. Saaty demonstrates [15] that a reciprocal matrix A is perfectly consistent if its largest eigenvalue is equal to the matrix order and defines the consistency index (CI) as:

CI =

λmax − n

(2.3)

n−1

where λmax is the largest eigenvalue of matrix A and n its order. This value is compared with an average over a large number of reciprocal random matrices of the same order, obtaining the consistency ratio (CR):

CR =

CI RI

(2.4)

where RI is the random index obtained from 500 positive reciprocal matrices randomly created with Saaty’s scale values (see Table 2.1). If the CR is small, about 10% or less, the principal eigenvector (or the priority vector), w, associated to the largest eigenvalue of A is accepted. In other words, to determine if a pairwise matrix is consistent, its CR index must be under 0.1 [15]. This ratio represents approximately 25% of matrices with order 3 but this percentage decreases drastically when the matrix order grows [42]. There are several methods to improve the CR index when its value is over 0.1 but it is necessary that the matrix is complete. The final matrix configuration will depend on the reconstruction method used to improve it or to complete it. Undoubtedly, the reliability of the results returned by these methods is crucial in order to help adequately in the decision making process. Table 2.1 Values of RI. n

3

4

5

6

7

8

9

RI

0.58

0.9

1.12

1.24

1.32

1.41

1.45

106

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

3. Ranking-based reliability measure As discussed before, in AHP it is necessary to assure that the decision maker judgments are consistent. Suppose that a decision maker expresses “I prefer A much better than B and I strongly prefer B better than C”. If he also says “I prefer C over A”, this last judgment is inconsistent respect to the previous two and it must be corrected. Now, the question that arises is: which judgment is wrong? Well, the answer is simple: it depends on the method used for the matrix reconstruction because different reconstruction methods would provide different results. That is, different consistent pairwise matrices would be obtained and different priority vectors would derivate from these matrices. The key factor here is the final ranking obtained from an inconsistent matrix. Certainly, this ranking is wrong due to the fact that it is based on inconsistent judgments and there is no way to determine what the correct solution for an incomplete or inconsistent matrix is. This occurs because there are many combinations to complete or improve it, and all might be considered adequate. Hence, all reconstruction methods return valid solutions as long as they return a consistent matrix. In short, several methods might return different rankings and they will be adequate because there is no way to contrast them against the correct ranking. Apparently, this problem seems to have no solution because it is not possible to obtain a correct priority vector from an inconsistent pairwise matrix. However, what would happen if the reconstruction methods are tested using consistent matrices (CR ≤ 0.1)? In this way it is possible to test the obtained ranking after reconstruction and evaluate it against the original one (the consistent one). Consequently, the methods that improve the CR maintaining the same priority vector can be considered more reliable. Putting into practice the idea expressed above requires a method to measure how different two priority vectors (rankings) are or, more generally, how different two vectors are. Suppose there are two priority vectors: one obtained from an original consistent AHP matrix, the original priority vector or OPV, and one obtained from this same original consistent AHP matrix but enhanced by some matrix reconstruction method, the enhanced priority vector or EPV. Both mentioned matrices have the same order (n) and, consequently, the two priority vectors have n elements (obviously they are the ranked alternatives). Here, two different situations can occur: either both rankings are equal or they are not. The first case is not problematic, the reconstruction method has improved the AHP pairwise matrix and the alternatives order is equal to the original consistent matrix. The second case is more complex because it requires the establishment of priority vectors similarity. There are several possibilities to determine the similarity of two rankings. One way is considering, position by position, all elements (or alternatives) into the priority vector. Whether any alternative occupies different positions in two rankings, evidently both are different. Clearly, this method is extremely restrictive and does not contain much information about the difference. Another way is to consider only the top-leading alternative, since it is the most appropriate solution to face the decision problem. This approach is based on the supposition that the best alternative is enough to solve the decision problem. Therefore, in this way, if a technique is good to keep the same alternative in the first place, no more exhaustive analysis is necessary. Besides the aforementioned points, it is certainly desirable to have a good measure for the rankings for comparing priority vector-based DSS. In this sense, to measure different levels of errors using a unique measure considering the whole set of compared alternatives in the rankings would be ideal. For that, the similarity within the rankings must be based on the position of the ranked alternatives. If an alternative does not change its position in both rankings there is no error. Otherwise, an error which is proportional to the distance between those positions is computed. Additionally, the error value must be increased if it involves alternatives on the top of original priority vector. This increment is based on the fact that decisions are taken considering the likelihood of each alternative to solve the problem, and with higher likelihood comes a better position into the ranking for that alternative. In order to explain our proposed measure, the next example is introduced and used through this section: Let A = {a1 , a2 , a3 , a4 } the set of alternatives associated to a decision making problem. Let OPV = [a1 , a3 , a4 , a2 ] the original priority vector associated to a consistent AHP pairwise matrix and EPV = [a2 , a1 , a4 , a3 ] the enhanced priority vector computed by a matrix reconstruction method. To establish the difference between OPV and EPV (both vectors of n order, and for the example n = 4) the distance value for every alternative must be computed. To do this, the expression used to compute the difference between the positions in both rankings for that same alternative ai (with i = 1, 2, ..., n) is:

dOPV,EPV (ai ) = | posOPV (ai ) − posEPV (ai )|

(3.1)

where dOPV,EPV (ai ) is the distance vector for all alternatives ai , posOPV (ai ) and posEPV (ai ) are the positions of the alternative ai within the rankings OPV and EPV respectively. For our example, the position and distance values for all alternatives are shown in Table 3.1. Notice that the ordered list in column 1 of Table 3.1 does not coincide with the OPV alternatives order. For notation reasons, the distance vector is ordered using the OPV alternatives order (this same criterion is used in the rest of tables of this section). Table 3.2 shows this modification. As different positions have different relevance or importance in the decision problem, a weight set needs to be associated to the positions of the alternatives. In this sense a normalized weight vector wm is defined and it contains the importance

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

107

Table 3.1 Distance vector for each alternative in OPV and EPV. Alternatives

posOPV (ai )

posEPV (ai )

dOPV,EPV (ai )

a1 a2 a3 a4

1 4 2 3

2 1 4 3

1 3 2 0

Table 3.2 Distance vector ordered using the OPV alternatives’ order. Alternatives in OPV

posOPV (ai )

posEPV (ai )

dOPV,EPV (ai )

a1 a3 a4 a2

1 2 3 4

2 4 3 1

1 2 0 3

Table 3.3 Weight system with exponential increment. p

w2 ( p )

1 2 3 4 n

0.5333 0.2666 0.1333 0.0666

p=1

w2 ( p )

1

Table 3.4 Terms needed to calculate the Absolute Divergence Ratio. dOPV,EPV (ai )

w2 ( posOPV (ai ))

w2 ( posOPV (ai )) · dOPV,EPV (ai )

1 2 0 3 ADR2 (OPV, EPV )

0.533 0.266 0.133 0.066

0.533 0.532 0.000 0.198 1.263

value assigned to every position in the ranking (remember, first places are more important than last ones). m is a natural number that represents the multiplicative factor that indicates how many times a position p is more important than the next one. Then, wm ( p) contains the importance value assigned to each position p in the ranking, where 1 ≤ p ≤ n and n p=1 wm ( p) = 1. Clearly, if the top elements must be invariant, a higher weight must be assigned when a change is made in those positions. Consequently, a position near the top of the ranking has more weight than those far away from the top. This relation holds for all positions, i.e., wm (1 ) > wm (2 ) > · · · > wm (n ). One way to achieve these requirements is to define wm ( p) exponentially:

m−p wm ( p) = n −q q=1 m

(3.2)

where m > 1. For the previous example, we consider m = 2 (as shown in Table 3.3). In this way, each position duplicates the importance respect to the next one. With the distance function defined in Eq. (3.1) and the weight vector defined in Eq. (3.2), the Absolute Divergence Ratio (ADRm ) can be defined as the divergence between the rankings using the following expression:

ADRm (OPV, EPV ) =

n

[wm ( posOPV (ai )) · dOPV,EPV (ai )]

(3.3)

i=1

where ADRm (OPV, EPV ) is the Absolute Divergence Ratio between OPV and EPV, n is the number of alternatives, wm ( posOPV (ai )) is the weight for the position of the alternative ai in OPV, and dOPV,EPV (ai ) is the distance for an alternative ai in both vectors. Since the reference ranking is OPV, the ADRm function considers this ranking as the goal that a technique is trying to achieve. The ADRm terms for our example are calculated and shown in Table 3.4.

108

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

Fig. 3.1. Process to generate E array for four alternatives.

The final value can be summarized in the following calculus:

ADR2 (OPV, EPV ) =

4

[w2 ( posOPV (ai )) · dOPV,EPV (ai )]

i=1

= 0.533 · 1 + 0.266 · 2 + 0.133 · 0 + 0.066 · 3 = 0.533 + 0.532 + 0.000 + 0.198 = 1.263 Therefore, this value can be useful to compare the rankings produced by two different techniques against the same OPV and see how well they perform. That is, given two matrix reconstruction methods it is possible to determine which returns the EPV more similar to OPV, and consequently, which one is more reliable. To do this the Absolute Divergence Ratio must be calculated for both EPVs (each one obtained from the respective matrix reconstruction method) and the minor ADRm (OPV, EPV ) value corresponds to the more reliable method (that is, the one whose EPV is more similar to the OPV). Even though the adequate mechanism provides for the Absolute Divergence Ratio, the value itself does not offer an intrinsic estimation of performance. In order to improve the technique and to obtain a single and intrinsic value representing the difference of EPV regarding OPV, it is necessary to compare the obtained ADRm (OPV, EPV ) against the worst possible case of divergence, i.e., the ADRm (OPV, WPV ) where WPV is the worst priority vector. The following expression resembles this comparison and returns the poor performance of certain techniques:

RDRm (OPV, EPV ) =

ADRm (OPV, EPV ) MDRm (n )

(3.4)

where RDRm (OPV, EPV ) is the Relative Divergence Ratio between the rankings and MDRm (n ) is the Maximum Divergence Ratio for n alternatives computed by:

MDRm (n ) = ADRm (OPV, WPV )

(3.5)

Notice that MDRm (n ) only depends on n and m. That is, for n alternatives it is possible to find the most dissimilar vector in terms of distance, by using the multiplicative factor m. To do this, a procedure that computes the maximum distances for all alternatives is necessary (i.e. the worst priority vector of order n). The procedure for getting these distances is iterative and uses two arrays, S and E, which stand for Start and End. The first one contains all the elements and the second one is empty. In order to simplify the explanation, S contains all alternatives ordered, but this condition is not obligatory (clearly, the calculus is referred to the position occupied by every alternative, not to the alternative itself). The iterative mechanism to find the Maximum Divergence Ratio is shown in Algorithm 3.1. Algorithm 3.1 works as follows: each iteration takes the topmost element in S (which is not yet placed in E), and places it in E in the farthest available position (lines 03 to 14). In case of two equally distant positions being available, the higher index in E is selected, making room in the top available positions of E for the bottom alternatives in S (this is indicated by the use of ≥ conditional in line 07). Lines 15–22 compute the weight vector using the multiplicative factor m. MDRm is computed in lines 23–26 and, finally, its value is returned in line 27. Fig. 3.1 shows the process to obtain E for four alternatives. In Step 0, the S vector is formed with all the alternatives (the order being irrelevant), and the E vector is empty. In Step 1, the top element in S, a1 , is placed in the farthest position in E (position 4). In Step 2, the second element is taken, and it has two alternatives: the first or the third position. In order to allow the next elements to go in farther positions, a2 is placed in the third position. In Step 3, a3 is placed at the top, which is the farthest position, and finally, in Step 4 the same happens to a4 . After E is completed, the same calculus for the MDRm (n ) is made using the maximum distances. The MDRm (n ) terms for S and E are calculated and shown in Table 3.5.

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

109

Algorithm 3.1 Maximum Divergence Ratio generation. Input: vector order n; Multiplicative factor m; ordered Start vector S = [a1 , a2 , ..., an ]; empty vectors E, w, dS,E ; Output:MDRm (n ); Process: 01. sIndex := 1; 02. occupied := 0 03. For sIndex := 1 To n 04. maxDist := 0; 05. For current: = n To 1 06. curDist := abs(sIndex − current ); 07. If (actDist ≥ maxDist ) And (E(current ) = = nil ) Then; 08. maxDist := curDist; 09. maxIndex := current; 10. EndIf 11. E(maxIndex ) := S(sIndex ); 12. dS,E (sIndex ) := maxDist; 13. EndFor 14. EndFor 15. sum := 0 16. For i: = 1 To n val := m−i ; 17. 18. sum := sum + val; 19. EndFor 20. For i: = 1 To n 21. w(i ) := m−i /sum; 22. EndFor 23. MDRm := 0 24. For i: = 1 To n 25. MDRm := MDRm + w(i ) × dS,E (i ); 26. EndFor 27. Return (MDRm );

Table 3.5 Terms needed to calculate Maximum Divergence Ratio. Alternatives in S

posS (ai )

posE (ai )

dS,E (ai )

w2 ( posS (ai ))

w2 ( posS (ai )) · dS,E (ai )

a1 a2 a3 a4 MDR2 (S, E )

1 2 3 4

4 3 1 2

3 1 2 2

0.533 0.266 0.133 0.066

1.599 0.266 0.266 0.132 2.263

The final value can be summarized in the following calculus:

MDR2 (4 ) = ADR2 (S, E ) =

4

[w2 ( posS (ai )) · dS,E (ai )]

i=1

= 0.533 · 3 + 0.266 · 1 + 0.133 · 2 + 0.066 · 2 = 1.599 + 0.266 + 0.266 + 0.132 = 2.263 Summarizing, the Relative Divergence Ratio for our example is:

RDR2 (OPV, EPV ) =

ADR2 (OPV, EPV ) 1.263 = 0.558 = 55.8% = MDR2 (4 ) 2.263

This value points out that the ranking EPV is 55.8% different from OPV, which is a clear consequence of the dissimilarities of the two top elements. Now, this value can be used to evaluate the intrinsic performance of matrix reconstruction methods. The RDRm is a helpful metric to establish priority vectors difference and, in consequence, how reliable the matrix reconstruction methods are. Also, in order to demonstrate the flexibility of this metric, several cases will be analyzed. With this aim, the same OPV of our example will be used to compare several EPVs. In addition, three different weight vectors are used:

w2 = [0.533, 0.266, 0.133, 0.066]; we = [0.644, 0.237, 0.087, 0.032]; w10 = [0.900, 0.090, 0.009, 0.001]. where the subindex indicates the m value of Eq. (3.2) (specifically, m equal to 2, e and 10 respectively).

110

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124 Table 3.6 Relative Divergence Ratio example with moderate different rankings. Alternatives in OPV

posOPV (ai )

posEPV (ai )

dOPV,EPV (ai )

w2 ( posOPV (ai ))

we ( posOPV (ai ))

w10 ( posOPV (ai ))

a1 a3 a4 a2 ADRm (OPV, EPV ) MDRm (4 ) RDRm (OPV, EPV )

1 2 3 4

2 4 3 1

1 2 0 3

0.533 0.266 0.133 0.066 1.263 2.263 55.8%

0.644 0.237 0.087 0.032 1.214 2.407 50.4%

0.900 0.090 0.009 0.001 1.083 2.810 38.5%

Table 3.7 Relative Divergence Ratio example with extremely different rankings. Alternatives in OPV

posOPV (ai )

posEPV (ai )

dOPV,EPV (ai )

w2 ( posOPV (ai ))

we ( posOPV (ai ))

w10 ( posOPV (ai ))

a1 a3 a4 a2 ADRm (OPV, EPV ) MDRm (4 ) RDRm (OPV, EPV )

1 2 3 4

4 2 3 1

3 0 0 3

0.533 0.266 0.133 0.066 1.797 2.263 79.4%

0.644 0.237 0.087 0.032 2.028 2.407 84.3%

0.900 0.090 0.009 0.001 2.703 2.810 96.2%

Table 3.8 Relative Divergence Ratio example with slightly different rankings. Alternatives in OPV

posOPV (ai )

posEPV (ai )

dOPV,EPV (ai )

w2 ( posOPV (ai ))

we ( posOPV (ai ))

w10 ( posOPV (ai ))

a1 a3 a4 a2 ADRm (OPV, EPV ) MDRm (4 ) RDRm (OPV, EPV )

1 2 3 4

1 2 4 3

0 0 1 1

0.533 0.266 0.133 0.066 0.087 2.263 8.8%

0.644 0.237 0.087 0.032 0.049 2.407 4.9%

0.900 0.090 0.009 0.001 0.010 2.810 0.4%

Under these conditions three cases are analyzed: Case 1. Rankings moderately different, the top element in OPV is not moved into a far position in EPV:

OPV = [a1 , a3 , a4 , a2 ]; EPV = [a2 , a1 , a4 , a3 ]; the results for the metrics are shown in Table 3.6. Case 2. Rankings extremely different, the top elements in OPV is moved into far positions in EPV:

OPV = [a1 , a3 , a4 , a2 ]; EPV = [a2 , a3 , a4 , a1 ]; the results can be seen in Table 3.7. Case 3. Rankings slightly different, the top elements in OPV are in the same positions in EPV:

OPV = [a1 , a3 , a4 , a2 ]; EPV = [a1 , a3 , a2 , a4 ]; the Relative Divergence Ratio achieves its smallest possible value. See Table 3.8. Clearly, it is shown that the ADRm and, consequently, RDRm are extremely affected by the selected weight vector, since it highlights the importance of the elements within the ranking. As more importance is assigned to the top-positioned elements (using a higher value for m), higher difference between the compared rankings exists. Notice that the use of RDRm can be extended to any ranking-based DSS method. In general, it is only necessary to compare the ranking returned by the DSS technique against the reference ranking obtained from any valid strategy (for the case addressed in this section consistent AHP matrices). 4. Applying the reliability measure on matrices’ reconstruction techniques In this section several reconstruction methods for AHP matrices are used to test the proposed reliability measure. The key idea is to improve consistent AHP matrices using diverse reconstruction methods and then to compare the original

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

111

Table 4.1 Comparison of consistency improvement techniques.

Matrix order

Original CR M1 CR

5 6 7 8 9

0.0986 0.0850 0.0779 0.0647 0.0591

0.0813 0.0782 0.0742 0.0629 0.0587

M2 CR M3 CR

M1 Ranking Coincid.

M2 Ranking Coincid.

M3 Ranking Coincid.

M1 Top Coincid.

M2 Top Coincid.

M3 Top Coincid.

M1 RDR2

M2 RDR2

M3 RDR2

0.0101 0.0121 0.0133 0.0129 0.0143

781 830 860 856 895

306 280 260 247 270

986 977 956 955 911

861 910 930 927 947

593 651 679 681 736

995 991 988 982 967

08.88% 05.85% 04.42% 04.19% 02.64%

28.15% 24.73% 21.54% 19.50% 15.83%

00.41% 00.54% 00.81% 01.19% 01.84%

0.0987 0.0851 0.0781 0.0648 0.0593

priority vector (OPV) against the obtained enhanced priority vectors (EPV). In order to show the proposed measure usability, firstly only reconstruction methods used for complete matrices are tested and after that, other methods used for incomplete matrices are checked. Notice that, with the aim to maintain the nomenclature of the described reconstruction methods, the OPV and EPV are indicated as w and w respectively. 4.1. Techniques for improving complete AHP matrices consistency The inconsistency of the pairwise matrix A grows when the elements ai j are deviate from wi /w j . For this reason, improving consistency implies the reduction of this deviation. There are two questions related to the replacement of matrix elements: (a) to detect the position that most inconsistency generates and (b) to suggest the best value to put in this place. Saaty describes three methods used for detecting the most inconsistent value in a pairwise matrix. The first one [34], called M1 for nomenclature issues, determines the variation of λmax with respect to ai j :

∂λmax = vi w j − a2i j v j wi ∂ ai j

(4.1)

where vi , v j are elements of eigenvector v associated to maximum eigenvalue of AT and wi , w j are elements of eigenvector w associated to maximum eigenvalue of the matrix A. For each ai j Eq. (4.1) is computed and the highest variation of that value indicates the position (i, j ) that produces more inconsistency in the matrix. The second method, M2 , attempts to establish the perturbation [34] of each matrix element ai j . That is, given the priority vector, matrix W = (wi /w j ) and the perturbation matrix, E = (εi j ), matrix A is defined by the Hadamar product:

A=W◦E

(4.2)

where each element of matrix A is defined by:

ai j = (wi /w j ) · εi j

(4.3)

Using Eq. (4.3) every perturbation εi j is calculated by:

εi j = ai j · (w j /wi )

(4.4)

If there is no perturbation then εi j = 1. For this reason, the ai j that generates more inconsistency is the one that has a related εi j which moves away from the value of 1. The third method [31], M3 , consists of replacing each element ai j (and its reciprocal a ji ) for zeros and put value 2 on their corresponding entries over main matrix diagonal (i.e. for replacing ai j and a ji ,“2” is put in aii and a j j ). Then, λmax is calculated for the modified matrix. This procedure is repeated for the (n 2 − n )/2 elements of the superior triangular portion of matrix A. The position to be modified is one that, by replacing zero value, generates the highest λmax . Any method explained above can be used for detecting the most inconsistent value in pairwise matrices, but none suggest the best value to replace it. To do that, Saaty proposes the detection of the most inconsistent value by using any method, and after that, to apply the third method [31] replacing ai j with zeros and its reciprocal and setting the corresponding main diagonal entries with value 2. Then the priority vector w = (w i ) is calculated using the modified matrix A and ai j is replaced with ai j which is a value of the Fundamental Scale proposed by Saaty closest to wi /wj (obviously a ji is replaced with 1/ai j ). Clearly, in order to implement any of the methods the original matrix A should have no missing values.

4.2. Measuring the reliability of consistency improvement techniques To test the reliability measure proposed in Section 3, one thousand of consistent pairwise matrices of each order (from five to nine) were used, and the mean CR value was calculated (original CR column in Table 4.1). Then, these matrices were reconstructed using the three methods mentioned above (M1 , M2 and M3 ) for detection of position that generates more inconsistency and the replacing technique described. Table 4.1 shows the mean CR obtained after reconstruction (M1 CR, M2 CR and M3 CR columns) for each method and each order matrix, the times that all ranking coincide after reconstruction

112

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

(M1 Ranking Coincid., M2 Ranking Coincid., and M3 Ranking Coincid. columns), the times that the top alternative coincides (M1 Top Coincid., M2 Top Coincid., and M3 Top Coincid. columns) and, finally, the Relative Divergence Ratio using the multiplicative factor m = 2 (M1 RDR2 , M2 RDR2 , and M3 RDR2 columns). In Table 4.1 it can be observed that M2 obtains the best consistency improvement, M1 reduces the matrix CR slightly and M3 is the worst option to improve matrix consistency (even inconsistency grows). Under this criterion (inconsistency reduction), clearly, M2 is the best method. Nevertheless, the times that the ranking and the top alternative coincide between OPV and EPV is significantly minor for this best method mentioned. Under ranking modification criterion, M3 is the best and M2 has much better results than M1 . These ranking differences are reflected in each RDR2 value. It can be observed that each divergence percentage evidences the ranking mismatches mentioned. Notice that the RDR2 percentage for M2 is relatively low; this occurs because the multiplicative factor (m = 2) implies that top alternative coincidence has more importance. In case that m tends to 1 the ranking coincidence criterion acquires more relevance. Now, here is the discussion point: what is the best reconstruction method? Evidently, a good balance between inconsistency reduction and original priority vector is necessary. If a method is good to reduce the matrix inconsistency but the EPV obtained does not reflect the expert’s opinion, it is not reliable. On the other hand, if a method is good at maintaining the original priorities but it does not reduce the matrix inconsistency, it is not efficient. Under the exposed circumstances for the case analyzed in this section, the best method is M2 . 4.3. Techniques for improving consistency of incomplete AHP matrices In Section 4.1, techniques for detecting and changing the value that produces the most inconsistency were evaluated using the proposed reliability measure. As seen, those methods require complete matrices to be used. There are others that work over incomplete matrices and the mentioned measure can be used to evaluate their reliability. Next, three pairwise matrix reconstruction methods are explained and evaluated using incomplete matrices: the connecting paths method (CP), the revised geometric mean method (RGM) and the multi-layer perceptron reconstruction method (MLP). 4.3.1. Connecting paths method (CP) The missing pairwise matrix values are judgments which should be related to other non-missing values in such a way that the matrix consistency is maintained or improved. A connecting path for the element ai j in a reciprocal multiplicative matrix is given by:

CP(ai j ) = aic1 × ac1 c2 × · · · × ack j

(4.5)

where c1 , c2 , . . . , ck are the matrix sub-indexes of elements that connect i with j. An elementary connecting path is defined as follows:

CP(ai j ) = aic1 × ac1 j

(4.6)

If a reciprocal matrix element, ai j , is missing, it is reasonable to think that its value can be calculated by using the average of all connecting paths from i to j. Aczél and Saaty [43] have demonstrated that the geometric mean must be used in judgment groups in order to preserve the reciprocal property. Harker [44] established that, given a set of incomplete comparisons, the missing matrix elements on the top triangular portion of the matrix A are found by taking the geometric mean of intensities of all elementary paths connecting the two attributes in A, that is:



ai j =

p 

1/p



CP ai j

k

(4.7)

k=1

where p is the number of elementary connecting paths. Finally, the lower triangular portion of A is calculated with the following reciprocal property:

a ji =

1 ai j

(4.8)

and w is obtained from this reconstructed matrix. 4.3.2. Revised geometric mean method (RGM) The revised geometric mean method [31] is used to fill in the missing elements in a reciprocal multiplicative matrix A. The process consists of adding one to the main diagonal for each missing element in the row. Then, missing elements are replaced with zeros, that is, the elements ri j of the transformed matrix, R, are calculated as follows:



ri j =

ai j 0 1 + ci

if ai j is not a missing value, with i = j if ai j is a missing value, with i = j with i = j

where ci is the number of missing values in row i.

(4.9)

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

113

Fig. 4.1. MLP matrix reconstruction process [29].

After that, the eigenvector associated to the largest eigenvalue of matrix R is calculated and every missing element of A is replaced with the corresponding ratio. That is to say, if wR = [wR1 , wR2 , . . . , wRn ] is the eigenvector associated to the largest eigenvalue of matrix R, then each missing element in matrix A, ai j , is filled using:

ai j =

wRi wR j

(4.10)

then, reciprocal value is computed using Eq. (4.8) and, finally, the w is obtained from this reconstructed matrix.

4.3.3. Multi-layer perceptron reconstruction method (MLP) Artificial neural networks (ANN) are widely used for many and diverse tasks such as function approximation problems, pattern recognition and clustering (what is more, the ANN have been used to represent hierarchical preferences [45]). Also, ANN have been used in AHP matrix reconstruction methods [29,41]. Both implement a multi-layer perceptron (MLP) that uses the backpropagation algorithm [46]. An MLP is a feed-forward neural network in which all neurons of each layer are connected with all neurons of the next layer, starting at the input layer and finishing at the output layer through one or several hidden layers. Gomez-Ruiz et al. [29] propose an MLP architecture that receives the top triangular portion of a pairwise matrix, including its missing elements, and returns the estimate missing and non-missing elements. The number of input and output neurons is equal to the number of elements of the top triangular portion of the matrix (for n order matrix, the MLP has (n2 − n )/2 inputs and outputs neurons). This MLP learns distinct configurations of consistent matrices and it can relate an incomplete pairwise matrix with the most similar to them. To reach this objective, the MLP network uses random reciprocal pairwise matrices with low CR (under 0.1) for the training process. When this process has finished, the network is ready to recognize incomplete pairwise matrices. An interesting issue proposed in this method is that it allows to complete missing values and, also, some non-missing values in order to improve the CR. To make this modification, a measure capable of controlling the magnitude of these changes was defined: Saaty’s Scale Distance (SSD) [29]. This magnitude is computed as the absolute difference between two values position in that scale:

SSDv1,v2 = | pos(v1 ) − pos(v1 ) |

(4.11)

where v1 and v2 belong to the Saaty’s scale. Note that the possible values for SSDv1,v2 are integers between 0 and 16, being 0 if the network output coincides with the original value. This mechanism permits to modify original values (missing and non-missing) having control over changes. Using the SSDv1,v2 concept, it is easy to create a matrix of distances which contains, in every position, the SSDv1,v2 between the original matrix value in that position and the estimated by the neural network. This matrix of distances provides an efficient control over the values to be modified; only the type of replacement must be indicated. That is, if the values will be replaced by the network output whose distances are lower, higher or equal to a specified value of SSDv1,v2 . In Fig. 4.1, the summarized process is shown. The top triangular portion of the normalized pairwise matrix is put on the trained network input and then, after the MLP recognition process, the original and the estimated top triangular portion matrices, the SSDv1,v2 value and the type of replacement are used in the reconstruction process in order to obtain the complete improved matrix. At last, the complete improved matrix is used to calculate the priority vector w .

114

Matrix order

Missing values

5 6 7 8 9 5 6 7 8 9 5 6 7 8 9

2

3

4

Original CR

CP CR

RGM CR

MLP CR

CP Ranking Coincid.

0.0986 0.0850 0.0779 0.0647 0.0591 0.0932 0.0864 0.0720 0.0614 0.0559 0.0992 0.0839 0.0768 0.0634 0.0570

0.0678 0.0664 0.0683 0.0596 0.0550 0.0573 0.0640 0.0593 0.0536 0.0515 0.0448 0.0551 0.0587 0.0521 0.0490

0.0668 0.0662 0.0686 0.0598 0.0551 0.0528 0.0629 0.0591 0.0536 0.0515 0.0519 0.0523 0.0581 0.0519 0.0490

0.0467 0.0442 0.0379 0.0399 0.0357 0.0395 0.0384 0.0378 0.0336 0.0364 0.0337 0.0339 0.0333 0.0324 0.0315

558 609 649 694 701 429 468 541 571 589 311 373 425 423 494

RGM Ranking Coincid.

MLP Ranking Coincid.

CP Top Coincid.

RGM Top Coincid.

MLP Top Coincid.

CP RDR2

RGM RDR2

MLP RDR2

558 573 638 679 682 420 452 548 565 574 305 391 413 450 496

506 403 256 201 151 394 311 256 178 95 315 246 205 127 98

804 875 894 916 933 762 807 858 880 905 661 761 806 848 873

807 863 887 915 931 756 807 862 888 906 664 773 802 851 871

780 747 700 628 632 720 709 676 629 583 635 648 648 569 571

15.52% 10.40% 08.02% 05.74% 04.76% 18.95% 15.43% 10.58% 09.08% 07.09% 26.18% 19.09% 14.12% 11.68% 09.15%

15.27% 11.27% 08.67% 06.02% 04.75% 19.75% 15.69% 10.56% 08.84% 07.09% 25.8%4 18.58% 14.54% 11.39% 09.24%

17.08% 18.57% 21.42% 23.37% 22.68% 21.26% 22.78% 21.80% 23.37% 24.45% 26.35% 25.47% 23.58% 26.39% 26.03%

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

Table 4.2 Comparison of missing values estimation techniques.

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

115

4.4. Measuring the reliability of missing values estimation techniques As it is explained above, to measure the reliability of reconstruction method, the original and enhanced priority vectors are used. For the simulations (as done in Section 4.1) one thousand of complete and consistent matrices for each matrix order (from 5 to 9) were created and the consistency ratio mean was calculated (original CR column in Table 4.2). The following process was performed in order to obtain three sets for each matrix order: 2 values were randomly removed from each matrix to obtain the first set, 3 values to obtain the second set and 4 values to obtain the third set (missing values column in Table 4.2). In this way, three sets with different level of incompleteness were used to test CP, RGM and MLP methods. Table 4.2 shows, for each matrix order and for each missing values quantity, the mean CR obtained by every method described in Section 4.3 (CP CR, RGM CR and MLP CR columns), the times that all rankings coincide after reconstruction (CP Ranking Coincid., RGM Ranking Coincid. and MLP Ranking Coincid. columns), the times that the top alternative coincides (CP Top Coincid., RGM Top Coincid., and MLP Top Coincid. columns) and, finally, the Relative Divergence Ratio using the multiplicative factor m = 2 (CP RDR2 , RGM RDR2 , and MLP RDR2 columns). The same conditions indicated in [29] for training MLP were used. That is, the same type of learning set (divided into training, validation and test sets), network topology, initial parameters and stopping conditions were used. Finally, the high value of SSDv1,v2 is used as type of replacement in order to change the original matrix values in every reconstructed matrix. The obtained results can be observed in Table 4.2. All tested methods reduce the original matrix inconsistency, but CP and RGM do not significantly improve the CR for orders 8 and 9. The highest consistency reduction is obtained with the MLP method. This occurs for all matrix orders and all missing values quantity. Under CR reduction criterion, the best method is MLP. However, as it happened with the tests in Section 4.1, the best method for improving matrix consistency generates the highest variation respect to the OPV. But this behavior occurs only for high order matrices (7, 8 and 9) while in low order matrices the ranking and top alternative coincidences are similar. An interesting point is that for CP and RGM, both of the coincidence quantities mentioned grow when the order matrix increases, and for MLP this quantity decreases for high order matrices. This situation is reflected in the mean CR evolution, i.e., when the CR improvement decreases the coincidence quantity increases. The OPV and EPV dissimilarity can be observed using RDR2 percentages. Clearly, the same behavior explained above is expressed by the RDR2 evolution. For orders 5 and 6, the difference between CP, RGM and MLP is not significant. In addition, when the RDR2 percentage decreases (because the ranking and top element coincidence number grows) the CR improvement decreases. This situation demonstrates that the Relative Divergence Ratio is an adequate parameter to measure the reconstruction methods reliability. Now, returning to the question: what is the best reconstruction method? Evidently, here MLP is more efficient and CP and RGM are more reliable, but this does not happen for all matrix orders. So, it is correct to say that MLP is the best method for low orders. In the same way CP and RGM are the best alternatives of reconstruction matrices of high orders. To answer the question above, another question comes up: will there be some way to find the balance between efficiency and reliability? In the next section, a Genetic Algorithm that considers both criteria is presented. 5. Bio-inspired models for pairwise matrix reconstruction In this section AHP matrices reconstruction is approached by using two types of bio-inspired models [47]: the Genetic Algorithms and the Swarm Intelligent Algorithms. They are both used to solve optimization problems and both of them use the concept of individuals belonging to populations. In this context the individuals’ characteristics are represented by data structures which contain the values of main variables used to solve the optimization problem. Also, every individual has a fitness value that indicates how good a solution such individual is. In order to find the best individual their values are modified by several operators or methods depending on the type of algorithm used. For AHP matrices reconstruction the same data structure and fitness function are used. Then, the performance of both approaches is tested in order to prove that the proposed reliability measure can be used in a wide range of reconstruction methods. 5.1. Genetic Algorithms generalities A Genetic Algorithm (GA) is an Artificial Intelligence technique used in large-space searches and optimization, inspired by natural selection and natural genetics [48]. As it occurs in any search task, the GA tries to efficiently find a solution to a problem in a space of candidate solutions [49]. As an optimization technique, the GA tries to approximate to the local maximal or minimal points of a target function. In general GA achieves good results, even when the function has several local maximal or minimal points. A GA simulates the behavior of natural species populations, evolving their individuals through several generations. In order to use a GA, the definition of two elements is required: the individual structure (or chromosome) and the fitness function. Each structure reflects a possible solution to the tackled problem. The fitness function guides the search process among iterations obtaining better solutions [50] and returns the value that designates how fit the individual’s structure to solve the problem is. After these two definitions, any combination of genetic operators can be used and tested until the GA provides satisfactory results.

116

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

START

Generate The Initial Population

Current Population

NOT MET

MET Stop Condition

Apply Selection Operator Apply Crossover Operator Apply Mutation Operaror

Return Solution

END

New Population

Fig. 5.1. GA process overview.

The fitness value is used to select the individuals for each operator: selection, crossover and mutation [49–51]. Also, this value can be used for stopping the GA process, especially when an individual reaches a required fitness level. Conceptual simplicity, broad applicability, hybridization with other methods, and parallelism are strong characteristics of GAs. They are robust to dynamical changes and they solve problems where the solution requires many resources if other techniques are applied. Additionally, they can solve problems which do not have a clear heuristic function to define the strength of a solution [52]. In Fig. 5.1 the GA process can be observed. The first step produces the initial population. This population is usually completely random, but some knowledge about the problem could be introduced. With each completed population, the GA checks if this population meets certain conditions to ending the process. If this condition is met, the best individual(s) is/are returned. Otherwise, this population passes through the three operators: selection, crossover and mutation. These operators produce new individuals, which are part of a new population, and the process starts again. Using the selection operator, a copy of the best individuals is put in the next generation. Crossover operator produces new individuals using two or more individuals of the current generation and mutation makes a copy of the structure changing some characteristic randomly. Selection allows to keep good "genetic material" across iterations while crossover and mutation bring diversity in [49,50]. 5.2. Swarm intelligence generalities: the Firefly Algorithm Swarm intelligence-based algorithms are a special class of bio-inspired algorithms. This kind of algorithms uses swarm behavior as strategy to obtain good solutions for many optimization problems. Some popular Swarm Algorithms are particle swarm optimization [53,54], ant colony optimization [55], cuckoo search [56,57], and Firefly Algorithm [58–60]. The last one is based on the fireflies’ light intensity and the attraction generated by their brightness. This type of algorithm is widely used, given its great performance and simplicity. The Firefly Algorithm (FA) was developed by Xin-She Yang and it has the following characteristics [59]: • All fireflies are unisex. • Their attractiveness is proportional to their brightness. • The landscape of the objective function influences the brightness of a firefly. In general, for an optimization problem, the brightness is proportional to the objective function (similar to the fitness function of GA). In the FA, there are two main issues: the light intensity variation and the attractiveness formulation. The

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

117

Algorithm 5.1 The Firefly Algorithm. Input: Objective function f (x ); x = [x1 , x2 , ..., xd ]; β0 ; γ ; α ; size population ni; maximum generation number MG; Output: Best vector solution x for f (x ); Process: 01. Initialize randomly the initial population of fireflies i at xi with i = 1, 2, ..., ni; t = 0; 02. Compute the light intensity Ii := f (xi ) for all fireflies using the objective function; 03. While (t < MG ) 04. For i := 1 To ni 05. For j := 1 To ni 06. If (Ii < I j ) 07. Compute ri2j and ε ; 08. xi := xi + β0 e−γ ri j (x j − xi ) + αεi ; 09. Update light intensity of firefly i, Ii := f (xi ); 10. EndIf 11. EndFor 12. EndFor 13. EndWhile 14. Return (best vector solution x); 2

light intensity is associated with the encoded objective function and the attractiveness of a firefly is determined by its light intensity (brightness). If the conditions above are fulfilled, the brightness is proportional to the objective function and the attractiveness depends on the firefly’s brightness and the scenery conditions (i.e. how other fireflies perceive that brightness). Intuitively, the attractiveness varies in function of the distance between two fireflies because part of the light is absorbed by the landscape. Yang proposes the form of the light intensity I as [59]:

I = I0 e−γ ri j 2

(5.1)

is the Cartesian distance between fireflies i and where I0 is the original light intensity, γ is the absorption coefficient and j. In general Ii = f (xi ) and since the attractiveness is proportional to the intensity, the attraction coefficient β is defined by using the absorption coefficient γ : ri2j

β = β0 e−γ ri j 2

(5.2)

where β0 is the attractiveness at ri j = 0. The behavior of the FA determines that the fireflies move to the brighter ones. This movement is computed by [59]:



xi = xi + β0 e−γ ri j x j − xi + αεi 2

(5.3)

where i and j are two fireflies at xi and x j positions respectively (firefly j is brighter than i), εi is a vector of random variables (drawn of Gaussian distribution) and α is a step-size parameter for εi . In Eq. (5.3) it can be observed that the movement to the brighter firefly is strongly influenced by the attractiveness (second term) and by a random variation (third term). The FA is summarized in Algorithm 5.1. 5.3. The GA-based and FA-based proposed models The reconstruction of AHP matrices can be viewed as an optimization problem and GA and FA can be used to solve it. Preliminary, both techniques are good methods to face the two reconstruction tasks simultaneously: (a) to complete missing values and (b) to improve matrix consistency. This simultaneity is equal to the MLP capacity explained in [29]. As explained before, the individual structure and the fitness/brightness function (fitness function for GA and brightness function for FA) must be defined taking optimization problem characteristics into account. Notice that the individual structure is the same for GA and FA. Also, the fitness function is equal to the brightness function. For these reasons, individual structure and evaluation function are explained for both of the methods which have been mentioned. As seen in Section 2, AHP pairwise matrices are reciprocal. This feature allows to use only the superior triangular portion of the matrix and to obtain the rest of its elements easily. Using this approach, redundancy is avoided and fewer computational resources are consumed. Therefore, the individual structure can be viewed as a set of values corresponding to the superior triangular portion of the matrix that contains missing and non-missing values. In this context, the individual structure represents the expert judgments (non-missing values) and the lack of preference information (missing values). The purpose of the GA-based and FA-based proposed models is to complete and make a n × n pairwise matrix A consistent. For this matrix an individual IA that represents such matrix and its elements are defined as an array of m = (n2 − n )/2 elements. Each array element is taken from the superior triangular portion of A. The individual’s structure is then an array defined by:

IA = [v1 , v2 , . . . , vi , ...vm ]

(5.4)

118

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

where:

  vi = a value of Saaty s Scale if there is a non-missing element in the ith position

(5.5)

0 if there is a missing element in the ith position

It can be observed that to obtain IA from A (as to obtain A from IA ) is a direct process. The main idea of the GA-based and FA-based proposed models is to generate populations of ni matrices Aj associated to A in order to fill incomplete positions and improve their CR (where ni is the number of individuals in the population and 1 ≤ j ≤ ni). The array representation of A can be extended to Aj as follows: 

IA j =



v 1 , v 2 , . . . , v i , . . . v m



(5.6)

where:

vi = a value of Saaty s Scale

(5.7)

Notice that no missing values exist for A j owing to the fact that the missing positions are filled with values of the Saaty’s Scale in order to obtain a complete matrix. Here arises one of several advantages of the GA-based and FA-based proposed models: flexibility. If only to complete the matrix A is required, the non-missing values must be fixed in A j . That is, vi = vi for all non-missing values in A. Also, if only the consistency improvement is required (there are no missing values in A) several combinations of values vi can be tested in order to obtain an adequate value of CR (these changes over non-missing values can be controlled by a distance measure as Saaty Scale Distance or SSDv ,v proposed in [29]). Now, the i

i

most important issue about this GA-based proposed model is that both requirements mentioned above can be addressed simultaneously. This is possible if the fitness/brightness function is designed taking both situations into account. The goal of the fitness/brightness function is to assign a score to the individual. This score represents how good the individual is (as problem solution). This function has to lead the GA and FA processes to return an individual that achieves two goals: (a) complete the matrix adequately and (b) improve matrix consistency. Also, this function must keep the reconstructed matrix A as similar as possible to the original matrix A, and consequently, the enhanced priority vector, EPV, as similar as possible to the original priority vector, OPV. Consequently, this fitness/brightness function must consider the following statements: 

• If the CR from the matrix A j decreases, then the fitness/brightness of IA

j

increases.

 • If the amount of change in A j with respect to A decreases, then the fitness/brightness of IA j increases.

In order to limit the values of individual fitness/brightness, the range [0, 10] is established (best individuals have high  fitness/brightness value). In this way the fitness/brightness(IA j ) function is defined as:



fitness/brightness(I

A j

)=





A j A A j 10  · [1 −CR(I )] · LH(I , I ) 

0.9

A j

CR(I

)

· LH(IA , IA j )

when



CR(IA j ) < 0.1

when 

CR(IA j ) ≥ 0.1

(5.8)

where:    CR(IA j ) is a function that reconstructs A j from IA j and returns the CR of this reconstructed matrix and LH(IA , IA j ) is a 

function which expresses the likelihood between IA and IA j . This function is defined by:





LH IA , IA j = 1 −







RDIST IA , IA j

(5.9)

where:  RDIST(IA , IA j ) and DIFF(vi , v i ) are defined by: 

RDIST(IA , IA j ) =



DIFF

vi

, v

i

 =

m

DIFF(vi , v i ) m

(5.10)

when vi =  v i when vi = v i

(5.11)

i=1

1 0 

Since the RDIST(IA , IA j ) returns a value in [0, 1] range, the square root increases it and, therefore, reduces the likelihood. This function is applied only to the known values of IA (i.e. vi = 0), and it divides the number of positions that have been  changed by the number of elements of the array (m). If more elements have been altered by the technique, RDIST(IA , IA j ) value increases. Finally DIFF(vi , v i ) is a binary function that expresses that there is a difference between vi and vi . In brief, this fitness/brightness function considers the CR improvement and the changes made into original matrix A to assign a fitness/brightness value to the reconstructed matrix A . When the CR and the changes tend to zero, better is the individual, achieving more consistency and reducing the difference between A and A . This produces a smaller difference between the Original and Enhanced Priorities Vectors.

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

119

Table 5.1 Simulation parameters. GA parameters Population (ni) Selection Crossover Mutation Maximum number of iterations (MG) Stop condition #1 Stop condition #2 FA parameters Population (ni)

β0 γ α

Maximum number of iterations (MG) Stop condition #1

200 individuals Elitist (10%) Roulette selection (adaptive percentage 80–60%) plus simple crossover Random selection (adaptive percentage 10–30%) plus simple mutation 150 Iteration limit reached Population convergence (98%) 50 individuals 1 1 − 0.5 × (iterationnumber/MG ) 1 − (iterationnumber/MG ) 100 Iteration limit reached

Table 5.2 GA-based and FA-based models results.

Matrix order

Missing values

5 6 7 8 9 5 6 7 8 9 5 6 7 8 9

2

3

4

Original CR

GA CR

GA Ranking Coincid.

0.0986 0.0850 0.0779 0.0647 0.0591 0.0932 0.0864 0.0720 0.0614 0.0559 0.0992 0.0839 0.0768 0.0634 0.0570

0.0697 0.0735 0.0709 0.0603 0.0560 0.0541 0.0670 0.0664 0.0579 0.0545 0.0422 0.0614 0.0634 0.0571 0.0536

718 729 671 671 679 588 592 584 559 585 476 523 478 494 490

GA Top Coincid.

GA RDR2

913 878 841 885 895 885 798 810 834 867 850 772 759 810 824

03.98% 04.03% 04.31% 02.92% 02.26% 05.81% 06.39% 05.35% 04.26% 02.92% 07.85% 07.25% 06.78% 04.94% 03.91%

FA CR

FA Ranking Coincid.

FA Top Coincid.

FA RDR2

0.0979 0.0938 0.0859 0.0711 0.0644 0.1224 0.1203 0.1053 0.0863 0.0765 0.1622 0.1536 0.1340 0.1093 0.0960

669 518 457 425 411 473 369 316 280 225 300 263 229 167 142

893 732 747 735 746 833 670 639 669 656 762 602 571 561 583

04.97% 08.26% 07.60% 06.77% 05.70% 08.67% 11.13% 10.69% 09.11% 08.66% 13.16% 14.15% 13.59% 12.55% 10.61%

5.4. Measuring the reliability of the GA-based and FA-based proposed models The simulation parameters used for the tests of each technique are shown in Table 5.1. GA-based model uses populations of 200 individuals in 150 iterations. Elitist selection, simple crossover and mutation are used as genetic operators. FA-based model uses populations of 50 individuals in 100 iterations and the values of parameters β0 , γ and α are shown in Table 5.1. The aforementioned values for both methods are resultant of the best balance between performance and computational cost (the GA was tested with populations of 150, 200 and 250 individuals in 100, 150 and 200 iterations and FA was tested with populations of 25, 50 and 100 individuals in 100, 150 and 200 iterations). In order to measure the reliability of GA-based and FA-based models, the original and enhanced priority vectors were used. In this case the OPV is obtained from original matrix A and the EPV is obtained from A (the best A j after simulation process). The same sets of complete and consistent matrices A of Section 4.2 were used. In this way, the consistency ratio means are the same in both sections (original CR column in Tables 4.2 and 5.2). Also, the three sets for each matrix order with 2, 3 and 4 random missing values were used (missing values column in Table 5.2). Table 5.2 shows, for each matrix order and for each missing values quantity, the mean CR obtained with the GA-based and FA-based models (GA CR and FA CR columns), the times that all rankings coincide after reconstruction (GA Ranking Coincid. and FA Ranking Coincid. columns), the times that the top alternative coincides (GA Top Coincid. and FA Top Coincid. column) and finally, Relative Divergence Ratio using the multiplicative factor m = 2 (GA RDR2 and FA RDR2 column). In Table 5.2 it can be observed that GA-based model reduces the original matrix inconsistency in all cases. On the other hand, this does not occur when the FA-based model is used. Even more, the FA does not reduce the original CR for any case. Regarding the ranking coincidence and top element coincidence GA and FA models decrease their performance when the missing values number grows. Beyond this, GA obtains better results than FA in all cases. All mentioned results can be observed in RDR2 percentages obtained by the two models. Clearly, the best performance of the GA model is reflected by the Relative Divergence Ratio. It is evident that there exists a better exploration of the search space by using GA model. This better exploration is based on the crossover and mutation operators because they generate new very dissimilar individuals throughout the search process. On the contrary, FA has a high dependence of best initial fireflies affecting the movement of

120

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

0.12

2 Missing Values

0.1

CR

0.08

Original CP

0.06

RGM MLP

0.04

GA 0.02

FA

0 4

5

6

0.14

7 8 Matrix Order

9

10

3 Missing Values

0.12 0.1 Original CP

0.06

RGM

CR

0.08

MLP 0.04

GA

0.02

FA

0 4

5

6

0.18

7 8 Matrix Order

9

10

4 Missing Values

0.16 0.14

CR

0.12

Original

0.1

CP

0.08

RGM

0.06

MLP

0.04

GA

0.02

FA

0 4

5

6

7 8 Matrix Order

9

10

Fig. 5.2. Comparison of CR improvement with different techniques.

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

30.00%

121

2 Missing Values

25.00% 20.00%

RDR2

CP 15.00%

RGM MLP

10.00%

GA FA

5.00% 0.00% 4

5

6

7

8

9

10

Matrix Order

30.00%

3 Missing Values

25.00% 20.00%

RDR2

CP RGM

15.00%

MLP 10.00%

GA FA

5.00% 0.00% 4

5

30.00%

6

7 8 Matrix Order

9

10

4 Missing Values

25.00% 20.00%

RDR2

CP RGM

15.00%

MLP 10.00%

GA FA

5.00% 0.00% 4

5

6

7 8 Matrix Order

9

10

Fig. 5.3. RDR2 comparison with different techniques.

122

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

the rest of them (i.e. they are very attracted by local maximums) and the random movements are insufficient to solve this problem. Even more, in Fig. 5.2 it can be observed that the FA does not reduce the CR adequately in comparison with the rest of the methods (CP, RGM and MLP). Beyond this, the Relative Divergence Ratio obtained by FA is better than MLP and closer to CP and RGM (Fig. 5.3). In Fig. 5.2 it can be observed that, using the GA-based model, the highest consistency reduction is obtained for a high number of missing values. Under CR reduction criterion, GA-based method has a similar performance to CP and RGM methods. However, it generates the lower variation respect to the OPV (see also Table 4.2). An interesting point is that the ranking coincidence is around 50% or higher for all order matrices (this coincidence decreases when the number of missing values grows). This behavior is logical because more CR reduction is observed in these cases. This also occurs with the top element coincidence, but the performance is more than 75%. Also, in Fig. 5.3 it can be observed that the GA-based proposed model is more reliable than CP, RGM and MLP, and its behavior can also be modified by changing the fitness function. In this way, some aspects of matrix reconstruction (CR improvement, non-missing values changes, priorities preservation, etc.) can be reinforced in the fitness function definition. As in previous sections the OPV and EPV dissimilarity can be observed using RDR2 percentages. Clearly, the top element coincidence influences the Divergence Ratio because the multiplicative factor (m = 2) assigns high importance to the first element of the ranking. These low values also show that, in most cases, only the less important alternatives in the rankings are altered. Again, this situation demonstrates that the Relative Divergence Ratio is an adequate parameter to measure the reconstruction methods reliability. 6. Conclusions In this paper a new efficient method to determine the confidence of matrix reconstruction methods has been proposed. Additionally, a GA-based reconstruction model for AHP matrices has been presented. The Relative Error Ratio, RDRm , defined in Section 3 can be efficiently used to indicate how different two rankings are. It is useful to compare priority vectors because it can handle the importance assigned to the top alternative and the total ranking coincidences. The multiplicative factor, m, provides a flexible way to reflect the importance of the two mentioned coincidence types over the comparison process. In simulations of Sections 4 and 5 the inversely proportional behavior of the CR and RDRm evolutions can be observed. This behavior demonstrates that the most efficient matrix reconstruction methods (i.e. those that achieve high inconsistency reduction) are less reliable. The proposed reliability measure can be extensively used by any ranking-based DSS. That is, the RDRm is not limited to AHP, because it can compare any pair of priority vectors. Furthermore, it can be used for comparing human experts’ ranking against the ranking returned by any DSS. The GA-based proposed model uses an easy-to-understand technique, which can be implemented with relatively low effort. The model is flexible, since any combination of operators and parameters can be used, and it can show a better performance with the use of parallelism. Even more, the model allows greater emphasis to CR reduction or ranking preservation (setting adequately the GA’s parameters). Simulations of Section 5 show that, in all cases, the GA-based proposed model reduces the CR keeping low levels of RDRm . This indicates that, in most cases, the expert’s judgments persist to the changes introduced by the proposed model. In addition, the GA-based proposed model (as multi-layer perceptron model) can also fill an incomplete pairwise matrix, being this a considerable optimization to lighten the data loading process while achieving a low CR. Finally, simulations show that under confidence conditions and this new way to compare rankings, the GA-based proposed model is more adequate than connecting paths, revised geometric mean techniques and even other bio-inspired algorithms (specifically the Swarm Intelligence Algorithm called Firefly Algorithm). Acknowledgements We would like to thank the referees for their valuable comments and suggestions. This work has been partially supported by the Project UTNRE-2106 of National Technological University (Argentina). References [1] D. Arnott, G. Pervan, Eight key issues for the decision support systems discipline, Decis. Support Syst. 44 (2008) 657–672, doi:10.1016/j.dss.2007.09.003. [2] D.J. Power, R. Sharda, Decision support systems, Springer Handbook of Automation, Springer Berlin Heidelberg, 2009, pp. 1539–1548, doi:10.1007/ 978-3-540-78831-7_87. [3] S. Sadagopan, Enterprise resource planning, Encyclopedia of Information Systems, Elsevier, 2003, pp. 169–184, doi:10.1016/B0-12-227240-4/00060-5. [4] C.W. Holsapple, M.P. Sena, ERP plans and decision-support benefits, Decis. Support Syst. 38 (2005) 575–590, doi:10.1016/j.dss.2003.07.001. [5] S. Gao, H. Wang, D. Xu, Y. Wang, An intelligent agent-assisted decision support system for family financial planning, Decis. Support Syst. 44 (2007) 60–78, doi:10.1016/j.dss.2007.03.001. [6] P. Gray, Group decision support systems, Decis. Support Syst. 3 (1987) 233–242, doi:10.1016/0167-9236(87)90178-3. [7] M.A. Nour, D. (Chi-Chung) Yen, Group decision support systems, Inform. Manage. 23 (1992) 55–64, doi:10.1016/0378-7206(92)90008-4. [8] G.P. Pervan, A review of research in Group Support Systems: leaders, approaches and directions, Decis. Support Syst. 23 (1998) 149–159, doi:10.1016/ S0167-9236(98)00041-4.

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

123

[9] R.W. Robbins, W.A. Wallace, Decision support for ethical problem solving: a multi-agent approach, Decis. Support Syst. 43 (2007) 1571–1587, doi:10. 1016/j.dss.2006.03.003. [10] P. Bharati, A. Chaudhury, An empirical investigation of decision-making satisfaction in web-based decision support systems, Decis. Support Syst. 37 (2004) 187–197, doi:10.1016/S0167-9236(03)00006-X. [11] H.K. Bhargava, D.J. Power, D. Sun, Progress in Web-based decision support technologies, Decis. Support Syst. 43 (2007) 1083–1095, doi:10.1016/j.dss. 2005.07.002. [12] D.J. Power, R. Sharda, Model-driven decision support systems: concepts and research directions, Decis. Support Syst. 43 (2007) 1044–1061, doi:10.1016/ j.dss.2005.05.030. [13] J. Figueira, V. Mousseau, B. Roy, Electre methods, Multiple Criteria Decision Analysis: State of the Art Surveys, Springer-Verlag, New York, 2005, pp. 133–153, doi:10.1007/0-387-23081-5_4. [14] J.-P. Brans, B. Mareschal, Promethee methods, Multiple Criteria Decision Analysis: State of the Art Surveys, Springer-Verlag, New York, 2005, pp. 163– 186, doi:10.1007/0-387-23081-5_5. [15] T.L. Saaty, The Analytic Hierarchy Process, McGraw-Hill, New York, 1980. [16] T.L. Saaty, Fundamentals of the analytic network process – Dependence and feedback in decision-making with a single network, J. Syst. Sci. Syst. Eng. 13 (2004) 129–157, doi:10.1007/s11518-006-0158-y. [17] T.O. Boucher, E.L. MacStravic, Multiattribute evaluation within a present value framework and its relation to the analytic hierarchy process, Eng. Econ. 37 (1991) 1–32, doi:10.1080/00137919108903055. [18] K. Mandic, B. Delibasic, S. Knezevic, S. Benkovic, Analysis of the financial parameters of Serbian banks through the application of the fuzzy AHP and TOPSIS methods, Econ. Modell. 43 (2014) 30–37, doi:10.1016/j.econmod.2014.07.036. [19] J. Yang, P. Shi, Applying analytic hierarchy process in firm’s overall performance evaluation: a case study in China, Int. J. Bus. 7 (2002) 31–46. [20] R. Tyagi, C. Das, A methodology for cost versus service trade-offs in wholesale location-distribution using mathematical programming and analytic hierarchy process, J. Bus. Logist. 18 (1997) 77–99. [21] S.H. Ghodsypour, C. O’Brien, A decision support system for supplier selection using an integrated analytic hierarchy process and linear programming, Int. J. Prod. Econ. 56–57 (1998) 199–212, doi:10.1016/S0925-5273(97)00009-1. [22] J. Korpela, K. Kyläheiko, A. Lehmusvaara, M. Tuominen, An analytic approach to production capacity allocation and supply chain design, Int. J. Prod. Econ. 78 (2002) 187–195, doi:10.1016/S0925-5273(01)00101-3. [23] G. Wang, S.H. Huang, J.P. Dismukes, Manufacturing supply chain design and evaluation, Int. J. Adv. Manuf. Technol. 25 (2005) 93–100, doi:10.1007/ s00170-003-1791-y. [24] X. Deng, Y. Hu, Y. Deng, S. Mahadevan, Supplier selection using AHP methodology extended by D numbers, Expert Syst. Appl. 41 (2014) 156–167, doi:10.1016/j.eswa.2013.07.018. [25] F. Lolli, A. Ishizaka, R. Gamberini, New AHP-based approaches for multi-criteria inventory classification, Int. J. Prod. Econ. 156 (2014) 62–74, doi:10. 1016/j.ijpe.2014.05.015. [26] S. Malladi, K.J. Min, Decision support models for the selection of internet access technologies in rural communities, Telematics Inform. 22 (2005) 201–219, doi:10.1016/j.tele.2004.10.001. [27] G.S. Kearns, A multi-objective, multi-criteria approach for evaluating IT investments: results from two case studies, Inform. Resour. Manage. J. 17 (2004) 37–62. [28] N. Kwak, C.W. Lee, Business process reengineering for health-care system using multicriteria mathematical programming, Eur. J. Oper. Res. 140 (2002) 447–458, doi:10.1016/S0377-2217(02)00082-6. [29] J.A. Gomez-Ruiz, M. Karanik, J.I. Peláez, Estimation of missing judgments in AHP pairwise matrices using a neural network-based model, Appl. Math. Comput. 216 (2010) 2959–2975, doi:10.1016/j.amc.2010.04.009. [30] Q. Chen, E. Triantaphyllou, Estimating data for multicriteria decision making problems: optimization techniques, Encyclopedia of Optimization, Springer US, Boston, MA, 2001, pp. 567–576, doi:10.1007/0-306-48332-7_123. [31] P.T. Harker, Alternative modes of questioning in the analytic hierarchy process, Math. Modell. 9 (1987) 353–360, doi:10.1016/0270-0255(87)90492-1. [32] S. Shiraishi, T. Obata, M. Daigo, Properties of a positive reciprocal matrix and their application to AHP, J. Oper. Res. Soc. Jpn. 41 (1998) 413–414. [33] M.T. Lamata, J.I. Pelaez, A method for improving the consistency of judgments, Int. J. Uncertain. Fuzziness Knowl.-Based Syst. 10 (2002) 677–686, doi:10.1142/S0218488502001727. [34] T.L. Saaty, Decision-making with the AHP: why is the principal eigenvector necessary, Eur. J. Oper. Res. 145 (2003) 85–91, doi:10.1016/S0377-2217(02) 00227-8. [35] C.A. Bana e Costa, J.-C. Vansnick, A critical analysis of the eigenvalue method used to derive priorities in AHP, Eur. J. Oper. Res. 187 (2008) 1422–1428, doi:10.1016/j.ejor.2006.09.022. [36] K. Kułakowski, Notes on order preservation and consistency in AHP, Eur. J. Oper. Res. 245 (2015) 333–337, doi:10.1016/j.ejor.2015.03.010. [37] T.L. Saaty, How to make a decision: the analytic hierarchy process, Eur. J. Oper. Res. 48 (1990) 9–26, doi:10.1016/0377-2217(90)90057-I. [38] T.L. Saaty, Axiomatic foundation of the analytic hierarchy process, Manage. Sci. 32 (1986) 841–855, doi:10.1287/mnsc.32.7.841. [39] R.L. Burden, J.D. Faires, Numerical analysis, 8th ed., Brookes Cole, 2004. [40] S.C. Chapra, R.P. Canale, Numerical Methods for Engineers, McGraw-Hill, 2012. [41] Y.-C. Hu, J.-F. Tsai, Backpropagation multi-layer perceptron for incomplete pairwise comparison matrices in analytic hierarchy process, Appl. Math. Comput. 180 (2006) 53–62, doi:10.1016/j.amc.2005.11.132. [42] J.I. Peláez, M.T. Lamata, A new measure of consistency for positive reciprocal matrices, Comput. Math. Appl. 46 (2003) 1839–1845, doi:10.1016/ S0898-1221(03)90240-9. [43] J. Aczél, T.L. Saaty, Procedures for synthesizing ratio judgements, J. Math. Psychol. 27 (1983) 93–102, doi:10.1016/0022-2496(83)90028-7. [44] P.T. Harker, Incomplete pairwise comparisons in the analytic hierarchy process, Math. Modell. 9 (1987) 837–848, doi:10.1016/0270-0255(87)90503-3. [45] A. Stam, M. Sun, M. Haines, Artificial neural network representations for hierarchical preference structures, Comput. Oper. Res. 23 (1996) 1191–1201, doi:10.1016/S0305-0548(96)00021-4. [46] S. Haykin, Neural Networks A Comprehensive Introduction, Prentice Hall, New Jersey, 1999. [47] I.J. Fister, X.-S. Yang, I. Fister, J. Brest, D. Fister, A brief review of nature-inspired algorithms for optimization, Neural Evol. Comput. 80 (2013) 1–7 http://arxiv.org/abs/1307.4186 . [48] D. Simon, Evolutionary Optimization Algorithms, Wiley, 2013. [49] M. Mitchell, An introduction to genetic algorithms, Comput. Math. Appl. 32 (1996) 133, doi:10.1016/S0898-1221(96)90227-8. [50] D.E. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning, first ed., Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA, 1989. [51] R. Haupt, S. Haupt, The continuous genetic algorithm, Practical Genetic Algorithms, John Wiley & Sons, Inc., Hoboken, NJ, USA, 2004, pp. 51–66, doi:10.1002/0471671746.ch3. [52] S.N. Sivanandam, S.N. Deepa, Introduction to Genetic Algorithms, first ed., Springer Publishing Company, Incorporated, 2007. [53] J. Kennedy, R. Eberhart, Particle swarm optimization, in: Proceedings of ICNN’95 – International Conference on Neural Networks, IEEE, 1995, pp. 1942– 1948, doi:10.1109/ICNN.1995.488968. [54] R. Poli, J. Kennedy, T. Blackwell, Particle swarm optimization, Swarm Intell. 1 (2007) 33–57, doi:10.1007/s11721-007-0002-0. [55] M. Dorigo, M. Birattari, T. Stutzle, Ant colony optimization, IEEE Comput. Intell. Mag. 1 (2006) 28–39, doi:10.1109/MCI.2006.329691. [56] X.S. Yang, S. Deb, Cuckoo search via Lévy flights, in: 2009 World Congress on Nature and Biologically Inspired Computing, NABIC 2009 – Proceedings, 2009, pp. 210–214, doi:10.1109/NABIC.2009.5393690.

124

M. Karanik et al. / Applied Mathematics and Computation 279 (2016) 103–124

[57] X.S. Yang, S. Deb, Engineering optimisation by cuckoo search, Math. Modell. Numer. Optim. (2010) 1–17, doi:10.1504/IJMMNO.2010.035430. [58] X.S. Yang, Firefly algorithms for multimodal optimization, Lecture Notes in Computer Science (including Subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), Springer Berlin Heidelberg, 2009, pp. 169–178, doi:10.1007/978-3-642-04944-6_14. [59] X.S. Yang, Firefly algorithm, stochastic test functions and design optimization, Int. J. Bio-Inspired Comput. 2 (2) (2010) 78–84, doi:10.1504/IJBIC.2010. 032124. [60] I. Fister, I. Fister, X.-S. Yang, J. Brest, A comprehensive review of firefly algorithms, Swarm Evol. Comput. 13 (2013) 34–46, doi:10.1016/j.swevo.2013.06. 001.