A MAGDM Method Based on Possibility Distribution Hesitant Fuzzy Linguistic Term Set and Its Application

: The sustainable third-party reverse logistics provider (3PRLP) selection, as the core of sustainable supply chain management, has become paramount in research nowadays. In the actual evaluation process, the decision makers may hesitate in a few linguistic terms and have di ﬀ erent partiality towards each term, hence the possibility distribution based hesitant fuzzy linguistic term sets (PDHFLTSs), as expressed by a consecutive or non-consecutive linguistic term set, is suitable for such an evaluation. The purpose of this paper is to solve sustainable 3PRLP selection problems with linguistic information by developing an e ﬀ ective and robust method. We ﬁrstly redeﬁne the covariance-based correlation coe ﬃ cient that can simplify the computation to calculate the consensus degree, and then introduce the hesitant degree in context of possibility distribution information, in order to enrich measures of PDHFLTSs. On this basis, the weights of experts are computed for expression aggregation. Secondly, to overcome attributes’ weights staying constant, the combination of group utility function and variable weight approach is introduced to give the weights of attributes. Most importantly, a decision method, called MULTIMOORA, is optimized by improving the ranking position method, and then, through the combination with PDHFLTS, we proposed a possibility distribution based hesitant fuzzy linguistic MULTIMOORA method with great robustness. At last, the presented method is applied to the ﬁeld of sustainable third-party reverse logistics provider selection in the e-commerce express industry and the e ﬀ ectiveness is veriﬁed by several comparative analyses. Method.


Introduction
Sustainable third-party reverse logistics provider selection, as a typical multi-attribute group decision making (MAGDM) issue, is the procedure of choosing the best satisfaction provider and is appropriately accomplished by a group, which is composed of experts with different knowledge and interests background. In addition, the MAGDM problems mainly include three processes, which are (1) experts make evaluation recommendations of every alternatives under different attributes (or criteria); (2) individuals reach consensus and opinion aggregation; and, (3) alternatives ranking and selecting. How to express opinions of experts accurately has become a paramount in past decades, since the first step of this problem is the processing of expert's evaluative expressions. In practice, some researchers have found that decision makers (DMs) prefer describing the performance of alternatives by linguistic variables. Hence, Zadeh [1] developed the fuzzy linguistic method, which explains and uses a single involved in MAGDM problems, the cost of modifying the decision results is much higher than that of single DM, it is also very difficult to identify or select the optimal scheme, which is mainly due to the different personal preferences and existence of the noninferior solution of the scheme. Therefore, it is important to choose a decision-making method with robustness in order to better cater to these situations. The MULTIMOORA method [18], further generated by extending the full multiplicative form into MOORA (The Multi-Objective Optimization by Ratio Analysis), is a relatively new method with robust solution. It refers to utilizing three ranking models with different effects to obtain three subordinate orders, and then synthesizing each sub-result to the final order by dominance theory. In [19], MOORA was proven to have explicit advantages over TOPSIS [20], VIKOR, Analytic Hierarchy Process (AHP), ELECTRE, PROMETHEE as easy mathematics, low computational time, stability, and straightforwardness for DMs. Not only that, the VIKOR and TOPSIS methods derive complete ranks that are only based on utility values; the ELECTRE method requires setting parameters, which are not always of clear economic significance, but it also needs to pass the concordance test and non-discordance test, and its process is quite complex; the PROMETHEE method does not contain strong robustness and objectivity, in which the preference function required being subjectively selected by DMs, and the same as the TOPSIS method, the increase or decrease of alternatives will affect final results. Improved MULTIMOORA method, which is simple in mathematics, does not need any subjective setting, but it also takes into account both the utility values and the sorting results. Accordingly, it can have strong robustness. Nevertheless, (3) the primitive MULTIMOORA and other extension forms still have some defects: 1. Dominance theory only involves three subordinate ranks but excludes their utility values. 2. Its maneuverability decreased under a large number of alternatives environment due to the times of pairwise comparisons skyrocketed significantly. 3. The optimal scheme might be easily derived, but the ranks of all schemes are difficult to obtain. As a result, we will introduce an improved aggregation operator to correct defects (3).
Moreover, the weight of experts and the weight of attributes are another important part of MAGDM problem. In most of the existing literature, decision makers are assigned completely the same weight or subjectively given in advance under the group decision-making (GDM) topic. However, in real-time application, the importance of experts' expression must be different up to their knowledge and experience. Wu [14] investigate the consensus degree from the perspective of distance measure to determine expert' weights, which is widely recognized and used, but, to our best knowledge, (4) few paper calculates the consensus degree from the perspective of correlation coefficient, and the amount of calculation using these formulas is relatively large. On the other hand, Liao [21] develops the concept of hesitant degree to compute expert' weights. These scholars apply the consensus degree or hesitant degree [22] separately, to show the different status information of DMs in the decision making process. Unfortunately, (5) nobody utilized both at same time, and this is meaningful, in the presented paper, consensus and hesitant degree are together decided for decision makers' weight, that is to say, the higher hesitant degree or the lower consensus degree, the smaller the weight of experts will be allocated. Secondly, in the vast majority of MAGDM methods with hesitant fuzzy linguistic information, (6) neither the most classic one nor novel one can achieve the effect, where the same alternative has different weights under the same attribute, also say as those methods [23][24][25][26] only produce a constant weight vector of all attributes in spite of complexity of realistic MAGDM problem. In addition, (7) decision makers are always influenced by their preference-induced behavioral characteristics when they make an evaluation, but this aspect is not taken into account by previous GDM researches. Hence, it is necessary to solve problems (6) and (7) with a suitable approach.
To sum up, this paper aims to accomplish the following creative contributions: (I) We presented the correlation coefficient formula of PDHFLTS by combining an improved rule that can reduce the number of calculations, and further defined consensus degree to fill gap (4). (II) Considering the consensus degree [27] can be capable of reflecting certainty and consistency from a group of experts, by contrary, hesitant degree [28] delivers the uncertainty and hesitant information of DMs. To some extent, those concepts, in our cognition, show DMs' two kinds of opposite psychological characteristics.
Thus, we proposed a novel function for determining the expert's weight that will overcome the gap (5). (III) By individuals' utility aggregation to gain group utility, we will develop the variable weight method for solving problems (6) and (7) (IV) We improved the ranking position approach, which makes up for the shortage (3), and combined it with MULTIMOORA, which aims at problems (1) and (2). (V) A case study regarding selecting sustainable third-party reverse logistics providers for e-commerce express industry is analyzed by the presented method and the effectiveness is proven by a comparative study.
The contents and structure of current paper are organized, as follows: Section 2 briefly reviews the basic concepts of PDHLFLTS, the main ideas of variable weight method, and primitive MULTIMOORA. Section 3 presents the correlation-coefficient based consensus degree and hesitant degree based on PDHFLTS to determining jointly weights of expert. Section 4, by integrating individual utility into group utility, gives the variable weight decision method and shows some properties. Section 5 proposes a PDHFL-MULTIMOORA approach with possibility distribution based hesitant fuzzy linguistic environment. Section 6 illustrates that the presented method is a suitable and effective way to solve sustainable third-party reverse logistics provider (3PRLP) selection problem.

Preliminaries
To make the presentation of subsequent research content clear, this section will illustrate core knowledge regarding the PDHFLTS, the variable weight, and the MULTIMOORA method.

Possibility Distribution Based HFLTS
The conventional linguistic presentation tool only applies a single term to assist DMs in evaluation. After that, Wu [14] presents a novel linguistic method, named possibility distribution based hesitant fuzzy linguistic term set (PDHFLTS), which conveys DMs' expression further precisely and conveniently. Its definition is as follows.
Definition 5 ([14]). Suppose S, h S , P are the same as defined above, then the variance value function is defined as follows: where NS s ψ denotes score function of s ψ .

The Traditional MULTIMOORA Method
MULTIMOORA is a useful and novel MADM method. Its final result originated from the aggregation of ternary subordinate models: Ratio System, Reference Point, and Full Multiplicative Form, and generally aggregating them through dominance theory.
A general decision matrix can be defined as D = x ij m×n , where x ij denotes the performance value of alternative a i under criteria c j . Thus, the procedures of the traditional MULTIMOORA method can be defined, as follows [18]: Step 1. Normalize all x ij of decision matrix to x N ij by formula: Generally, m ≥ 2.
Step 2. Calculate the utility values of Ratio System (RS) by arithmetical weighted aggregation operator of RS model defined, as follows: x N ij (6) where g denotes the cardinality of beneficial attributes, in turn, the cardinality of cost attributes is (n − g). The top-ranking position corresponds to maximum utility value y i . Afterwards, the subordinate result in descending order is denoted as: Step 3. Calculate the utility values of Reference Point (RP) model. The main purpose of this approach is to find the worst performance of alternatives with respect to every attribute, and then choosing the best result from those worst rating. Hence, in this step, the final chosen alternative is the best of the worst ones.
where g denotes the cardinality of beneficial attributes, in turn, the cardinality of cost attributes is (n − g). The ranking of final results in ascending order is denoted, as follows: Step 4. Calculate the utility values of the Full Multiplicative Form (FMF) model, which is the basis of weighted geometric integration operator, as follows: where g denotes the cardinality of beneficial attributes, in turn, the cardinality of cost attributes is (n − g). The ranking of final results in descending order denotes, as follows: Step 5. Derive final ranking of all alternatives by using the dominance theory to aggregate three subordinate ranks into one.
In the aspect of application, the MULTIMOORA approach has been developed in many fields, including industry, economics, healthcare management, information and communication technology, environmental policy making, and so on. Specifically speaking, industry receives the most attention, which can be summarized as construction [32], automotive [47], mining [48], logistics [49], and manufacturing system [36]. Economics includes economic evaluation [50] and sustainable development [51]. Healthcare management is divided into biomedical service [33] and medical service [52]. ICT involves information technology [53] and telecommunication technology [54]. Environmental policy making contains climate change policy making [55] and supplier selection [36].

The Core of Variable Weight Theory
Since Wang [56] proposed the concept of variable weight and Li [57] and Yu [58] further discussed some important properties of this method, the variable weight approach has received a lot of attention due to its flexibility, which mainly manifested in changing with schemes' configuration.
Definition 8 ([57,59]). A n-dimensional variable weight vector is a n-ary mapping. ω j : , where x j is the j-th criteria of alternative a i and ω j is corresponding weight. Besides, it should satisfy four properties: normality, nonnegativity, monotonicity, and continuity.
Definition 9 ([57,59]). Suppose the vector S X is a mapping: . . , s n (x n ) . If S X satisfies following conditions: (1) s j (x 1 , . . . , x n ) is continuous with respect to x j .
(2) S X σ ij x j = S X x j , where σ ij is the function to exchange the position of the i-th element and the j-th element. . , x n ), . . . , w n (x 1 , . . . , x n )) is the Hadamard product of w c and S X . (4) s j (x 1 , . . . , x n ) is monotonous regarding every variable x j .
After that, S X is the variable weight state vector of alternative a i .
To our best knowledge, few scholars [59][60][61] have studied this approach in MAGDM or MADM background. However, as can be seen from the above definition, the greatest advantage of the variable weight method is that the same attribute has different weights under different scheme configuration, which is not possible for other weighting methods [34][35][36][37][38][39].

Nomenclature Section
For the more convenience of readers, Table 1 shows all acronyms.

Determination of Expert Weight Based on Consensus and Hesitant Degree
In a realistic situation, a wise decision was usually made by a group of individuals with professional opinions instead of a single person. Differentiated decision power in group aggregation process among experts should be considered, which is defined as differentiated potential capabilities of an individual or a group to exert own influence or control over other individual or group, since it is the cooperative problem of GDM background. Accordingly, it can be quantified as their weights based on different position, expertise, percentages of shares, etc. [62]. However, it is precisely because different levels and fields of each expert's expertise, resulting in every expert conveying different information, which creates conflicts. Hence, in the vast majority of MAGDM problems under hesitant fuzzy linguistic environment, the weight of expert is assumed to be the same, or sets subjectively in advance is unreasonable.
In this section, we develop a method for determining experts' weights from the perspective of group opinion aggregation and personal information expression. Furthermore, the consensus degree denotes the level of agreement between every decision-maker in a group, while the hesitant degree denotes the uncertainty level of information evaluation in each expert. We link those two concepts to detect the level to which an expert is close to other experts and uncertainty in his/her heart. Thus, the main idea of this method is to utilize consensus degree and hesitant degree of each decision-maker, that is, the higher hesitant degree or the lower consensus degree, the lower the weight of expert will be.

The Consensus Degree of the PDHFLTS
In recent consensus degree computation research, most of the main approaches are based on preference relation [63][64][65] and distance measure [14,15]. These papers rarely involve methods that are based on correlation coefficient, and, even if there are, the calculation is complicated. As shown in [3], as an useful tool for describing expert information and for studying a limited or small number of samples, HFLTS is more suitable for defining covariance from the perspective of information energy within unit interval [0, 1]. Besides, in most cases, the cardinality of elements in two PDHFLTS is always different, which is the main reason for complex computation. To our best knowledge, there is no literature on extending two PDHFLTS with different lengths to the same length of PDHFLTS. Therefore, it is necessary to define the PDHFLTS extension rule first for filling the research gap.

Definition 10.
Let , h a S 2 satisfy following rule: Step 1. Compare s 1 and s 2a Step 2. Continue to compare s 1 It is important to note that p 1i = p 2i + p s 1i+1 + p ss 1i+2 + . . . + p s...s 1i+k and p 2i = p 1i + p s 2i+1 + p ss 2i+2 + . . . + p s...s 2i+k , the possibilities corresponding to each linguistic variable remains the same. In addition, when the possibilities of two comparative elements are the same, it can be handled by any two conditions. The possibility permutation set PP * is different from P of PDHFLTS. Figure 1 shows the detailed calculation procedure. The covariance and correlation coefficient of possibility distribution are thus defined, as follows: However, without Definition 10, we are likely to calculate the covariance by following formula: We can see from above two covariance measures that      The covariance and correlation coefficient of possibility distribution are thus defined, as follows: Definition 11. Suppose S, H S , P are the same as defined before. Adjusted h 1a indicated by PP * is defined as: where NS s ψ denotes score function of s ψ , and s 1a be the same as example 1, suppose NS s ψ = ψ, we can obtain E h 1a (1). Subsequently, the covariance can be calculated by Equation (13) and the result is as follows: However, without Definition 10, we are likely to calculate the covariance by following formula: Additionally, the calculation process and result are: We can see from above two covariance measures that Cov h 1a which demonstrate that the computation is simplified by Definition 10, and the correctness of the results is guaranteed.
On basis of Definition 10, it is obvious for two PDHFLTSs h S 1 , h S 2 to satisfy those good properties: By defining the covariance of two PDHFLTs, we further develop the correlation coefficient formula.

Definition 11.
Let h S 1 , h S 2 be two PDHFLTs, then the correlation coefficient between h S 1 and h S 2 will be derived by following equation.
Based on the correlation coefficient measure, the consensus degree of the experts can be defined. , where d ij ⊂ PDHFLTSs denotes the performance of alternative X i under attribute C j . Afterwards, the consensus degree of alternative X i under attribute C j for the T-th decision maker can be defined as: The consensus degree of alternative X i for the T-th decision maker is defined as follows: Subsequently, the consensus degree for the T-th decision maker is defined: For the convenience of expression, we will adopt NS s ψ = ψ in following contents.

The Hesitant Degree of the PDHFLTS
In this section, the basic task is working on decision makers' information when they are hesitant in conveying their preferences. Since Liao et al. [28] first proposed the concept of hesitant degree, many scholars have studied it, such as [66,67]. However, the above hesitant degree is only related to the number of HFLTS or EHFLTS elements, for our cognition, that is because DMs tend to assign equal possibilities to every linguistic variable under hesitant situation, at beginning of evaluation, it is difficult for DMs to differentiate the possibilities of linguistic variables in a state of self-confidence.
Therefore, when hesitant fuzzy linguistic information added the possibility distribution, the hesitant degree is also related to the difference in the possibilities of linguistic variables. In this sense, we further define hesitant degree with possibility distribution. Definition 13. Let S = s δ δ = 1, 2, . . . , g be a LTS, X, C, E are the same as defined before and d ij is a PDHFLTS on S that denotes the performance of alternative X i under attribute C j . P ij is the corresponding possibility distribution of d ij and L ij is the cardinality of d ij , p a , p b ⊂ P ij , the hesitant degree of alternative X i under attribute C j for the T-th decision maker is thus given by: The consensus degree of alternative X i for T-th decision maker is defined as: The consensus degree for T-th decision maker is denoted, as follows: There is only single linguistic term, the consensus degree is 0. in turn, there are g linguistic variables and their possibilities remain equal, then the consensus degree is 1, that is 0 ≤ HD ≤ 1.
Based on consensus and hesitant degree, the weight of experts can be defined as: Suppose the consensus degree is CD T and the hesitant degree is HD T for the T-th expert, there are K(K ≥ 2) experts E = (E 1 , . . . , E K ) participated. Hence, we obtain: where ω T denotes the weight of expert T and 0 ≤ ω F ≤ 1,

Determination of Attribute Weight Based on the Variable Weight Theory
Each expert has various knowledge background and problem recognition ability regarding expressing their evaluation on different attributes; these factors work together to influence the behavior of evaluation of decision makers. In our cognition, the utility function depicts and quantifies the psychological cognition of DMs on the evaluation results. Thus, in this section, we will introduce the utility function into the variable weight method for determining the weight of the attributes.
In management decision problems, the weight of an attribute keeps unchanged with respect to different alternatives, which leads to unreasonable decision-making results. For instance, the safety risk assessment of amusement facilities is based on the multi-factor comprehensive results. If an abnormality occurs in one link, that is, the attribute value is abnormal, and other attribute values are maintained at normal levels, there is still a high probability of a safety accident. However, the existing papers [14,15,[23][24][25][26] cannot illustrate and cope with this situation. A more intuitive assumption is that there are attributes a and b, and their weight vector is (0.5, 0.5), so the score synthesis function is S = 0.5x a + 0.5x b , where x a , x b denote the performance of alternatives under the attribute a or b, respectively. If alternative 1 has x a = 0.05, x b = 0.95 while alternative 2 has x a = 0.5, x b = 0.5, according to common sense in this case, alternative 2 is far better than alternative 1. However, the fact is that their scores both are 0.5, which is unreasonable. Unless the attribute weight changes with the change of the scheme, it can overcome this defect. Hence, we should consider not only the preference for the relative importance of attributes, but also the configuration equilibrium of the attributes in each scheme.
In [59], Yu proposed a compromise-typed variable weight method to solve Hybrid MADM problems. However, it is also necessary to discuss a single-typed variable weight, and the utility-based variable weight has not been introduced into the GDM context. If combing the variable weight method with group utility, which is aggregating from individual utility, we can solve the decision problem with the GDM environment. Besides, if we transform all of the cost attributes into benefit attributes, then the single-typed variable weight function can deal with most of the MAGDM problems. Therefore, in this section, we necessarily proposed a novel incentive-typed variable weight method for determining the weight of attribute under the MAGDM environment.
∂x j 2 = 0, where x j denotes j-th attribute of alternative a i and0 ≤ x j ≤ 1, j = 1, . . . , n, then W(X) is called as the neutral-incentive type variable weight vector (NIVW).
, where x j denotes j-th attribute of alternative a i and 0 ≤ x j ≤ 1, j = 1, . . . , n,r j ⊂ [0, 1] is a reference point, then W(X) is called as the S-shaped-incentive type variable weight vector (SIVW).
Proof. Since U(X) is the S-shaped utility function, then U (X) > 0 and U " (X) < 0, X ⊂ [0, r],U " (X) > 0, X ⊂ [r, 1], thus we have: where x j is the j-th attribute of alternative a i , W j (X) = ω j x j denotes the j-th element of W(X), w j is the j-th element of w c .

Theorem 2. If U(X) is the inverse S-shaped utility function and
is the ISIVW.
is the CIVW.

Theorem 5. If U(X) is the neutral-typed utility function and
The proof process of Theorems 2-5 is similar to the Theorem 1.

Theorem 6.
If W(X) is the SIVW, then W −1 (X) is the ISIVW, where W −1 (X) is the inverse function of W(X).

Proof. Since W(X) is the SIVW, then
Besides, according to the properties of inverse function, we can derive symmetrically about y = x, and The proof process of Theorems 7-8 is similar to the Theorem 6. The rest of the work is to obtain the group utility.

The MULTIMOORA Method Based on PDHFLTS
In this section, we proposed the PDHFL-MULTIMOORA method. The weight of expert is determined by the approach that is based on the consensus and hesitant degree, while the weight of attribute is defined by the variable weight method in which the group utility is taken into consideration. The MULTIMOORA method is adjusted by introducing the improved ranking position method under PDHFL environment. It is necessary to give a brief introduction to MAGDM problems before we get down to the core.
A multi-attribute group decision-making problem consists of three parts: expert set, alternative set, and attribute set. The expert set is expressed as E = (E 1 , . . . , E K ), K ≥ 2, alternative set can be defined as A = (A 1 , . . . , A m ), m ≥ 2, and attribute set is denoted as C = (C 1 , . . . , C n ), n ≥ 2. The assessment of expert K with respect to alternative A i under attribute C j is expressed as

the PDHFL-MULTIMOORA Method with the Improved Ranking Position Method
We aim to improve the traditional MULTIMOORA approach, in this part, by introducing the ranking position method to define the novel MULTIMOORA method. First of all, owing to the advantages of the vector normalization, we use the vector normalization formula to normalize the group decision matrix D m×n : Further, we derive the normalized group decision matrix D N = . After this foundational work, the calculation of three models can be defined. (

1) The Possibility Distribution Hesitant Fuzzy Linguistic Ratio System (PDHFLRS)
To calculate the utility value of Ratio System, the arithmetic weighted operator computes the normalized value d N ij . Subsequently, the subordinate utility value of alternative A i is given by: where ω j denotes the weight of attribute C j ( j = 1, . . . , n). C j ( j = 1, . . . , g) belong to the benefit attribute, while C j ( j = g + 1, . . . , n) belong to cost attribute. The ranking results of each alternative are arranged in descending order by utility value U 1 (A i ), so the first subordinate orders are shown as: (2) The Possibility Distribution Hesitant Fuzzy Linguistic Reference Point (PDHFLRP) The subordinate utility value of alternative A i is given by: where r j is calculated by r j = maxd N ij , j ≤ g; mind N ij , j ≥ g . Ranking alternative by utility value U 2 (A i ) in ascending order, and then the second subordinate orders are shown as: (

3) The Possibility Distribution Hesitant Fuzzy Linguistic Full Multiplicative Form (PDHFLFMF)
According to geometric weighted operator, the subordinate utility value of alternative A i is calculated by: Alternatives are ranked by utility value U 3 (A i ) in descending order, and then the third subordinate ranking is denoted as:

(4) The final ranking obtained by improved Ranking Position method
In final stage, the rest of the work is the aggregation of three models. If we regard those models as attributes that represent different aspects of information in decision-making problems, then we can transform it into a MADM problem, in which every alternative A i has a pairwise values (utility value U k (A i ), ranking result R k (A i )) with respect to three attributes C k (k = 1, 2, 3). Therefore, there are two decision matrices that consist of utility value matrix D(U) and ranking matrix D(R), they denote, as follows: However, the ranking position method [30] only focuses on the position factor, but it ignores the subordinate utilities and the relative importance of the subordinate order, which do not reflect the real performance of the alternative in the subordinate ranking. In other words, this method only involves the ranking matrix, but not both the ranking matrix and the utility value matrix. A simple example can illustrate the shortcomings of this method: there are three alternatives, which are denoted as A 1 , A 2 , A 3 . Their utility value matrix of three models (RS, RP, FMF) is denoted as: Obviously, the ranking lists can be given, as follows: The computation results of each alternatives denote, as follows: Thus, the final ranking result is: A 1 = A 2 = A 3 , if the relative importance of the three sub-models is the same. Actually, the performance of A 1 and A 3 in RP is similar, but A 1 is much better than A 3 in RS and FMF. Comprehensively speaking, it must be more reasonable that A 1 > A 3 , which is more in line with our cognition. When considering that the PDHFLRP model is the cost-based attribute, while PDHFLRS and PDHFLFMF models are benefit-based attribute, the score function of the ranking position method is thus improved by an arithmetic weighted operator, as follows: By the Equation (27), we derive a new result of 141. This result defeats some unreasonable defects.

The Procedure of the PDHFL-MULTIMOORA Method
The PDHFL-MULTIMOORA methodology will be performed by the following step in order to clearly illustrate the structure of our proposed method and facilitate the reader's application, as shown in Figure 2.
Step 1. Collect the linguistic evaluation information from expert E T (T = 1, . . . , K) on alternative A i (i = 1, . . . , m) with respect to attribute C j ( j = 1, . . . , n) to derive the PDHFLTE d T ij = p T ij,1 , . . . p T ij,l , . . . , p T ij,g , and then establish the individual decision matrix D T = d T ij m×n .
Step 2. Calculate the weight of expert that detailed operational procedures have been described in Section 3. On basis of obtaining expert weight vector ω = ω 1 , . . . , ω L , . . . , ω K , we aggregate the individual decision matrix into the group decision matrix at step 3.
Step 3. Let D = d ij m×n be the group decision matrix consisting of elements d ij = p ij,1 , . . . p ij,l , . . . , p ij,g . By Equation (3), we have: where p ij,l is given by: Step 4. Normalize the group decision matrix D = d ij m×n into D N = d N ij m×n by Equation (23), and then go to step 5.
Step 5. Gain the utility function U T (X) of each expert and the constant vector w c = (w 1 , . . . , w n ) with respect to attribute set, we use Equation (22) to obtain the group utility function U(X), and then the weight vector of attribute W = (w 1 , . . . , w n ) is derived by the following formula: Noteworthily, the value of cost type attribute should be converted into 1 − d N ij .
Step 6. Compute the utility value U 1 (A i ) of PDHFLRS model by Equation (24), and then obtain the first subordinate ranking set: . . , A i|i⊂minU 1 (A i ) and utility value set: Step 7. Calculate the utility value U 2 (A i ) of PDHFLRP model by Equation (25), thus obtaining the second subordinate ranking set: . . , A i|i⊂maxU 2 (A i ) and utility value set: Step 8. Derive the utility value U 3 (A i ) of PDHFLFMF model by Equation (26), hence we receive the third subordinate ranking set: . . , A i|i⊂minU 3 (A i ) and utility value set: Step 9. Based on the utility values and subordinate ranking sets of three model, establish the utility value matrix D(U) and ranking matrix D(R). After that, we apply the improved ranking position rule as Equation (27)

Case Description
In China, express industry has shown explosive growth, owing to the tremendous success of ecommerce. However, at the same time, a lot of resources are wasted and the contradiction between environment and society is increasingly prominent. As the government and the public increasing attention to sustainable development, the enterprise's actions in the environment are playing a crucial role. Thus, whether or not the express industry has the economic, social, and environmental capabilities in the process of reverse logistics becomes the key to solve the problem. In fact, it is According to presentation on algorithmic steps, Figure 2 shows the general framework.

Case Description
In China, express industry has shown explosive growth, owing to the tremendous success of e-commerce. However, at the same time, a lot of resources are wasted and the contradiction between environment and society is increasingly prominent. As the government and the public increasing attention to sustainable development, the enterprise's actions in the environment are playing a crucial role. Thus, whether or not the express industry has the economic, social, and environmental capabilities in the process of reverse logistics becomes the key to solve the problem. In fact, it is commonplace for enterprises to entrust reverse logistics management to third-party reverse logistics providers (3PRLP) in order to ensure that the reverse flow of products is handled in an efficient way [68], because enterprises can focus on their core business by outsourcing reverse logistics business to appropriate suppliers. If not, their image and interests will be greatly damaged. Hence, sustainable 3PRLP selection bears the vital importance of express companies, especially the sustainability.
The process of choosing 3PRLP is very complex and full of uncertainty in order to maximize profits and minimize risks. Thus, a good decision-making method is particularly important. To data, different types of methods have been applied to this field, but the most concerned and special method is MADM, because of its effectiveness in dealing with complex problems. Table 2 summarizes previous literatures on 3PRLP selection methods.  [79,80] Existing papers demonstrate that the attributes of selecting a sustainable 3PRLP are mainly concentrated on four aspects: society, economy, risk, and environment. Further, we summarize the core attributes in recent studies, as shown in Table 3.
However, existing literatures hardly concentrate on the expert evaluation system. In practice, increasing enterprises tend to invite a group of experts to make an evaluation instead of single individual. Moreover, due to the complexity and fuzziness of the problem, experts cannot start with an accurate description. Meanwhile, the objective data of quantitative indicators are not accurately available, only being replaced by experts' linguistic expertise. Therefore, the proposed MAGDM method, which focuses on processing linguistic variables, in this paper, is suitable and reasonable for the sustainable 3PRLPs selection issue.  [68,72,77,87] Suppose that there are four attributes, namely, economic ability (C 1 ), environmental ability (C 2 ), risk level(C 3 ) [88], and social ability (C 4 ).
Subsequently, four experts E = (E 1 , E 2 , E 3 , E 4 ) are invited to participate in assessing four providers X = (A 1 , A 2 , A 3 , A 4 ). The LTS for assessing the providers with respect to the four attributes is defined as S = {s 1 = VeryBad(VB), s 2 = Bad(B), s 3 = SlightBad(SB), s 4 = Medium(M), s 5 = SlightGood(SG), s 6 = Good(G), s 7 = VeryGood(VG)}, and then experts provide their linguistic evaluations as PDHFLTSs form respectively. The results are shown as the following Tables 4-7. Table 4. Possibility distribution hesitant fuzzy linguistic decision matrix by expert E 1 .       We calculate the case according to the procedures of the proposed method. As the first step has already been done in previous part, we start with step 2.
Step 2. Based on the linguistic decision matrix of each expert, the weight of expert can be determined by Equation (21), where the consensus degree and hesitant degree are computed by Equations (17) and (20). The weight vector, the consensus degree vector, and hesitant degree vector of experts are established, respectively, as: Step 3. Calculate the collective expressions by Equation (3), and then the group decision matrix is denoted by: Step 4. By Equation (23), derive the normalized group decision matrix as: Step 5. Suppose that the utility functions of expert E 1 , E 2 , E 3 , E 4 are given by, respectively: According to Equation (22), the group utility function is: where ω i (i = 1, 2, 3, 4) denotes the weight of expert E i ,r = 0.5. Subsequently, we have variable weights, as follows: where w c j = 0.25 is a constant weight of attribute C j and the normalized value of cost attribute should be updated to 1 − d N ij . Subsequently, W(X) is expressed as: Step 6. Based on the weight of attributes and normalized value of group decision matrix, Obtain the utility value of PDHFLTRS by Equation (24): Thus, the first subordinate ranks are Step 7. Derive the utility value of PDHFLRP with Equation (25): The second subordinate orders are shown as: Step 8. By analogical computation with Equation (26), the utility value of PDHFLFMF is: The third subordinate ranks are: Step 9. Finally, the IRPM value can be computed by Equation (27) Afterwards, the final ranks are In conclusion, the sustainable provider 1 is the best option for selection.

Comparative Analysis with the PDHFL-MULTIMOORA Method Based on Existing Weight Methods
In this part, we make a comparison of the proposed variable weight-based weight-determining method and other existing methods.
Firstly, the entropy-weighting method is one of the most widely used weight-determining method; its basic idea is that the more information energy contained in an attribute, the greater the weight will be allocated. Since the entropy measure of hesitant fuzzy linguistic elements [89], defined by Farhadinia, the weights of attribute can be obtained by following formulas: The PDHFL-MULTIMOORA with the entropy-weighting method is summarized, as follows: Step 1. For the above case, the content of procedures is same as Step 1-3 of the PDHFL-MULTIMOORA with the variable weight method.
Step 3. Reference Step 4 and 6-7 of the PDHFL-MULTIMOORA with the variable weight method.
Subsequently, we get the new utility value and subordinate ranks matrix of three models, as follows: Step 4. Calculate the final results by Equation (27): Secondly, the standard deviation method's main idea is similar to the entropy weight method. Usually, the larger the standard deviation of an attribute, the greater the variation of the attribute value, the greater the amount of information provided, and the greater its weight. The PDHFL-MULTIMOORA with the standard deviation method is thus summarized, as follows: Step 1. It is same as Step 1-4 of the proposed PDHFL-MULTIMOORA method.
Step 2. Derive the weight of attributes by Equation (34), while using standard deviation formula Equation (33).
Hence, we get: σ = (σ 1 , σ 2 , σ 3 , σ 4 ) = (0.125, 0.087, 0.122, 0.174) ω = (ω 1 , ω 2 , ω 3 , ω 4 ) = (0.246, 0.171, 0.24, 0.343) Step 3. Obtain the utility value and subordinate orders: Step 4. Derive final ranks by Equation (27) Thirdly, derive the results by the PDHFL-MULTIMOORA with the same weight method, where all attributes' weights are equal. After an analogical analysis, we can receive: Thus, the final orders are From above weighting method, the final results that were obtained by the proposed method are the same as entropy-weighting and same weight methods, while they are different from the standard deviation method. However, there is little difference in Table 8. Hence, this can prove the effectiveness and rationality of our proposed weighting approach, to a large extent. However, we can also find that the results are sensitive to attributes' weights, that is, different weights of attributes may lead to different results, as Figure 3 clearly shows, only the weight of the variable weight method is inconstant with respect to the same attribute, while the other methods are constant. Through the example of the Section 4, we know that, if the weight is unchanged, it will probably produce unreasonable results. Therefore, it can be concluded that the variable weight method can maintain the reasonableness of the results to a greater extent than the other weight methods. Table 8. The comparison of mentioned weight methods.

Comparative Analysis with the PDHFL-TOPSIS and the PDHFL-VIKOR Methods
This section mainly works on the reliability of the proposed method by making a comparison with the PDHFL-TOPSIS and the PDHFL-VIKOR methods.
Motivated by Wu [15], The steps of the TOPSIS based approach to solve the above case are

Comparative Analysis with the PDHFL-TOPSIS and the PDHFL-VIKOR Methods
This section mainly works on the reliability of the proposed method by making a comparison with the PDHFL-TOPSIS and the PDHFL-VIKOR methods.
Motivated by Wu [15], The steps of the TOPSIS based approach to solve the above case are shown, as below.
Step 5. Find the positive ideal solution (PIS) p + j and negative ideal solution p − j regarding attribute C j by following formulas: where C j , j = 1, . . . , g belong to beneficial attributes and C j , j = g + 1, . . . , n belong to cost attributes. Similar to 15, p + j and p − j are defined by Definition 6, and then we have: Step 6. Compute the distance of each provider from the PIS and NIS: Subsequently, we get: Step 7. Determine the relative closeness of each provider to ideal solution by following-up equation: Afterwards, we have: RC 1 = 0.883, RC 2 = 0.246, RC 3 = 0.526, RC 4 = 0.528. The final ranks are thus denoted as: Analogically, solve the above case based on VIKOR, and the steps are expressed as: Step 1. Consult Step 1 of the PDHFL-TOPSIS method.
Step 6. Determine the values of S i and R i by: Then we get: Step 7. Compute the values of Q i by: where θ denotes the proportion of group utility maximization strategy, 1 − θ denotes the proportion of individual regret. Besides, If we pursue the maximization of group utility and the minimization of individual regret, then put θ = 0.7, and we have: Step 8. Rank all of the schemes by sorting S i , R i , Q i in ascending order, as follows: Step 9. According to verification by conditions 1 and 2 in VIKOR model that was proposed by Wu [15], the providers 1, 3, and 4 are compromise solutions and final results can be denoted as: The case of θ = 0.3 is shown in Table 9.
All three methods show that sustainable 3PRLP 1 is the best option in Table 9. Hence, the original express incorporation can cooperate with the provider 1 to maximize its value. However, there are certain different values for each alternative obtained by different MAGDM methods, we find that the utility values of provider 3 and 4 are close to each other in the TOPSIS and VIKOR based methods. The TOPSIS based method illustrates that provider 3 is worse than provider 4, because of the ignorance of individual regret, and just provider 3 performs badly with respect to attribute C 4 . The VIKOR based method considers this aspect that the greater value of θ, the greater the proportion of individual regret, the better perform provider 3 is than provider 4, but it is difficult for experts to subjectively determine the threshold θ in practical application. If we revise the weight of individual regret, the results of them are likely to change. This situation can be seen in Table 9. From the theoretical level, the operation essence of the TOPISIS based method is arithmetic weighting aggregation, which is similar to operation and sort results of PDHFLRS. Although the VIKOR based method has considered individual regret, the result is very sensitive to its proportion, which is difficultly determined. The PDHFL-MULTIMOORA method considers both the group utility values and individual regret values, which controlled by three models, and results will depend on utility values and subordinate ranks by introducing improved ranking position method. The utility value of the improved ranking position method, from the theoretical level of speaking, is equivalent to controlling the relative importance of three models. This approach overcomes the shortcoming of VIKOR based approach's threshold θ, which is difficult to estimate.
Further, we see no change in the result of the PDHFL-MULTIMOORA approach, which proves the robustness and rationality of presented approach in Table 10. In detail, operation 1 means that normalized utility value of the PDHFLRS increased by 50%; Operation 2 means that the normalized utility value of the PDHFLRP increased by 50%; Operation 3 means that normalized utility value of the PDHFLFMF increased by 50%; Operation 4 means that normalized utility value of the PDHFLRS decreased by 50%; Operation 5 means that normalized utility value of the PDHFLRP decreased by 50%; and, Operation 6 means that normalized utility value of the PDHFLFMF decreased by 50%.

Conclusions and Further Directions
The PDHFLTS is a suitable tool for expressing complex linguistic information to cope with MAGDM problems. This paper combines PDHFLT with the MULTIMOORA method for deriving solution of sustainable 3PRLP selection problem. The major contributions can be summarized, as follows: (1) A novel method is presented to compute the weight of expert both considering hesitant degree and consensus degree, and a covariance-based correlation coefficient measure is defined to determine the consensus degree, which is proven to be able to simplify calculations and it is specifically designed for hesitant fuzzy linguistic environment based on possibility distribution. A hesitant degree formula under PDHFL context is also introduced.
(2) The incentive-type variable weight method considering group utility function is developed to calculate the weight of attributes in which group utility derived by aggregation of individual utility. Some properties of utility function, combined with the variable weight method, are examined. (3) A new aggregation operator, named the improved ranking position method, is introduced into the MULTIMOORA method. When compared to existing MAGDM methods, the proposed approach has the following advantages: (a). the HFLTS with possibility distribution can carry more abundant information of experts. (b). the weight of attributes avoids some unreasonable situations, the weight of expert reflects each individual's uncertainty and closeness to the group. (c). the final ranks retain robustness and effectiveness, because three models with different effects are applied to calculate the group decision matrix. (4) A sustainable 3PRLP selection problem is solved by the proposed method to verify its rationality and suitability. Two comparative analyses are made to illustrate the difference of presented methods and other existing methods. (5) In the future, the research can focus on several aspects: the MULTIMOORA with improved ranking position method will be combined with another useful information expression tool to promote it much wider. Secondly, the PDHFLTs based on other MAGDM methods will be introduced in later work. What is more, aggregation of heterogeneous evaluation data of sustainable 3PRLP selection could also continue to mine deeply in the field of MAGDM background. Funding: This research received no external funding.