*Article* **Research on Hybrid Multi-Attribute Three-Way Group Decision Making Based on Improved VIKOR Model**

**Jiekun Song \*, Zeguo He \*, Lina Jiang, Zhicheng Liu and Xueli Leng**

School of Economics and Management, China University of Petroleum, Qingdao 266580, China

**\*** Correspondence: songjiekun@upc.edu.cn (J.S.); z21080318@s.upc.edu.cn (Z.H.)

**Abstract:** In the era of internet connection and IOT, data-driven decision-making has become a new trend of decision-making and shows the characteristics of multi-granularity. Because three-way decision-making considers the uncertainty of decision-making for complex problems and the cost sensitivity of classification, it is becoming an important branch of modern decision-making. In practice, decision-making problems usually have the characteristics of hybrid multi-attributes, which can be expressed in the forms of real numbers, interval numbers, fuzzy numbers, intuitionistic fuzzy numbers and interval-valued intuitionistic fuzzy numbers (IVIFNs). Since other forms can be regarded as special forms of IVIFNs, transforming all forms into IVIFNs can minimize information distortion and effectively set expert weights and attribute weights. We propose a hybrid multiattribute three-way group decision-making method and give detailed steps. Firstly, we transform all attribute values of each expert into IVIFNs. Secondly, we determine expert weights based on interval-valued intuitionistic fuzzy entropy and cross-entropy and use interval-valued intuitionistic fuzzy weighted average operator to obtain a group comprehensive evaluation matrix. Thirdly, we determine the weights of each attribute based on interval-valued intuitionistic fuzzy entropy and use the VIKOR method improved by grey correlation analysis to determine the conditional probability. Fourthly, based on the risk loss matrix expressed by IVIFNs, we use the optimization method to determine the decision threshold and give the classification rules of the three-way decisions. Finally, an example verifies the feasibility of the hybrid multi-attribute three-way group decision-making method, which provides a systematic and standard solution for this kind of decision-making problem.

**Keywords:** hybrid multi-attribute; three-way group decision making; VIKOR model; grey correlation analysis; interval-valued intuitionistic fuzzy numbers

**MSC:** 90B50; 03E72

#### **1. Introduction**

With the rapid popularization of the internet and the internet of things, the generation and collection speed of various decision-making data in economic production and life is rapidly increasing. Due to the limitations of data collection technology and expert judgment [1,2], the decision-making data show the characteristics of incompleteness, uncertainty, incongruity, fuzziness and hesitation [3,4]. For this kind of decision-making problem with complex decision data and uncertain evaluation information, the traditional optimization mechanism model based on function relationship becomes more difficult in decision analysis and problem-solving. In fact, there is a large amount of effective decision information hidden in the decision data. Based on the existing decision data, we use scientific data processing technology to objectively analyze and evaluate them and transform them into effective decision indicators and knowledge, which can provide reliable and reasonable suggestions and decision support for decision-makers. This data-driven decision-making has become a new trend in modern decision-making [5–7].

**Citation:** Song, J.; He, Z.; Jiang, L.; Liu, Z.; Leng, X. Research on Hybrid Multi-Attribute Three-Way Group Decision Making Based on Improved VIKOR Model. *Mathematics* **2022**, *10*, 2783. https://doi.org/10.3390/ math10152783

Academic Editors: Wen Zhang, Xiaofeng Xu, Jun Wu and Kaijian He

Received: 18 June 2022 Accepted: 1 August 2022 Published: 5 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Multi-attribute decision making (MADM) is the most common decision-making problem in practice. Objects are evaluated and sorted according to the comprehensive performance of multi-attribute. In order to reflect the uncertainty of human cognition, Zadeh proposed fuzzy set theory [8], linguistic variable [9–11] and possibility measure and introduced them into the MADM problem [12]. Nowadays, fuzzy set theory has been developed and produced in many forms. Because the fuzzy set only has a membership index of fuzzy objects, it is difficult to describe people's subjective understanding of fuzzy concepts completely. Atanassov proposed intuitive fuzzy sets by adding a non-membership degree and hesitation degree to the relationship between objects and sets [13], which can more truly reflect the subject's understanding of the fuzzy nature of objective things when expressing uncertain information [14]. Since the membership degree and non-membership degree may also be uncertain, Atanassov and Gargov further extended them into the form of interval numbers and proposed the interval-valued intuitive fuzzy set (IVIFS) [15]. Intuitionistic fuzzy sets and interval-valued intuitionistic fuzzy sets have been introduced into many traditional decision models to solve MADM problems, such as the combination with AHP (Analytic Hierarchy Process) [16,17], TOPSIS (Technique for Order Preference by Similarity to an Ideal Solution) [18,19], VIKOR (VlseKriterijumska Optimizacija I Kompromisno Resenje) [20], ELECTRE (Elimination et Choice Translating Reality) [21,22], PROMETHEE (preference ranking organization methods for enrichment evaluations) [23], etc.

We can sort and select different schemes by MADM. However, in practice, we often encounter the following situation: we plan to select the top 10 of the 15 suppliers as the access suppliers, but after a comprehensive evaluation, the evaluation results of the ninth to 11th suppliers may be slightly different. There are certain risks in accepting or rejecting these three suppliers, and further field visits may be required. This means that the 15 suppliers can be divided into three types, i.e., accepted, rejected and to be further determined. The three-way decision theory can make exactly three kinds of decisions in this situation. The three-way decision is a new theory proposed by Yao on the basis of the rough set theory. A rough set applies the lower and upper approximations of equivalence relation to divide the universe of objects into three pair-wise disjoint regions, i.e., positive region, negative region and boundary region [24]. In a classical rough set, the positive region and the associated positive rules are the focus of attention, as these rules produce consistent acceptance and rejection decisions. However, the decisions are made without any tolerance of uncertainty, which is too strict for dealing with incomplete and noisy data and is insensitive to the cost of classification errors. In order to overcome these deficiencies, Yao et al. introduced the Bayesian minimum risk decision theory and proposed the decision-theoretic rough set models [25], which can allow a tolerance of inaccuracy in lower and upper approximations and define three regions including probabilistic positive, boundary, and negative regions. However, there is difficulty in interpreting rules in the decision-theoretic rough set models. For example, an object in the probabilistic positive region does not certainly belong to the decision class, but with a high probability (i.e., the probability value is above a certain threshold) [26], so a rule from the probabilistic positive region may be uncertain and nondeterministic. In order to better interpret the rules qualitatively, Yao et al. introduced the notion of three-way decision rules [27]. Positive, negative, and boundary rules are constructed from the corresponding regions, and they represent the results of a three-way decision of acceptance, rejection, or abstaining. In addition, the decisions of acceptance and rejection are made with certain levels of tolerance for errors, which reflects the cost sensitivity of decision-makers to incorrect classification decisions. Obviously, the semantics of three-way decisions are consistent with the thinking of human beings in dealing with complex decision-making problems. At present, three-way decision has been widely used in the field of MADM and produced good results [28–31]. In reality, the various indicators of evaluation objects have different expression forms. Some indicators can be expressed in exact real numbers, some can be expressed in uncertain interval numbers, some can be expressed as the fuzzy values of qualitative linguistic variables, and some can be expressed in the forms of fuzzy numbers, intuitive fuzzy numbers, IVIFNs, etc. Therefore, it is of great

significance to discuss the three-way decisions under a hybrid multi-attribute environment, especially in the case of attributes represented by intuitionistic fuzzy numbers or IVIFNs with more fuzzy information.

The representative studies on the three-way decisions under intuitionistic fuzzy or interval-valued intuitionistic fuzzy multi-attribute environments are shown in Table 1.

**Table 1.** The representative three-way decision methods under intuitionistic fuzzy or interval-valued intuitionistic fuzzy multi-attribute environment.


The main methods for determining conditional probability in three-way decisions include TOPSIS [32], grey correlation analysis [33] and VIKOR [34]. Two methods are used to determine the decision thresholds: one is to use the optimization method to determine the thresholds based on the subjective risk loss matrix [40,41]; the second is to determine the losses based on the preference coefficient and the distance from the ideal points and then use the formula derived from Bayesian decision to determine the thresholds [33]. In addition, some scholars put forward the method of weight determination based on deviation [34], and some scholars put forward the method of grey correlation analysis to determine the weights of experts in group decision-making [39].

The above literature provides a good foundation for this study. However, the existing studies still have the following contents that may be deepened. Firstly, there is a lack of discussion on the hybrid multi-attribute three-way decision, even the study on the intervalvalued intuitionistic fuzzy three-way decision is relatively lacking. Secondly, there are few discussions about expert weight and attribute weight in the interval-valued intuitionistic fuzzy three-way group decisions. In fact, the interval-valued intuitionistic fuzzy group decision matrix contains a lot of information. It is of great significance to make effective use of the information and give the scientific weights of experts and attributes for decision results. Thirdly, the determination method of conditional probability in the three-way decision can be further improved. For example, the advantages of VIKOR, TOPSIS and grey correlation analysis can be fully integrated to form a grey correlation improved VIKOR model, which can give the conditional probability more objectively. In order to make up for the above deficiencies, we will discuss the hybrid multi-attribute three-way group decisionmaking method. The attribute values of different forms are unified into IVIFNs with the least information distortion. Based on the IVIFNs group decision matrix, the expert weight and attribute weight are determined. Then the conditional probability is determined by using the improved VIKOR model based on grey correlation, and the three-way decision rules can be formed by comparing with the threshold pair based on optimization.

The rest of this paper is organized as follows. Section 2 proposes research preliminaries, including interval-valued intuitionistic fuzzy sets and three-way decisions. Section 3 proposes a hybrid multi-attribute three-way group decision method based on an improved VIKOR model. Section 4 provides an illustrative example to verify the validity of the method. Section 5 summarizes the conclusions of this study.

#### **2. Preliminaries**

*2.1. Interval-Valued Intuitionistic Fuzzy Sets*

**Definition 1 [15].** *Let X be a non-empty set and an IVIFS A in X is expressed as follows:* 2

$$\widetilde{A} = \left\{ \langle \mathbf{x}, \,\widetilde{\boldsymbol{\mu}}\_{\widetilde{A}}(\mathbf{x}), \,\widetilde{\boldsymbol{v}}\_{\widetilde{A}}(\mathbf{x}) \rangle \middle| \mathbf{x} \in X \right\} = \left\{ \langle \mathbf{x}, \,\left[ \boldsymbol{\mu}\_{\widetilde{A}}^{L}(\mathbf{x}), \,\boldsymbol{\mu}\_{\widetilde{A}}^{R}(\mathbf{x}) \right], \,\left[ \boldsymbol{v}\_{\widetilde{A}}^{L}(\mathbf{x}), \,\boldsymbol{v}\_{\widetilde{A}}^{R}(\mathbf{x}) \right] \rangle \middle| \mathbf{x} \in X \right\} \tag{1}$$

*where, μ<sup>L</sup> <sup>A</sup>*2(*x*) *and <sup>μ</sup><sup>R</sup> <sup>A</sup>*2(*x*)*, respectively, represent the upper and lower boundaries of the membership degree <sup>μ</sup>*2*A*2(*x*) *of the element x in X belonging to <sup>A</sup>*2*; <sup>v</sup><sup>L</sup> <sup>A</sup>*2(*x*) *and <sup>v</sup><sup>R</sup> <sup>A</sup>*2(*x*)*, respectively, represent the upper and lower boundaries of the non-membership degree <sup>v</sup>*2*A*2(*x*) *of the element x that belong to <sup>A</sup>*2*. For each x* ∈ *X, it satisfies the conditions:* 0 ≤ *μ<sup>L</sup> <sup>A</sup>*2(*x*) <sup>≤</sup> *<sup>μ</sup><sup>R</sup> <sup>A</sup>*2(*x*) <sup>≤</sup> <sup>1</sup>*,* <sup>0</sup> <sup>≤</sup> *<sup>v</sup><sup>L</sup> <sup>A</sup>*2(*x*) <sup>≤</sup> *<sup>v</sup><sup>R</sup> <sup>A</sup>*2(*x*) <sup>≤</sup> <sup>1</sup>*,* 0 ≤ *μ<sup>R</sup> <sup>A</sup>*2(*x*) + *<sup>v</sup><sup>R</sup> <sup>A</sup>*2(*x*) <sup>≤</sup> <sup>1</sup>*,* <sup>∀</sup>*<sup>x</sup>* <sup>∈</sup> *<sup>X</sup>*.

**Definition 2 [15].** *For an IVIFS A, the hesitation degree of element x in* 2 *A is:* 2

$$\pi\_{\vec{A}}(\mathbf{x}) = 1 - \left[\boldsymbol{\pi}\_{\vec{A}}(\mathbf{x}) - \boldsymbol{\overline{\boldsymbol{\nu}}}\_{\vec{A}}(\mathbf{x}) = \left[\boldsymbol{\pi}\_{\vec{A}}^{\mathsf{L}}(\mathbf{x}), \ \boldsymbol{\pi}\_{\vec{A}}^{\mathsf{R}}(\mathbf{x})\right] = \left[1 - \boldsymbol{\mu}\_{\vec{A}}^{\mathsf{R}}(\mathbf{x}) - \boldsymbol{\nu}\_{\vec{A}}^{\mathsf{R}}(\mathbf{x}), \ 1 - \boldsymbol{\mu}\_{\vec{A}}^{\mathsf{L}}(\mathbf{x}) - \boldsymbol{\nu}\_{\vec{A}}^{\mathsf{L}}(\mathbf{x})\right] \tag{2}$$

**Definition 3 [42].** *For an IVIFS <sup>A</sup>*2*, the fuzzy degree* <sup>Δ</sup>*A*2(*x*) *of element x belonging to <sup>A</sup>*<sup>2</sup> *is given as follows:*

$$\Delta\_{\bar{A}}(\mathbf{x}) = \sqrt{\frac{\left(\Delta\_{\bar{A}}^{L}(\mathbf{x})\right)^{2} + \left(\Delta\_{\bar{A}}^{R}(\mathbf{x})\right)^{2}}{2}} \tag{3}$$

*where:*

$$
\Delta^L\_{\bar{A}}(\mathbf{x}) = \left| \mu^L\_{\bar{A}}(\mathbf{x}) - v^L\_{\bar{A}}(\mathbf{x}) \right| , \Delta^R\_{\bar{A}}(\mathbf{x}) = \left| \mu^R\_{\bar{A}}(\mathbf{x}) - v^R\_{\bar{A}}(\mathbf{x}) \right| \tag{4}
$$

**Definition 4 [15].** *The complement of an IVIFS A is given as follows:* 2

$$\tilde{A}^{\mathbb{C}} = \left\{ \langle \mathbf{x}, \overline{v}\_{\tilde{A}}(\mathbf{x}), \overline{\mu}\_{\tilde{A}}(\mathbf{x}) \rangle | \mathbf{x} \in X \right\} \tag{5}$$

**Definition 5 [15].** *Given three IVIFNs* <sup>2</sup>*<sup>α</sup>* <sup>=</sup> ([*a*, *<sup>b</sup>*], [*c*, *<sup>d</sup>*])*,* <sup>2</sup>*α*<sup>1</sup> <sup>=</sup> ([*a*1, *<sup>b</sup>*1], [*c*1, *<sup>d</sup>*1]) *and* <sup>2</sup>*α*<sup>2</sup> <sup>=</sup> ([*a*2, *b*2], [*c*2, *d*2])*, their basic operations are summarized as follow:*

*(1)* <sup>2</sup>*α*<sup>1</sup> <sup>+</sup> <sup>2</sup>*α*<sup>2</sup> <sup>=</sup> ([*a*<sup>1</sup> <sup>+</sup> *<sup>a</sup>*<sup>2</sup> <sup>−</sup> *<sup>a</sup>*1*a*2, *<sup>b</sup>*<sup>1</sup> <sup>+</sup> *<sup>b</sup>*<sup>2</sup> <sup>−</sup> *<sup>b</sup>*1*b*2], [*c*1*c*2, *<sup>d</sup>*1*d*2])*; (2)* <sup>2</sup>*α*<sup>1</sup> <sup>⊗</sup> <sup>2</sup>*α*<sup>2</sup> <sup>=</sup> ([*a*1*a*2, *<sup>b</sup>*1*b*2], [*c*<sup>1</sup> <sup>+</sup> *<sup>c</sup>*<sup>2</sup> <sup>−</sup> *<sup>c</sup>*1*c*2, *<sup>d</sup>*<sup>1</sup> <sup>+</sup> *<sup>d</sup>*<sup>2</sup> <sup>−</sup> *<sup>d</sup>*1*d*2])*; (3) <sup>λ</sup>*2*<sup>α</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> (<sup>1</sup> <sup>−</sup> *<sup>a</sup>*) *<sup>λ</sup>*, 1 − (1 − *b*) *λ* , " *<sup>c</sup>λ*, *<sup>d</sup><sup>λ</sup>*# , *λ* ≥ 0*; (4)* <sup>2</sup>*α<sup>λ</sup>* <sup>=</sup> " *aλ*, *bλ*# , 1 − (1 − *c*) *<sup>λ</sup>*, 1 − (1 − *d*) *λ* , *<sup>λ</sup>* <sup>≥</sup> 0.

**Definition 6 [42].** *Let IVIFS(X) be the set of all IVIFSs in X*, *a real-valued function E: IVIFS(X) [0, 1] is called an entropy measure for IVIFSs if it satisfies the following axiomatic requirements:*


*(4) For a constant* a *in (0, 1), let* Δ*<sup>L</sup> <sup>a</sup> ,* Δ*<sup>R</sup> <sup>a</sup> , π<sup>L</sup> <sup>a</sup> and π<sup>R</sup> <sup>a</sup> be the sets of all IvIFSs* . *x*, *<sup>μ</sup>*2*A*2(*x*), *<sup>v</sup>*2*A*2(*x*) / *in X with* Δ*<sup>L</sup> <sup>A</sup>*2(*x*) = *a,* <sup>Δ</sup>*<sup>R</sup> <sup>A</sup>*2(*x*) = *a, <sup>π</sup><sup>L</sup> <sup>A</sup>*2(*x*) = *<sup>a</sup> , <sup>π</sup><sup>R</sup> <sup>A</sup>*2(*x*) = *<sup>a</sup> , respectively. <sup>E</sup> A*2 *is strictly monotone decreasing with respect to* Δ*<sup>L</sup> <sup>A</sup>*2(*x*) *on* <sup>Δ</sup>*<sup>L</sup> <sup>a</sup> and* Δ*<sup>R</sup> <sup>A</sup>*2(*x*) *on* <sup>Δ</sup>*<sup>R</sup> <sup>a</sup> respectively and is strictly monotone increasing with respect to π<sup>L</sup> <sup>A</sup>*2(*x*) *on <sup>π</sup><sup>L</sup> <sup>a</sup> and π<sup>R</sup> <sup>A</sup>*2(*x*) *on <sup>π</sup><sup>R</sup> <sup>a</sup> , respectively.*

**Definition 7.** *In [43], for an IVIFS A*2 *in X, X = {x1, x2,* ... *, xn}, the authors define the following entropy function:*

$$E\left(\bar{A}\right) = \frac{1}{n} \sum\_{i=1}^{n} \frac{2 - 2\left(\Lambda\_{\bar{A}}(\mathbf{x}\_{i})\right)^{3} + \left(\pi\_{\bar{A}}^{L}(\mathbf{x}\_{i})\right)^{3} + \left(\pi\_{\bar{A}}^{R}(\mathbf{x}\_{i})\right)^{3}}{4} \tag{6}$$

*It is not difficult to prove that the above entropy function satisfies the axiomatic condition of interval-valued intuitionistic fuzzy entropy in Definition 6.*

**Definition 8 [44].** *Given two IVIFSs* 3*A and B*2 *in X, X = {x1, x2,* ... *, xn}, the cross entropy of them is defined as follows:*

$$D\left(\widehat{A}\_{r}^{\top},\widehat{B}\right) = \sum\_{i=1}^{n} \left[ d\_{\widehat{A}}\left(\mathbf{x}\_{i}\right) \ln \frac{d\_{\widehat{A}}\left(\mathbf{x}\_{i}\right)}{\frac{1}{2}\left(d\_{\widehat{A}}\left(\mathbf{x}\_{i}\right) + d\_{\widehat{B}}\left(\mathbf{x}\_{i}\right)\right)} + \left(1 - d\_{\widehat{A}}\left(\mathbf{x}\_{i}\right)\right) \ln \frac{1 - d\_{\widehat{A}}\left(\mathbf{x}\_{i}\right)}{\frac{1}{2}\left(2 - d\_{\widehat{A}}\left(\mathbf{x}\_{i}\right) - d\_{\widehat{B}}\left(\mathbf{x}\_{i}\right)\right)} \right] \tag{7}$$

*where:*

$$d\_{\stackrel{\frown}{A}}(\mathbf{x}\_{i}) = \frac{1}{2} \left[ \frac{\mu\_{\stackrel{\frown}{A}}^{L}(\mathbf{x}\_{i}) + \mu\_{\stackrel{\frown}{A}}^{R}(\mathbf{x}\_{i})}{2} + 1 - \frac{\upsilon\_{\stackrel{\frown}{A}}^{L}(\mathbf{x}\_{i}) + \upsilon\_{\stackrel{\frown}{A}}^{R}(\mathbf{x}\_{i})}{2} \right] = \frac{\mu\_{\stackrel{\frown}{A}}^{L}(\mathbf{x}\_{i}) + \mu\_{\stackrel{\frown}{A}}^{R}(\mathbf{x}\_{i}) + 2 - \upsilon\_{\stackrel{\frown}{A}}^{L}(\mathbf{x}\_{i}) - \upsilon\_{\stackrel{\frown}{A}}^{R}(\mathbf{x}\_{i})}{4} \tag{8}$$

$$d\_{\bar{B}}\left(\mathbf{x}\_{i}\right) = \frac{1}{2} \left[ \frac{\mu\_{\bar{B}}^{L}(\mathbf{x}\_{i}) + \mu\_{\bar{B}}^{R}(\mathbf{x}\_{i})}{2} + 1 - \frac{\upsilon\_{\bar{B}}^{L}(\mathbf{x}\_{i}) + \upsilon\_{\bar{B}}^{R}(\mathbf{x}\_{i})}{2} \right] = \frac{\mu\_{\bar{B}}^{L}(\mathbf{x}\_{i}) + \mu\_{\bar{B}}^{R}(\mathbf{x}\_{i}) + 2 - \upsilon\_{\bar{B}}^{L}(\mathbf{x}\_{i}) - \upsilon\_{\bar{B}}^{R}(\mathbf{x}\_{i})}{4} \tag{9}$$

*Obviously,* 0 ≤ *D* 3*A* , *B*2 ≤ *n* ln 2*, and D* 3*A* , *B*2 <sup>=</sup> <sup>0</sup>*, if and only if <sup>μ</sup>*2*A*2(*x*) = *<sup>μ</sup>*2*B*2(*x*)*, <sup>v</sup>*2*A*2(*x*) = *<sup>v</sup>*2*B*2(*x*)*. The cross entropy can also be called the relative entropy or divergence measure,* *which indicates the discrimination degree of IVIFS* 3*A fromB*2*. Since the cross entropy formula does not satisfy the symmetry, we rewrite it as follows:*

$$D^\*\left(\overleftarrow{A},\,^\*\overline{B}\right) = D\left(\overleftarrow{A},\,^\*\overline{B}\right) + D\left(\overline{B},\,^\*\overline{A}\right) \tag{10}$$

*It is not difficult to prove that the following relationships hold: D*<sup>∗</sup> 3*A* , *B*2 = *D*<sup>∗</sup> 3*A C* , *<sup>B</sup>*2*<sup>C</sup> , D*∗ 3*A* , *B*2 = *D*<sup>∗</sup> *B*2, 3*A and* 0 ≤ *D*<sup>∗</sup> 3*A* , *B*2 ≤ 2*n* ln 2.

**Definition 9 [15].** *Let* <sup>2</sup>*α<sup>j</sup>* <sup>=</sup> "*aj*, *bj* # , " *cj*, *dj* # ( *j* = 1, 2, ··· , *n*) *be a set of IVIFNs, intervalvalued intuitionistic fuzzy weighted averaging operator is as follows:*

$$\text{InvFlow}\_{\omega}(\overline{a}\_1, \overline{a}\_2, \dots, \overline{a}\_n) = \left( \left[ 1 - \prod\_{j=1}^n (1 - a\_j)^{\omega\_j}, \newline 1 - \prod\_{j=1}^n (1 - b\_j)^{\omega\_j} \right], \left[ \prod\_{j=1}^n c\_j^{\omega\_j}, \newline \prod\_{j=1}^n d\_j^{\omega\_j} \right] \right) \tag{11}$$

*where ω* = (*ω*1, *ω*2, ··· , *ωn*) *Tis the weighting vector of the IVIFNs* <sup>2</sup>*α<sup>j</sup>* (*<sup>j</sup>* <sup>=</sup> 1, 2, ··· , *<sup>n</sup>*).

**Definition 10 [45].** *Given two IVIFNs* <sup>2</sup>*α*<sup>1</sup> <sup>=</sup> ([*a*1, *<sup>b</sup>*1], [*c*1, *<sup>d</sup>*1]) *and* <sup>2</sup>*α*<sup>2</sup> <sup>=</sup> ([*a*2, *<sup>b</sup>*2], [*c*2, *<sup>d</sup>*2])*, the distance of them is as follows:*

$$d(\overline{a}\_1, \ \overline{a}\_2) = \frac{1}{6}(|a\_1 - a\_2| + |b\_1 - b\_2| + |c\_1 - c\_2| + |d\_1 - d\_2| + |e\_1 - e\_2| + |f\_1 - f\_2|) \tag{12}$$

*where <sup>π</sup>*2*α*<sup>1</sup> <sup>=</sup> [*e*1, *<sup>f</sup>*1] *and <sup>π</sup>*2*α*<sup>2</sup> <sup>=</sup> [*e*2, *<sup>f</sup>*2] *are the hesitation degree of* <sup>2</sup>*α*<sup>1</sup> *and* <sup>2</sup>*α*2*, respectively.*

#### *2.2. Three-Way Decision*

Assuming *U* is a finite nonempty set, *R* is an equivalence relation defined on *U*, and *apr*(*α*,*β*) = (*U*, *R*) is a probabilistic rough approximation space, then for *X* ⊆ *U*, let 0 ≤ *β* ≤ *α* ≤1, the upper and lower (*α*, *β*)- approximation sets of *apr*(*α*,*β*) can be expressed as [25]:

$$\begin{cases} apr\_{(a,\emptyset)}(X) = \{ x \in \mathcal{U} | \text{Pr}(X|[x]) \ge a \} \\ \overline{apr}\_{(a,\emptyset)}(X) = \{ x \in \mathcal{U} | \text{Pr}(X|[x]) > \emptyset \} \end{cases} \tag{13}$$

where [*x*] is the equivalence class of *X* with respect to *R*.

In the above formula, Pr(*X*|[*x*]) = |[*x*] ∩ *X*|/|[*x*]| represents the conditional probability of classification, and |·| represents the cardinality of elements in the set. (*α*, *β*)-upper and lower approximation sets divide the domain into three parts, i.e., positive domain *POS*(*α*,*<sup>β</sup>*)(*X*), negative domain *NEG*(*α*,*<sup>β</sup>*)(*X*) and boundary domain *BND*(*α*,*<sup>β</sup>*)(*X*) [27]:


The thresholds *α* and *β* are often given artificially in advance, and so are too subjective and difficult to obtain. Decision rough set introduces Bayesian theory into probability rough set and uses loss function to construct the division strategy of three-way decision with the minimum overall risk, which promotes the development of rough set theory. The decision rough set describes three-way decision processes through the state set Ω = {*X*, ¬*X*} and the action set *A* = {*aP*, *aB*, *aN*}. The state set Ω = {*X*, ¬*X*} represents two states of events, that is, belonging to concept *X* and not belonging to concept *X*. The action set *A* = {*aP*, *aB*, *aN*} indicates that three action strategies of acceptance, delay and rejection can be adopted for different states. Considering that different actions will cause different losses, we record that *λPP*, *λBP* and *λNP*, respectively, represent the losses of actions *aP*, *aB* and *aN* when *x* ∈ *X*, and *λPN*, *λBN* and *λNN*, respectively, represent the losses of actions *aP*, *aB* and *aN* when *x* ∈/ *X*. Therefore, the expected losses of actions *aP*, *aB* and *aN* can be expressed as:

$$\begin{cases} L(a\_P | [x]) = \lambda\_{PP} \Pr(X | [x]) + \lambda\_{PN} \Pr(\neg X | [x]) \\ L(a\_B | [x]) = \lambda\_{BP} \Pr(X | [x]) + \lambda\_{BN} \Pr(\neg X | [x]) \\ L(a\_N | [x]) = \lambda\_{NP} \Pr(X | [x]) + \lambda\_{NN} \Pr(\neg X | [x]) \end{cases} \tag{14}$$

According to Bayesian decision criteria, we select the action set with the minimum expected loss as the best action scheme, and obtain the following three decision criteria [27]:

(**P**): Both *L*(*aP*|[*x*]) ≤ *L*(*aB*|[*x*]) and *L*(*aP*|[*x*]) ≤ *L*(*aN*|[*x*]) are satisfied, then x ∈ POS(X);

(**B**): Both *L*(*aB*|[*x*]) ≤ *L*(*aP*|[*x*]) and *L*(*aB*|[*x*]) ≤ *L*(*aN*|[*x*]) are satisfied, then x ∈ BND(X);

(**N**): Both *L*(*aN*|[*x*]) ≤ *L*(*aP*|[*x*]) and *L*(*aN*|[*x*]) ≤ *L*(*aB*|[*x*]) are satisfied, then x ∈ NEG(X).

Because Pr(*X*|[*x*]) + Pr(¬*X*|[*x*]) = 1, the above rules (**P**)~(**N**) are only related to the conditional probability Pr(*X*|[*x*]) and the loss function *λ*••(• = **P, B, N**). Generally, the loss of accepting the right thing is not greater than that of delaying to accept it, and both of them are less than the loss of rejecting the right thing. The loss of rejecting the wrong thing is not greater than that of delaying rejecting it, and both of them are less than the loss of accepting the wrong thing. Therefore, these loss parameters satisfy the following relationships: 0 ≤ *λPP* ≤ *λBP* < *λNP*, 0 ≤ *λNN* ≤ *λBN* < *λPN*, and the decision rules (**P**)~(**N**) can be rewritten as [27]:

(**P1**): If Pr(*X*|[*x*]) ≥ *α*, *x* ∈ *POS*(*X*); (**B1**): If *β* < Pr(*X*|[*x*]) < *α*, x ∈ BND(X); (**N1**): If Pr(*X*|[*x*]) ≤ *β*, *x* ∈ *NEG*(*X*). where:

$$\begin{cases} \mathfrak{a} = \frac{\lambda\_{PN} - \lambda\_{BN}}{(\lambda\_{PN} - \lambda\_{BN}) + (\lambda\_{BP} - \lambda\_{PP})}\\ \mathcal{B} = \frac{\lambda\_{BN} - \lambda\_{NN}}{(\lambda\_{BN} - \lambda\_{NN}) + (\lambda\_{NP} - \lambda\_{PP})} \end{cases} \tag{15}$$

#### **3. Hybrid Multi-Attribute Three-Way Group Decision Based on Improved VIKOR Model**

Several experts evaluate multiple programs based on multiple indicators. Quantitative indicators may be expressed as exact real numbers, or as interval numbers with minimum and maximum boundaries. Qualitative indicators may be expressed by proper linguistic expressions (values of some linguistic variables), fuzzy numbers, intuitionistic fuzzy numbers or IVIFNs. In accordance with the actual situation, all experts adopt the same expression for the same indicator of each scheme. For this hybrid multi-attribute group decision-making problem, scholars have proposed two different methods. One is to directly construct a hybrid multi-attribute decision matrix and apply TOPSIS, prospect theory, or other methods to make decisions [46,47]. Another is to transform different forms of attributes into the same form and construct a decision model based on a single form of attributes [48–51]. IVIFNs are more flexible and practical in dealing with fuzziness and uncertainty, and other forms of expression can be regarded as special forms of IVIFNs. Therefore, transforming hybrid multi-attribute values into IVIFNs can minimize information distortion. Moreover, after being transformed to the same form, we can effectively calculate the expert weight and attribute weight. Therefore, we choose the latter method for the hybrid multi-attribute group decision-making. The overall decision-making steps are shown in Figure 1.

**Figure 1.** The steps of hybrid multi-attribute three-way group decision making.

#### *3.1. IVIFN Conversion of Different Forms of Attributes*

Let scheme set *G* = {*G*1, *G*2, ... , *Gn*}, attribute set *A* = {*A*1, *A*2, ... , *Am*} and decision maker set *D* = {*D*1, *D*2, ... , *Dl*}. The decision maker *Dk* applies real numbers, interval numbers, values of linguistic variables, intuitionistic fuzzy numbers and IVIFNs to give evaluation value *rij* (*k*) for the attribute *Aj* (*j* = 1, 2, ... , *m*) of the scheme *Gi* (*I* = 1, 2, ... , *n*), thus forming a hybrid multi-attribute decision-making matrix: *R(k) =* [*rij (k)*]*n*×*m*. Where, *rij* (*k*) = *xij* (*k*) is expressed by an exact real number, *rij* (*k*) = [*xij L*(*k*) , *xij R*(*k*) ] by an interval number, *rij* (*k*) = *sij* (*k*) by a linguistic variable value, *rij* (*k*) = (*μij* (*k*) , *vij* (*k*) ) by an intuitionistic fuzzy number, and *r* (*k*) *ij* = *μ*2(*k*) *ij* , *<sup>v</sup>*<sup>2</sup> (*k*) *ij* = *μL*(*k*) *ij* , *<sup>μ</sup>R*(*k*) *ij* , *vL*(*k*) *ij* , *v R*(*k*) *ij* by an IVIFN.

For the intuitionistic fuzzy number (*uij* (*k*) , *vij* (*k*) ), we can transform it to an IVIFN as follows:

$$
\overrightarrow{r}\_{ij}^{(k)} = \left( \left[ \mu\_{ij}^{(k)}, \mu\_{ij}^{(k)} \right], \left[ v\_{ij}^{(k)}, v\_{ij}^{(k)} \right] \right) \tag{16}
$$

For a real number *xij* (*k*) , we first use the linear proportion, vector normalization, extreme value transformation, or other methods to make dimensionless processing. For example, the calculation formula of the linear proportion method is as follows:

$$\mathcal{Y}\_{ij}^{(k)} = \begin{cases} \begin{array}{c} \text{x}\_{ij}^{(k)} \\ \frac{\max\limits\_{k=1,2,\cdots,n} \text{x}\_{kj}^{(k)}}{\max\limits\_{k=1,2,\cdots,n} \text{x}\_{kj}^{(k)} - \text{x}\_{ij}^{(k)} \\ \frac{\max\limits\_{k=1,2,\cdots,n} \text{x}\_{kj}^{(k)} - \text{x}\_{ij}^{(k)} \end{array}}{\max\limits\_{k=1,2,\cdots,n} \text{x}\_{kj}^{(k)} - \text{x}\_{kj}^{(k)} \end{array}, \quad j \in J\_2$$

where *J*<sup>1</sup> is an indicator of benefit type that the larger the better, and *J*<sup>2</sup> is an indicator of cost type that the smaller the better. Then we transform *yij* (*k*) into an intuitionistic fuzzy number *rij* (*k*) = (*yij* (*k*) , 1 − *yij* (*k*) ), and transform *rij* (*k*) into an IVIFN <sup>2</sup>*<sup>r</sup>* (*k*) *ij* <sup>=</sup> *<sup>y</sup>* (*k*) *ij* , *y* (*k*) *ij* , 1 − *y* (*k*) *ij* , 1 − *y* (*k*) *ij* .

For an interval number [*xij L*(*k*) , *xij R*(*k*) ], we first carry out dimensionless processing. For example, the calculation formula of the linear proportion method is as follows:

$$\begin{split} \mathcal{Y}\_{ij}^{L(k)} &= \begin{cases} \frac{\mathbf{x}\_{ij}^{L(k)}}{\max\limits\_{h=1,2,\dots,m} \mathbf{x}\_{hj}^{R(k)}}, & j \in \mathcal{J}\_{1} \\\\ \frac{\max\limits\_{h=1,2,\dots,m} \mathbf{x}\_{hj}^{R(k)} - \mathbf{x}\_{ij}^{R(k)} \\\\ \frac{\max\limits\_{h=1,2,\dots,m} \left(\max\limits\_{h=1,2,\dots,m} \mathbf{x}\_{hj}^{R(k)} - \mathbf{x}\_{jj}^{L(k)}\right)}, & j \in \mathcal{J}\_{2} \end{cases} \\\\ \mathcal{Y}\_{ij}^{R(k)} &= \begin{cases} \frac{\mathbf{x}\_{ij}^{R(k)}}{\max\limits\_{h=1,2,\dots,m} \mathbf{x}\_{hj}^{R(k)}}, & j \in \mathcal{J}\_{1} \\\\ \frac{\max\limits\_{h=1,2,\dots,m} \mathbf{x}\_{hj}^{R(k)} - \mathbf{x}\_{ij}^{L(k)}}{\max\limits\_{h=1,2,\dots,m} \mathbf{x}\_{hj}^{R(k)} - \mathbf{x}\_{ij}^{L(k)}}, & j \in \mathcal{J}\_{2} \end{cases} \end{split} \tag{19}$$

Then we transform [*yij L*(*k*) , 1 − *yij R*(*k*) ] into an IVIFN <sup>2</sup>*<sup>r</sup>* (*k*) *ij* = ([*yL*(*k*) *ij* , *<sup>y</sup>L*(*k*) *ij* ], [1 − *y R*(*k*) *ij* , 1 − *y R*(*k*) *ij* ]).

Let a linguistic evaluation set *S<sup>q</sup>* = *Sq i i* ∈ −*q*−<sup>1</sup> <sup>2</sup> , ··· , <sup>−</sup>1, 0, 1, ··· , *<sup>q</sup>*−<sup>1</sup> 2 , here *<sup>q</sup>* is an odd positive number, the IVIFN corresponding to the *q* linguistic evaluation granularity can be expressed as [52]:

$$\begin{split} \tilde{\mu}^{q} &= \left( \tilde{\mu}^{q} \mathop{\!\!\!}\_{-\frac{q-1}{2}}, \tilde{\mu}^{q} \mathop{\!\!\!}\_{-\frac{q-1}{2}+1}, \cdot \text{ \textquotedblleft}, \tilde{\mu}^{q} \mathop{\!\!\!}\_{0}, \cdots, \"\tilde{\mu}^{q} \mathop{\!\!\!}\_{\frac{q-1}{2}-1}, \tilde{\mu}^{q} \mathop{\!\!\!\!}\_{\frac{q-1}{2}-1} \right) \\ &= \left( \left( \tilde{\mu}^{q}\_{0} \right)^{\frac{q+1}{2}}, \left( \tilde{\mu}^{q}\_{0} \right)^{\frac{q+1}{2}-1}, \cdot \text{ \textquotedblright} \left( \tilde{\mu}^{q}\_{0} \right)^{1/\left( \frac{q+1}{2}-1 \right)}, \left( \tilde{\mu}^{q}\_{0} \right)^{1/\left( \frac{q+1}{2}-1 \right)} \right) \end{split} \tag{20}$$

$$\begin{split} \overline{v}^{\underline{q}} &= \begin{pmatrix} \overline{v}^{\underline{q}} & \overline{v}^{\underline{q}} & \cdots & \overline{v}^{\underline{q}} & \cdots & \overline{v}^{\underline{q}}\_{0} & \cdots & \overline{v}^{\underline{q}}\_{\underline{q}-1} & \overline{v}^{\underline{q}}\\ -\frac{q-1}{T} & \overline{v}^{\underline{q}-1} + 1 & \cdots & \overline{v}^{\underline{q}} & \overline{v}^{\underline{q}} & \cdots & \overline{v}^{\underline{q}} & \cdots & 1 \end{pmatrix} \\ &= \left( 1 - \left( 1 - \overline{v}\_{0}^{\underline{q}} \right)^{\frac{q+1}{2}}, 1 - \left( 1 - \overline{v}\_{0}^{\underline{q}} \right)^{\frac{q+1}{2} - 1}, \cdots, \; \overline{v}\_{0}^{\underline{q}} & \cdots, \; 1 - \left( 1 - \overline{v}\_{0}^{\underline{q}} \right)^{1/\left( \frac{q+1}{2} - 1 \right)}, 1 - \left( 1 - \overline{v}\_{0}^{\underline{q}} \right)^{1/\left( \frac{q+1}{2} \right)} \right) \end{split} \tag{21}$$

where *<sup>μ</sup>*2*<sup>q</sup>* <sup>0</sup> <sup>=</sup> *<sup>v</sup>*<sup>2</sup> *q* <sup>0</sup> = 0.5 − <sup>1</sup> <sup>2</sup>*<sup>q</sup>* , 0.5 . Then, for a linguistic variable value *sij* (*k*) , we determine the linguistic evaluation value of the corresponding level in the *q* granularity, and then express it with the corresponding IVIFN.

In this way, we can transform the hybrid multi-attribute decision-making matrix *R*(*k*) into an interval-valued intuitionistic fuzzy decision matrix *<sup>R</sup>*2*<sup>k</sup>* = 2*r* (*k*) *ij n*×*m* , *k* = 1, 2, ... , *l*, where <sup>2</sup>*<sup>r</sup>* (*k*) *ij* = *μ*2(*k*) *ij* , <sup>2</sup>*<sup>v</sup>* (*k*) *ij* = *μL*(*k*) *ij* , *<sup>μ</sup>R*(*k*) *ij* , *vL*(*k*) *ij* , *v R*(*k*) *ij* .

#### *3.2. Determination of Expert Weight Based on Entropy and Cross Entropy*

In multi-attribute group decision-making, the smaller the difference between the evaluation value of a decision-maker and other decision-makers, the greater weight should be given to this decision-maker. At the same time, the higher the effectiveness of information in a decision-maker's evaluation matrix, that is, the smaller the redundancy, the greater the weight of this decision-maker. In evaluating the redundancy and difference of information, we introduce entropy and cross-entropy to measure them, respectively, and then build a model to determine the weights of experts.

For the evaluation matrix of a single decision maker, we use entropy *E*(*k*) to express the redundancy of evaluation information, and the formula is as follows:

$$E^{(k)} = \frac{1}{m} \sum\_{j=1}^{m} E\_j^{(k)} \tag{22}$$

where *E*(*k*) *<sup>j</sup>* represents the entropy of the *j*th indicator obtained from the decision matrix of the *k*th expert. According to Definition 7, its expression is as follows:

$$E\_j^{(k)} = \frac{1}{n} \sum\_{i=1}^{n} \frac{2 - 2\left(\Delta\_{\tilde{A}}\left(\overline{r}\_{ij}^{(k)}\right)\right)^3 + \left(\pi\_{\tilde{A}}^L\left(\overline{r}\_{ij}^{(k)}\right)\right)^3 + \left(\pi\_{\tilde{A}}^R\left(\overline{r}\_{ij}^{(k)}\right)\right)^3}{4},\tag{23}$$

$$j = 1, 2, \dots, m$$

Based on the entropy of each expert, we can calculate the expert weight as follows:

$$w\_1^{(k)} = \frac{1 - E^{(k)}}{\sum\_{h=1}^{l} \left[1 - E^{(h)}\right]}, \; k = 1, \; 2, \cdots, \; l \tag{24}$$

To reflect the difference between a single decision-making matrix and the other decision-making matrices, we define the cross entropy as follows:

$$D^{(k)} = \frac{1}{(l-1)mn} \sum\_{t=1, t \neq k}^{l} D^\* \left( r^{(k)}, \ r^{(t)} \right) \tag{25}$$

According to Definition 8, the formula of *D*<sup>∗</sup> *r*(*k*),*r*(*t*) is as follows:

$$\begin{split} D^\* \left( r^{(k)}, r^{(t)} \right) \\ &= \sum\_{i=1}^n \sum\_{j=1}^m \left[ d \left( \overline{r}\_{ij}^{(k)} \right) \ln \frac{d \left( \overline{r}\_{ij}^{(k)} \right)}{\frac{1}{2} \left( d \left( \overline{r}\_{ij}^{(k)} \right) + d \left( \overline{r}\_{ij}^{(t)} \right) \right)} \right. \\ &\left. + \left( 1 - d \left( \overline{r}\_{ij}^{(k)} \right) \right) \ln \frac{1 - d \left( \overline{r}\_{ij}^{(k)} \right)}{\frac{1}{2} \left( 2 - d \left( \overline{r}\_{ij}^{(k)} \right) - d \left( \overline{r}\_{ij}^{(t)} \right) \right)} \right. \\ &\left. + d \left( \overline{r}\_{ij}^{(t)} \right) \ln \frac{d \left( \overline{r}\_{ij}^{(t)} \right)}{\frac{1}{2} \left( d \left( \overline{r}\_{ij}^{(k)} \right) + d \left( \overline{r}\_{ij}^{(t)} \right) \right)} \right. \\ &\left. + \left( 1 - d \left( \overline{r}\_{ij}^{(t)} \right) \right) \ln \frac{1 - d \left( \overline{r}\_{ij}^{(t)} \right)}{\frac{1}{2} \left( 2 - d \left( \overline{r}\_{ij}^{(k)} \right) - d \left( \overline{r}\_{ij}^{(t)} \right) \right)} \right] \end{split} \tag{26}$$

Because 0 ≤ *D*<sup>∗</sup> *r*(*k*), *r*(*t*) ≤ 2*mn* ln 2, 0 ≤ *D*(*k*) ≤ 2 ln 2. Then, based on the cross-entropy, we can calculate the expert weight as follows:

$$w\_2^{(k)} = \frac{2\ln 2 - D^{(k)}}{\sum\_{h=1}^{l} \left[2\ln 2 - D^{(h)}\right]}, \; k = 1, \; 2, \cdots, \; l \tag{27}$$

By aggregating *w*(*k*) <sup>1</sup> and *<sup>w</sup>*(*k*) <sup>2</sup> with weight coefficients *γ* and (1-*γ*), respectively, we can calculate the final expert weight as follows:

$$w\_k = \gamma w\_1^{(k)} + (1 - \gamma) w\_2^{(k)} \tag{28}$$

#### *3.3. Determination of Group Comprehensive Evaluation Matrix*

Combined with all the experts' weights, we apply the interval-valued intuitionistic fuzzy weighted averaging operator to calculate the group comprehensive evaluation matrix *X* = " *x*2*ij*# *<sup>n</sup>*×*m*, where:

$$
\begin{split}
\widetilde{\mathbf{x}}\_{ij} &= \left( \left[ \mu\_{ij}^{L}, \mu\_{ij}^{R} \right], \left[ v\_{ij}^{L}, v\_{ij}^{R} \right] \right) = \textit{IIIFWA} \left( \widetilde{r}\_{ij}^{(1)}, \widetilde{r}\_{ij}^{(2)}, \dots, \widetilde{r}\_{ij}^{(l)} \right) \\
&= \left( \left[ 1 - \prod\_{k=1}^{l} \left( 1 - \mu\_{ij}^{L(k)} \right)^{w\_{k}}, 1 - \prod\_{k=1}^{l} \left( 1 - \mu\_{ij}^{R(k)} \right)^{w\_{k}} \right], \left[ \prod\_{k=1}^{l} \left( v\_{ij}^{L(k)} \right)^{w\_{k}}, \prod\_{k=1}^{l} \left( v\_{ij}^{R(k)} \right)^{w\_{k}} \right] \right)
\end{split} \tag{29}
$$

#### *3.4. Determination of Attribute Weight Based on Entropy*

Based on the group comprehensive evaluation matrix, we apply the entropy value method to determine the weight of each attribute:

$$
\omega\_{\dot{j}} = \frac{1 - E\_{\dot{j}}}{\sum\_{\text{lr}=1}^{m} (1 - E\_{\text{lr}})}, \; j = 1, \; 2, \cdots, \; m \tag{30}
$$

where:

$$E\_j = \frac{1}{n} \sum\_{i=1}^n \frac{2 - 2\left(\Delta\_{\bar{A}}\left(\tilde{\mathbf{x}}\_{\bar{i}\bar{j}}\right)\right)^3 + \left(\pi\_{\bar{A}}^L\left(\tilde{\mathbf{x}}\_{\bar{i}\bar{j}}\right)\right)^3 + \left(\pi\_{\bar{A}}^R\left(\tilde{\mathbf{x}}\_{\bar{i}\bar{j}}\right)\right)^3}{4}, \; j = 1, \; 2, \cdots, \; m \tag{31}$$

#### *3.5. Determination of Conditional Probability*

The determination of conditional probability is the key to a three-way decision. The VIKOR method originates from TOPSIS and can take group utility and individual regret into account. Grey correlation analysis can make full use of sample information to reflect the internal law of sample data. We use the VIKOR method improved by grey correlation analysis to determine the conditional probability, and the concrete steps are as follows:

Step 1: According to the evaluation matrix *X*, the positive and negative ideal solutions are as follows:

$$\mathbf{x}^+ = (\overline{\mathbf{x}}\_1^+, \overline{\mathbf{x}}\_2^+, \dots, \overline{\mathbf{x}}\_m^+), \overline{\mathbf{x}}^- = (\overline{\mathbf{x}}\_1^-, \overline{\mathbf{x}}\_2^-, \dots, \overline{\mathbf{x}}\_m^-) \tag{32}$$

where:

$$\begin{cases} \overline{\boldsymbol{x}}\_{j}^{+} = \left( \left[ \max\_{\overline{i}} \mu\_{\overline{i}\overline{j}}^{L}, \max\_{\overline{i}} \mu\_{\overline{i}\overline{j}}^{R} \right], \; \left[ \min\_{\overline{i}} \boldsymbol{v}\_{\overline{i}\overline{j}}^{L}, \min\_{\overline{i}} \boldsymbol{v}\_{\overline{i}\overline{j}}^{R} \right] \right) \\ \overline{\boldsymbol{x}}\_{j}^{-} = \left( \left[ \min\_{\overline{i}} \mu\_{\overline{i}\overline{j}}^{L}, \min\_{\overline{i}} \mu\_{\overline{i}\overline{j}}^{R} \right], \; \left[ \max\_{\overline{i}} \boldsymbol{v}\_{\overline{i}\overline{j}}^{L}, \max\_{\overline{i}} \boldsymbol{v}\_{\overline{i}\overline{j}}^{R} \right] \right) \end{cases} \tag{33}$$

Step 2: Calculate the group utility value *Si* and the individual regret value *Ri* of the *i*th scheme:

$$S\_{i} = \sum\_{j=1}^{m} \frac{\omega\_{j} d\left(\tilde{\boldsymbol{\tilde{x}}}\_{j}^{+}, \tilde{\boldsymbol{\tilde{x}}}\_{ij}\right)}{d\left(\tilde{\boldsymbol{\tilde{x}}}\_{j}^{+}, \tilde{\boldsymbol{\tilde{x}}}\_{j}^{-}\right)}, \ R\_{i} = \max\_{j=1,2,\cdots,m} \frac{\omega\_{j} d\left(\tilde{\boldsymbol{\tilde{x}}}\_{j}^{+}, \tilde{\boldsymbol{\tilde{x}}}\_{i}\right)}{d\left(\tilde{\boldsymbol{\tilde{x}}}\_{j}^{+}, \tilde{\boldsymbol{\tilde{x}}}\_{j}^{-}\right)}, \ i = 1, 2, \text{n} \tag{34}$$

where *d*(*x, y*) represents the distance between two IVIFNs *x* and *y*, which can be calculated according to Definition 10. The smaller the value of *Si*, the higher the group utility. The smaller the value of *Ri*, the smaller the individual regret.

Step 3: Determine the best and the worst group utility values as follows:

$$\mathcal{S}^{+} = \min\_{i=1,2,\cdots,n} \mathcal{S}\_{i\prime} \mathcal{S}^{-} = \max\_{i=1,2,\cdots,n} \mathcal{S}\_{i} \tag{35}$$

The best and the worst individual regret values are:

$$R^{+} = \min\_{i=1,2,\cdots,n} R\_i \quad R^{-} = \max\_{i=1,2,\cdots,n} S\_i \tag{36}$$

Step 4: Calculate the grey correlation degree between the *i*th scheme and the positive and negative ideal solutions as follows:

$$
\varepsilon\_i^+ = \frac{1}{m} \sum\_{j=1}^m \varepsilon\_{ij\prime}^+ \cdot \varepsilon\_i^- = \frac{1}{m} \sum\_{j=1}^m \varepsilon\_{ij\prime}^- \cdot i = 1, 2, \cdots, m \tag{37}
$$

where:

$$\varepsilon\_{ij}^{+} = \frac{\min\_{\mathbf{g} = 1, 2, \dots, n \mathbf{h} = 1, 2, \dots, m} \,\omega\_{h} d\left(\overline{\mathbf{x}}\_{h}^{+}, \overline{\mathbf{x}}\_{\mathbf{g} \mathbf{h}}\right) + \rho \max\_{\mathbf{g} = 1, 2, \dots, n \mathbf{h} = 1, 2, \dots, m} \,\omega\_{h} d\left(\overline{\mathbf{x}}\_{h}^{+}, \overline{\mathbf{x}}\_{\mathbf{g} \mathbf{h}}\right)}{w\_{j} d\left(\overline{\mathbf{x}}\_{j}^{+}, \overline{\mathbf{x}}\_{i \overline{j}}\right) + \rho \max\_{\mathbf{g} = 1, 2, \dots, n \mathbf{h} = 1, 2, \dots, m} \max\_{\mathbf{h}, \mathbf{h} = 1, 2, \dots, m} \omega\_{h} d\left(\overline{\mathbf{x}}\_{h}^{+}, \overline{\mathbf{x}}\_{\mathbf{g} \mathbf{h}}\right)}\tag{38}$$

$$\varepsilon\_{ij}^{-} = \frac{\min\_{\mathcal{S} = 1, 2, \cdots, n \text{h} = 1, 2, \cdots, m} \, \omega\_{\mathcal{h}} d\left(\overline{\boldsymbol{x}}\_{\boldsymbol{h}}^{-}, \overline{\boldsymbol{x}}\_{\mathcal{S}^{\mathcal{H}}}^{\mathcal{A}}\right) + \rho \, \max\_{\mathcal{S} = 1, 2, \cdots, n \text{h} = 1, 2, \cdots, m} \, \omega\_{\mathcal{h}} d\left(\overline{\boldsymbol{x}}\_{\boldsymbol{h}}^{-}, \overline{\boldsymbol{x}}\_{\mathcal{S}^{\mathcal{H}}}^{\mathcal{A}}\right)}{w\_{\mathcal{f}} d\left(\overline{\boldsymbol{x}}\_{\boldsymbol{f}}^{-}, \overline{\boldsymbol{x}}\_{\mathcal{U}}\right) + \rho \, \max\_{\mathcal{S} = 1, 2, \cdots, n \text{h} = 1, 2, \cdots, m} \, \omega\_{\mathcal{h}} d\left(\overline{\boldsymbol{x}}\_{\boldsymbol{h}}^{-}, \overline{\boldsymbol{x}}\_{\mathcal{S}^{\mathcal{H}}}^{\mathcal{A}}\right)} \tag{39}$$

In the above formula, *ρ* ∈ [0, 1] is the distinguishing coefficient. The smaller the value of *ρ*, the greater the distinguishing ability. Generally, *ρ* is taken as 0.5.

Step 5: Calculate the group utility value and individual regret value of the *i*th scheme based on grey correlation analysis as follows:

$$\mathcal{L}\_{i} = \frac{\mathcal{E}\_{i}^{-}}{\mathcal{E}\_{i}^{+}}, \mathcal{J}\_{i} = \max\_{j} \frac{\mathcal{E}\_{ij}^{-}}{\mathcal{E}\_{ij}^{+}}, i = 1, 2, \cdots, n \tag{40}$$

Both the group utility value and the individual regret value are indicators that the smaller the better. Then the best and the worst group utility values are, respectively:

$$\zeta^{+} = \min\_{i=1,2,\cdots,n} \zeta\_{i\prime} \text{ } \zeta^{-} = \max\_{i=1,2,\cdots,n} \zeta\_{i} \tag{41}$$

The best and the worst individual regret values are:

$$\xi^{+} = \min\_{i=1,2,\cdots,n} \xi\_{i\prime}^{\tau} \xi^{-} = \max\_{i=1,2,\cdots,n} \xi\_{i}^{\tau} \tag{42}$$

Step 6: Determine the benefit ratio of the *i*th scheme based on the VIKOR-grey correlation analysis method as follows:

$$Q\_{i} = \sigma \left( \frac{S\_{i}\mathbb{Z}\_{i} - S^{+}\mathbb{Z}^{+}}{S^{-}\mathbb{Z}^{-} - S^{+}\mathbb{Z}^{+}} \right) + (1 - \sigma) \left( \frac{R\_{i}\mathbb{Z}\_{i} - R^{+}\mathbb{Z}^{+}}{R^{-}\mathbb{Z}^{-} - R^{+}\mathbb{Z}^{+}} \right), \ i = 1, 2, \cdots, n \tag{43}$$

where *σ* represents the compromise coefficient between group utility and individual regret, 0 ≤ *σ* ≤ 1. If *σ* > 0.5, it represents the principle of conformity.

Step 7: The smaller the benefit ratio of the *i*th scheme, the greater the probability that it belongs to the acceptable state *Z*. The conditional probability can be calculated as follows:

$$\Pr(Z|\left[G\_{i}\right]) = 1 - Q\_{i} \tag{44}$$

#### *3.6. Determination of Decision Thresholds*

The threshold pair (*α*, *β*) is another key parameter of a three-way decision, which is determined by the loss function. In practice, it is difficult for decision-makers to give the exact value of risk loss of each action under different states. They prefer to use uncertain expressions, such as interval number, fuzzy number, linguistic variable value, intuitionistic fuzzy number and IVIFN. According to the linear or nonlinear ordering rules of various uncertain forms, scholars proposed the corresponding determination methods of the threshold pair [40,41,53,54]. Considering the deficiency of large information distortion in linear ordering, Liu et al. proposed a generalized scalable and nonlinear sorting method to determine the threshold pair for the risk loss matrix represented by IVIFNs from the perspective of optimization [41].

The expert group expresses the risk loss values of three actions *aP* (acceptance), *aB* (delay) and *aN* (rejection) under two states *Z* (acceptable) and *Z<sup>C</sup>* (unacceptable) as IVIFNs, as shown in Table 2.

**Table 2.** Risk loss matrix.


Then the optimization model for solving *α* and *β* is as follows [41]:

*α* = *min g s*.*t*. ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ 2−(*v<sup>L</sup> PZ*) *g vL PZC h* −(*v<sup>R</sup> PZ*) *g vR PZC h* 2+(1−*μ<sup>L</sup> PZ*) *g* 1−*μ<sup>L</sup> PZC h* +(1−*μ<sup>R</sup> PZ*) *g* 1−*μ<sup>R</sup> PZC h* −(*v<sup>L</sup> PZ*) *g vL PZC h* −(*v<sup>R</sup> PZ*) *g vR PZC <sup>h</sup>* ≤ 2−(*v<sup>L</sup> BZ*) *g vL BZC h* −(*v<sup>R</sup> BZ*) *g vR BZC h* 2+(1−*μ<sup>L</sup> BZ*) *g* 1−*μ<sup>L</sup> BZC h* +(1−*μ<sup>R</sup> BZ*) *g* 1−*μ<sup>R</sup> BZC h* −(*v<sup>L</sup> BZ*) *g vL BZC h* −(*v<sup>R</sup> BZ*) *g vR BZC h* 2−(*v<sup>L</sup> PZ*) *g vL PZC h* −(*v<sup>R</sup> PZ*) *g vR PZC h* 2+(1−*μ<sup>L</sup> PZ*) *g* 1−*μ<sup>L</sup> PZC h* +(1−*μ<sup>R</sup> PZ*) *g* 1−*μ<sup>R</sup> PZC h* −(*v<sup>L</sup> PZ*) *g vL PZC h* −(*v<sup>R</sup> PZ*) *g vR PZC <sup>h</sup>* ≤ 2−(*v<sup>L</sup> NZ*) *g vL NZC h* −(*v<sup>R</sup> NZ*) *g vR NZC h* 2+(1−*μ<sup>L</sup> NZ*) *g* 1−*μ<sup>L</sup> NZC h* +(1−*μ<sup>R</sup> NZ*) *g* 1−*μ<sup>R</sup> NZC h* −(*v<sup>L</sup> NZ*) *g vL NZC h* −(*v<sup>R</sup> NZ*) *g vR NZC h g* + *h* = 1, *g*, *h* ≥ 0 (45)

*β* = *max g s*.*t*. ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ 2−(*v<sup>L</sup> NZ*) *g vL NZC h* −(*v<sup>R</sup> NZ*) *g vR NZC h* 2+(1−*μ<sup>L</sup> NZ*) *g* 1−*μ<sup>L</sup> NZC h* +(1−*μ<sup>R</sup> NZ*) *g* 1−*μ<sup>R</sup> NZC h* −(*v<sup>L</sup> NZ*) *g vL NZC h* −(*v<sup>R</sup> NZ*) *g vR NZC <sup>h</sup>* ≤ 2−(*v<sup>L</sup> PZ*) *g vL PZC h* −(*v<sup>R</sup> PZ*) *g vR PZC h* 2+(1−*μ<sup>L</sup> PZ*) *g* 1−*μ<sup>L</sup> PZC h* +(1−*μ<sup>R</sup> PZ*) *g* 1−*μ<sup>R</sup> PZC h* −(*v<sup>L</sup> PZ*) *g vL PZC h* −(*v<sup>R</sup> PZ*) *g vR PZC h* 2−(*v<sup>L</sup> NZ*) *g vL NZC h* −(*v<sup>R</sup> NZ*) *g vR NZC h* 2+(1−*μ<sup>L</sup> NZ*) *g* 1−*μ<sup>L</sup> NZC h* +(1−*μ<sup>R</sup> NZ*) *g* 1−*μ<sup>R</sup> NZC h* −(*v<sup>L</sup> NZ*) *g vL NZC h* −(*v<sup>R</sup> NZ*) *g vR NZC <sup>h</sup>* ≤ 2−(*v<sup>L</sup> BZ*) *g vL BZC h* −(*v<sup>R</sup> BZ*) *g vR BZC h* 2+(1−*μ<sup>L</sup> BZ*) *g* 1−*μ<sup>L</sup> BZC h* +(1−*μ<sup>R</sup> BZ*) *g* 1−*μ<sup>R</sup> BZC h* −(*v<sup>L</sup> BZ*) *g vL BZC h* −(*v<sup>R</sup> BZ*) *g vR BZC h g* + *h* = 1, *g*, *h* ≥ 0 (46)

#### *3.7. Classification and Sorting of Schemes*

According to the value of the threshold (*α*, *β*), we can classify schemes:


In addition, the larger the value of Pr(*Z*|[*Gi*] ), the greater the possibility of selecting the scheme *Gi*. If *α* = *β*, the three-way decision model degenerates into a two-way decisionmaking model. If Pr(*Z*|[*Gi*] ) ≥ *α*, we accept the scheme *Gi*; Otherwise, we reject the scheme *Gi*.

#### **4. An Illustrative Example**

We use the latent dirichlet allocation topic model to mine customers' demand factors for mobile phone performance, and extract six features, namely appearance (*A*1), fast response (*A*2), endurance (*A*3), screen definition (*A*4), running fluency (*A*5) and battery heating (*A*6). We organize four experts *D*1, *D*2, *D*<sup>3</sup> and *D*<sup>4</sup> from China Mobile Communications, China United Network Communications and China Telecommunications in the field of mobile communication performance evaluation to evaluate the above characteristics of the five mobile phone brands *G*1~*G*5. In order to verify the feasibility of the method proposed in this paper, after discussion with experts, the forms of different indicators are set as follows:


The evaluation matrices of the four experts are shown in Tables 3–6, respectively. We will select the brands that can be agented, rejected and pending from the five mobile brands.


**Table 3.** The evaluation matrix of expert *D*1.

**Table 4.** The evaluation matrix of expert *D*2.



**Table 5.** The evaluation matrix of expert *D*3.

**Table 6.** The evaluation matrix of expert *D*4.


We transform all the elements in the above four evaluation matrices into IVIFNs. The results are shown in Tables 7–10.


**Table 7.** The transformed evaluation matrix of expert *D*1.



**Table 9.** The transformed evaluation matrix of expert *D*3.


**Table 10.** The transformed evaluation matrix of expert *D*4.


According to Formulas (22)–(23), we calculate that the values of entropy *E*(1), *E*(2), *E*(3) and *E*(4) are 0.341134, 0.360570, 0.331364 and 0.336861, respectively. Then, according to (24), the four experts' weights *w*<sup>1</sup> (1), *w*<sup>1</sup> (2), *w*<sup>1</sup> (3) and *w*<sup>1</sup> (4) are 0.250513, 0.243123, 0.254227 and 0.252137, respectively. According to (25)–(26), the values of cross entropy *D*(1), *D*(2), *D*(3) and *D*(4) are 0.03027, 0.033022, 0.035811 and 0.037154, respectively, and according to (27), the expert' weights *w*<sup>2</sup> (1), *w*<sup>2</sup> (2), *w*<sup>2</sup> (3) and *w*<sup>2</sup> (4) are 0.222170, 0.242347, 0.262814 and 0.272669, respectively. Taking the weight coefficient *γ* as 0.5 and substituting it with (28), the final expert weights *w*1, *w*2, *w*<sup>3</sup> and *w*<sup>4</sup> are 0.236341, 0.242735, 0.258521 and 0.262403, respectively. Combined with expert weights, we apply the formula (29) to calculate the group comprehensive evaluation matrix, as shown in Table 11.


**Table 11.** The group comprehensive evaluation matrix.

According to (31), we calculate that the entropy values of six attributes are 0.105597, 0.187031, 0.439343, 0.378227, 0.423710 and 0.469141, respectively. Then, according to (30), the six attributes' weights are 0.223771, 0.203397, 0.140271, 0.155562, 0.144183 and 0.132816, respectively. According to (32)–(33), we obtain the positive and negative ideal solutions are:

*<sup>x</sup>*2<sup>+</sup> = (([1, 1], [0, 0]), ([0.9364, 0.9364], [0, 0]), ([0.682, 0.731], [0.2298, 0.269]), ([0.8315, 0.8315], [0.0783, 0.0783]), ([0.7313, 0.8067], [0.1003, 0.1309]), ([0.7305, 0.7722], [0.1934, 0.2278])) *<sup>x</sup>*2<sup>−</sup> = (([0.9079, 0.9079], [0.0921, 0.0921]), ([0.8293, 0.8293], [0.1069, 0.1069]), ([0.1329, 0.1901], [0.7506, 0.8099]), ([0.4602, 0.4602], [0.4589, 0.4589]), ([0.1329, 0.1901], [0.7506, 0.8099]), ([0.4602, 0.4602], [0.4589, 0.4589]), ([0.5418, 0.6303], [0.2308, 0.2731]), ([0.4013, 0.4675], [0.4709, 0.5325]))

> According to (34)–(44), We calculate the group utility value *Si*, the individual regret value *Ri*, the group utility value *ζ<sup>i</sup>* of grey correlation analysis, the individual regret value *ξ<sup>i</sup>* of grey correlation analysis, the benefit ratio *Qi* and the conditional probability Pr(*Gi*) of each mobile phone brand in turn. The results are shown in Table 12.

**Table 12.** The conditional probability value of each scheme based on improved VIKOR model.


The four experts jointly give the risk loss matrix represented by IVIFNs, as shown in Table 13.

**Table 13.** Risk loss matrix results.


We substitute the data in the above table into the nonlinear programming models (45) and (46) and obtain that *α* = 0.608646 and *β* = 0.122339. It can be seen that the conditional probabilities of *G*<sup>3</sup> and *G*<sup>4</sup> are greater than *α*, indicating that the two mobile phone brands can be chosen as an agent. If the conditional probability of *G*<sup>1</sup> is less than *β*, this mobile phone brand shall be excluded. The conditional probabilities of *G*<sup>2</sup> and *G*<sup>5</sup> are between *α* and *β*, so they need to be further investigated

In order to reflect the difference between the improved VIKOR model and other conditional probability models, we calculated the conditional probability results and the three classification results under TOPSIS, grey correlation analysis and traditional VIKOR models. The results are shown in Table 14.

**Table 14.** The results under TOPSIS, grey correlation analysis and VIKOR models.


It can be seen that the conditional probability results of grey correlation analysis are too close to effectively distinguish the differences between brands. TOPSIS results of different brands are different to some extent, but brands G1, G2 and G3 are all pending, which indicates that the distinction is not obvious enough. Of course, this is related to the risk loss matrix given by decision-makers. However, considering only the proximity to positive and negative ideal points, it is difficult to reflect the intrinsic characteristics of data. Nor does it capture decision-makers' attitudes to utility and regret. The results of the VIKOR method are similar to those of improved VIKOR, but there are differences in brand G1, which is greatly related to the addition of grey correlation analysis results reflecting the inherent characteristics of data. In general, the improved VIKOR model can not only reflect the proximity to the ideal points, but also reflect the inherent characteristics of data and decision-makers' trade-offs on utility and regret, and the results of it are relatively objective.

#### **5. Conclusions**

For the hybrid multi-attribute decision-making problem, we propose a three-way group decision-making method based on the improved VIKOR model. Based on the transformed interval-valued intuitionistic fuzzy decision matrix, we apply entropy and cross-entropy to determine the expert weights and obtain the group comprehensive evaluation matrix. Then, we use entropy to obtain attribute weights. By using the improved VIKOR method by grey correlation analysis, we determine conditional probability. By

comparing the conditional probability with the decision threshold pair based on optimization, we obtain the classification rules of the three-way decision. The example analysis shows that the method has good three-way classification and can provide support for actual management decision-making. This study has the following features and benefits: First, it considers the hybrid multi-attribute environment, especially the interval-valued intuitionistic fuzzy environment containing more fuzzy information, which is closer to the actual decision-making and has better universality. Second, considering the group decision-making environment, the hybrid multi-attribute evaluation matrix is given by each expert, which is more consistent with reality. Moreover, the proposed expert weight determination method can not only reflect the differences among experts' opinions, but also reflect the uncertainty degree of each expert's evaluation opinion, and the obtained weights are more reasonable and objective. Different from scholars' studies, this paper mainly has three aspects of innovation. First, from the perspective of research, it expands the research of hybrid multi-attribute decision-making and three-way group decision-making. Second, it deepens the research on expert weights and attribute weights in interval-valued intuitionistic fuzzy group decision making and improves the objectivity of weights. Thirdly, an improved VIKOR model based on grey correlation analysis is proposed to determine the conditional probability, which improves the scientificity of the conditional probability.

There are some shortcomings in this study. First, in the determination of expert weights and attribute weights, only one form of interval-valued intuitionistic fuzzy entropy is considered. In fact, there are many forms of interval-valued intuitionistic fuzzy entropy that meet the axiom conditions. How do they affect the weight results and final results, and whether there will be contradictory conclusions? These are not tested. Second, for the risk loss matrix represented by IVIFNs, we use the threshold determination method based on the optimization model, but there is another interactive threshold determination method, that is, to determine the losses based on the preference coefficient and the distance from the ideal points, and then calculate the thresholds. How much is the difference between the results of these two methods? In addition, is it more advantageous to combine the two, that is, to first determine the threshold loss matrix in an interactive way and then determine the thresholds by an optimization method? These aspects are also not explored. Third, we adopt the method of conditional probability determination of the improved VIKOR model. In fact, the prospect theory based on an ordinary utility curve is being gradually introduced to determine conditional probability. Limited by the fact that the prospect theory based on the IVIFN decision matrix is not perfect, we have not conducted research on this aspect. Based on the shortcomings of the method, further research can be conducted in the following aspects. First, we can analyze the influence of other forms of interval-valued intuitionistic fuzzy entropy on expert weights, attribute weights, and the final results. Second, based on the risk loss matrix expressed by IVIFNs, we can discuss the impact of other threshold determination methods on the decision results. Third, we can further improve the prospect theory based on the IVIFN decision matrix and introduce it into the determination of the conditional probability of a three-way decision.

**Author Contributions:** Conceptualization, J.S.; methodology, J.S.; writing—original draft preparation, Z.H.; investigation, L.J.; writing—review and editing, Z.L.; supervision, X.L.; project administration, Z.L. and L.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China. (NO.71871222).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We greatly appreciate the associate editor and the anonymous reviewers for their insightful comments and constructive suggestions, which have greatly helped us to improve the manuscript and guide us forward to the future research.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Jun Wu 1, Xin Liu 1, Yuanyuan Li 1, Liping Yang 1, Wenyan Yuan <sup>2</sup> and Yile Ba 1,\***


**Abstract:** Distribution centers are quite important for logistics. In order to save costs, reduce energy consumption and deal with increasingly uncertain demand, it is necessary for distribution centers to select the location strategically. In this paper, a two-stage model based on an improved clustering algorithm and the center-of-gravity method is proposed to deal with the multi-facility location problem arising from a real-world case. First, a distance function used in clustering is redefined to include both the spatial indicator and the socio-economic indicator. Then, an improved clustering algorithm is used to determine the optimal number of distribution centers needed and the coverage of each center. Third, the center-of-gravity method is used to determine the final location of each center. Finally, the improved method is compared with the traditional clustering method by testing data from 12 cities in Inner Mongolia Autonomous Region in China. The comparison result proves the proposed method's effectiveness.

**Keywords:** multi-facility location problem; clustering algorithm; center-of-gravity method

**MSC:** 68T20

#### **1. Introduction**

Selecting a proper site for a distribution center can effectively save costs, increase profits, improve customer satisfaction and reduce circulation time. Therefore, the location problem has been one of the most important decisions in logistics system planning [1]. Among these, the selection of distribution center sites has become one of the most popular research interests in logistics management. Readers interested in this field could refer to some recent research [2–7].

As the logistics in real-world become increasingly complicated due to factors that significantly increase the uncertainty in logistics planning, such as economy development, the prevalence of online retailing, increasing natural disasters, etc. It is necessary to introduce more indicators to deal with the uncertainty during the site selection process. This research includes demand uncertainty in the model described in a proper mathematical function, which takes both the spatial condition and the economic condition into consideration when selecting the proper locations for distribution centers. Hopefully, this research could provide a planning strategy for the selection of distribution center sites.

In this paper, all data is extracted from a real case of 12 cities in Inner Mongolia Autonomous Region, China. The author proposes a two-stage model based on an improved clustering algorithm and the center-of-gravity method to solve the location problem of distribution centers and provides a strategy for the logistics transportation in this case. The improved model proposed in this paper is compared with the traditional clustering method

**Citation:** Wu, J.; Liu, X.; Li, Y.; Yang, L.; Yuan, W.; Ba, Y. A Two-Stage Model with an Improved Clustering Algorithm for a Distribution Center Location Problem under Uncertainty. *Mathematics* **2022**, *10*, 2519. https:// doi.org/10.3390/math10142519

Academic Editors: Stelios Psarakis and Andreas C. Georgiou

Received: 10 May 2022 Accepted: 15 July 2022 Published: 20 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

by testing data from 12 cities and is proved to be effective. The main work of this paper is as follows:


The remainder of this paper is organized as follows. Section 2 reviews the literature on location problems. Section 3 presents the methodology of a two-stage model, including the improved clustering algorithm and the center-of-gravity method. Section 4 describes the obtained results using the proposed model. Section 5 compares the results between proposed model and different clustering algorithms. Section 6 concludes the research.

#### **2. Literature Review**

The location problem (LP) is a classical problem. The original problem can be traced back to the Fermat problem. The first industrial application of the Fermat problem was proposed by Alfred Weber in 1909 and called the Weber problem [8]. The purpose of the Weber problem is to find the location of a warehouse to minimize the total distance among all customers. Ever since Weber's seminal work, LP has been booming in both the academical world and the industrial world. In previous researches, LP is well classified on the basis of model assumptions, constraints and objectives functions [9].

Based on the number of facilities involved, the location problem could also be classified into two categories, the single-facility location problem (SFLP) or the multi-facility location problem (MFLP). The SFLP is to determine the location of a single facility using the criteria of distance and cost. The common methods include the Weiszfeld method [10] and the center-of-gravity method [11]. However, one single facility might not be able to meet the increased service demand for its capacity is limited. This is where it becomes a multi-facility location problem, which is the problem to be studied in this paper.

The MFLP is commonly seen in humanitarian relief [12,13], emergency response [14], health care facility locations [15,16] and other fields. The objective function and criteria for facility location are considered from different perspectives in many studies, such as total cost, total transportation time and satisfaction [17]. Because MFLP covers a wide range and need to consider many factors, it can be divided into different problems, for example, it could be divided as the deterministic problem if the information in the problem is determined, and the probabilistic and stochastic problem if not. The former mainly focuses on coverage maximization [18,19] and resource minimization [20]. The latter had been largely discussed in earlier researches on the uncertainties of location problems, such as uncertain demand [21,22] and disruption risk [23]. Probabilistic and stochastic location problems can be divided into stochastic programming problems [24] and robust optimization problems [25]. For the detailed classification of MFLPs, interested readers can refer to Wang et al. [26].

There are two commonly methodologies used in MFLP researches: exact methods and heuristic algorithms [27]. In early researches, the scale of the location problem is usually small. Branch-and-bound [28] and cutting plane [29] are quite popular. However, with the expansion of the problem scales and the increase of the constraints, it becomes difficult to obtain the accurate solution. Finding a proper solution of a large-scale problem in a reasonable time has become the mainstream of LP researches, among which, the heuristic algorithm has been proposed for this purpose. Commonly used heuristic algorithms include genetic algorithm [30], particle swarm optimization [31], tabu search [32], cuckoo search [33] and Lagrangian relaxation [34]. However, the most distinctive disadvantage of heuristic algorithms is the computationally burdensome to find the best solution.

How to deal with the computational burden has then become many researchers' interests. One widely adopted perspective is the clustering-based algorithms. By defining the location problem as a kind of clustering problem, the researcher would divide customers or demand points into different clusters, each of which is supplied by one facility. Then there would be only one SFLP in each cluster needs to be solved. Thus, by decomposing the multifacility location problem (MFLP) into multiple single-facility location problems (SFLP), the researcher could simplify the complexity of calculation [35]. Among those researchers, Esnaf and Kucukdeniz [36] proposed a hybrid method using the fuzzy clustering algorithm to divide the whole region into certain number of areas, each supplied by one very facility, then they used the center-of-gravity method to select the location point of each facility. Kuecuekdeniz et al. [37] proposed a method which integrates fuzzy c-means and convex programming for solving a capacitated MFLP, where fuzzy c-means is used to convert an MFLP to an SFLP. Gupta et al. [38] used fuzzy c-means and particle swarm optimization to optimize the locations of public service centers. Researches mentioned above proved that the clustering-based algorithms are applicable in practice [39]. Nonetheless, few researchers have dedicated to improving the distance function in clustering.

In this paper, an MFLP of distribution centers arising from a real-world case is studied with clustering-based algorithms. The author uses the clustering algorithm to decompose the MFLP into multiple SFLPs and then uses the center-of-gravity method to select the locations of distribution centers by solving each SFLP. The work differs clearly from the previous literature. In this paper, both spatial and socio-economic indicators are well embedded in the clustering method. An improved distance function is proposed to explain the clustering results affected by the two types of indicators. Furthermore, the effectiveness of the proposed method is proved by a comparison on different methods.

#### **3. Methodology**

The proposed two-stage model (2SM) aims to compute the optimal number and locations of distribution centers with minimal total transportation cost. The structure of the model is composed of two sequential stages.

1. Determine the optimal number of distribution centers by improved clustering for an MFLP.

In this stage, the author proposes an improved *K*-Medoids clustering algorithm to classify the demand points and determine the optimal clustering strategy. The improved clustering algorithm can be summarized as follows:


After clustering, the location of each cluster distribution center is solved as a SFLP by employing the center-of-gravity method, which optimizes transportation cost. The framework of the model is given in Figure 1.

**Figure 1.** Framework of the proposed model.

#### *3.1. Improved Clustering Algorithm*

The clustering algorithm [40] is a kind of unsupervised machine learning method. Its main idea is to divide the objects into different clusters to maximize the object similarity in each cluster and minimize the object similarity between any two clusters. To deal with noise data and outliers effectively, the author employs the *K*-Medoids clustering algorithm and modifies it to fit the work. The *K*-Medoids clustering algorithm [41,42] is an adjustment of the classical *K*-Means clustering algorithm. The pseudocodes for the *K*-Medoids are shown in Algorithm 1.

**Algorithm 1** *K*-Medoids Clustering Algorithm

**Input**: Data points {*xi*}*<sup>N</sup> <sup>i</sup>*=1, number of clusters *K* Choose different *xi* as initial clustering center *ck* for *k* = 1 to *K* Set *Ck* ← {*ck*} for *k* = 1 to *K* **for** *i* = 1 to *K* **do** *Cp* ← *Cp* ∪ {*xi*} , where p = argmin*p*∈*I<sup>K</sup>* 1 *d*(*xi*, *cP*) **end while** allocation result of any *xi* changed **do for** *k* = 1 to *K ck* ← argmin*xi*∈*Ck* ∑ *xt*∈*Ck d*(*xt*, *xi*) **end for** *i* = 1 to *N* **do if** *xi* ∈ *Ck* and *d*(*xi*, *ck*) < *d*(*xi*, *ck*) **then** *Ck* ← *Ck* ∪ {*xi*} and *Ck* ← *Ck* \{*xi*} **end end end Output**: Clustering results {*Ck*}*<sup>K</sup> k*=1

Different from the *K*-Means, in the iterative process, the *K*-Medoids selects the data point closest to the data's mean in the cluster as the new clustering center instead of choosing the average of all data points.

In the clustering algorithms, there are three factors affect the performance of clustering directly: distance measurement, algorithm efficiency and number of clusters. In consideration of the factors, an improved clustering method is proposed to meet the real-world situation.

#### 3.1.1. Redefine the Distance Function

Euclidean distance is a common distance function in clustering algorithm that measures the distance between two points in *m*-dimensional space. The two-dimensional distance function considering the longitude and the latitude is as follows:

$$D\_{ij} = \sqrt{\left(\mathbf{x}\_i - \mathbf{x}\_j\right)^2 - \left(y\_i - y\_j\right)^2} \tag{1}$$

where *Dij* measures the distance between point *i* and point *j*; *xi* and *xj* are the longitudes of point *i* and point *j*, respectively; and *yi* and *yj* are the latitudes of point *i* and point *j*, respectively.

Euclidean distance is simple and effective when dealing with homogeneous indicators. But when the indicators are different, it misses some practical interpretability. To cope with this setback, the author redefines the distance function considering both spatial indicators and socio-economic indicators. The spatial indicators represent the actual distance between two points. The socio-economic indicators reflect the logistics-level score of each point. The detailed computing process is as follows.

#### (1) Compute the logistics-level score.

The logistics-level score proposed in this research mainly considers three dimensions: economy development, traffic congestion and total logistics demand, which is represented respectively by the per capita disposable income of urban residents, the population density and the permanent urban population, as shown in Table 1. Note that there is no available data on population density. Instead, the ratio of permanent urban residents to urban built-up area is used as in Equation (2).

$$population\ density = \frac{permanent\,\,urbam\,\,population}{unbank\,\,bulk\,\,up\,\,area} \tag{2}$$

**Table 1.** The indicator composition of the logistics-level score.


As the actual data of the corresponding indicators is obtained, the author computes each city's logistics-level score using the entropy weight method, which is used to determine the weight coefficient of each indicator by computing the information entropy. To a certain extent, the entropy weight method can avoid subjective judgment when determining weight coefficients [43]. The smaller the information entropy value of one indicator is, the larger the weight coefficient assigned to it becomes [44]. The logistics-level score is computed as Equations (3)–(8):

First, the data are normalized by Equations (3) and (4):

$$Y\_{ij} = \frac{X\_{ij} - \min\_{i} \{ X\_{ij} \}}{\max\_{i} \{ X\_{ij} \} - \min\_{i} \{ X\_{ij} \}} \tag{3}$$

$$Y\_{ij} = \frac{\max\_{i} \left( X\_{ij} \right) - X\_{ij}}{\max\_{i} \left( X\_{ij} \right) - \min\_{i} \left( X\_{ij} \right)} \tag{4}$$

Equation (3) is the normalized formula of the utility indicator, and Equation (4) is the normalized formula of the cost indicator. *Xij* is the *j*th (*j* = 1, 2 ... , *m*) indicator of city *i* (*i* = 1, 2 ... , *n*). *Yij* is the indicator data after normalization. Since the greater the population density is, the greater the traffic congestion becomes, population density should be considered as the cost indicator. The other two indicators are utility indicators.

Then, the information entropy is computed by Equations (5) and (6):

$$P\_{ij} = \mathbf{Y}\_{ij} / \sum\_{i=1}^{n} \mathbf{Y}\_{ij} \tag{5}$$

$$E\_{\bar{l}} = -\frac{1}{\text{Im}} \sum\_{i=1}^{n} P\_{\bar{i}\bar{j}} \ln P\_{\bar{i}\bar{j}} \tag{6}$$

where *Pij* is the weight of city *i* in the *j*th indicator and *Ej* is the information entropy of the *j*th indicator.

After that, the weights of the indicators are computed by Equation (7):

$$\mathcal{W}\_{\dot{\jmath}} = \begin{pmatrix} 1 - E\_{\dot{\jmath}} \end{pmatrix} / \left( n - \sum\_{j=1}^{m} E\_{\dot{\jmath}} \right) \tag{7}$$

where *Wj* is weight of the *j*th indicator.

Finally, the logistics-level score of each city is computed by Equation (8):

$$Z\_i = \sum\_{j=1}^{m} Y\_{ij} \mathcal{W}\_j \tag{8}$$

where *Zi* is the logistics-level score of city *i*.

(2) Define the distance between two points:

$$D\left(\mathbf{X}\_{i\prime}\mathbf{X}\_{j}\right) = \frac{D\_{ij}^{2}}{\left(Z\_{i}Z\_{j}\right)^{u}}\tag{9}$$

In Equation (9), *D Xi*, *Xj* measures the distance between any two cities, which is affected by two distance indicators: spatial and socio-economic. *D*<sup>2</sup> *ij* is the spatial distance between any two cities, which represents the spatial indicators. The greater the spatial distance is, the greater the distance between the two cities becomes. *ZiZj* is the logistics score distance between any two cities, which represents the socio-economic indicators. *Zi* and *Zj* are the logistics-level scores of city *i* and city *j*, respectively. The greater the logistics score distance is, the closer the connection between the two cities at the logistics level lies, and the smaller the distance between the two cities in the clustering process becomes. *u* represents the degree of influence of the logistics-level score in the distance function between any two cities. The larger *u* is, the greater the influence of the logistics-level score becomes. When *u* is 0, the distance function is equivalent to the Euclidean distance.

#### 3.1.2. Select the Initial Clustering Center

To reduce the time complexity of the clustering algorithm effectively and avoid the interference of noise data, the author introduces the idea of density into the clustering algorithm to determine the initial clustering center in this research. The steps of selecting the initial clustering center are as follows.

First, the domain radius *ra* of the demand points *Xi* is computed by Equation (10):

$$r\_a = \frac{1}{2} \max\{ \|X\_i - X\_k\| \}\tag{10}$$

The domain radius *ra* measures the maximum distance between *Xi* and other demand points.

Then, the density factor *Di* of *Xi* is computed by Equation (11):

$$D\_i = \sum\_{j=1}^{n} \exp\left[ -\frac{||X\_i - X\_j||^2}{\left(2r\_a\right)^2} \right] \tag{11}$$

A larger *Di* means that more demand points surround *Xi*, i.e., *Xi* has a closer relationship with other points. The point with a greater density factor is more likely to be a clustering center. The exponential function could make sure that the point surrounded by more demand has a greater density factor.

After figuring out the density factors of all demand points, the factors are arranged in the descend order. Then, the demand points are selected corresponding to the top *k* density factors as the initial clustering center.

#### 3.1.3. Determine the Optimal *k*

The clustering algorithm aims to divide the similar sample points into the same cluster and effectively distinguish the dissimilar sample points. In this process, *k* will significantly affect the final clustering performance. To obtain the optimal clustering results, it is necessary to determine the optimal selection of *k* before clustering. In this research, the elbow method is used to obtain the optimal *k* by computing the sum of the squares of error (*SSE*) between the sample points and their respective clustering centers. The equation for computing *SSE* is Equation (12):

$$SSE = \sum\_{i=1}^{k} \sum\_{\mathbf{x} \in \mathcal{C}\_i} (\mathbf{x} - P\_i)^2 \tag{12}$$

where *Ci* represents the *i*th cluster, *Pi* is the clustering center of *Ci* and *x* is the point belonging to the *Ci*.

The core idea of the elbow method lies in the following: With the increase in cluster number *k*, sample division will be more refined, and the aggregation degree of each cluster will be gradually improved; thus, *SSE* will gradually become smaller. In addition, when *k* is less than the optimal cluster number, *SSE* decreases greatly because the increase of *k* will greatly increase the degree of aggregation of each cluster. However, when *k* reaches the optimal cluster number, the return of aggregation degree obtained by increasing *k* decreases rapidly. This means decrease of *SSE* will tend to be gentler as *k* continues to increase.

#### 3.1.4. The Flowchart of Clustering Algorithm

Based on the improvements discussed above, the clustering algorithm flow is as follows:

Step 1: Obtain the data of each city required by the model, normalize the data for each indicator to eliminate the dimensional influence.

Step 2: Compute the density factor of each city and arrange each city in the order of density factor from large to small.

Step 3: Take the top *k* points in Step 2 as clustering centers and take the remaining *n* − *k* points as sample points for clustering. Corresponding clustering results are obtained.

Step 4: Compute the sum of distances between each point in the cluster and the remaining points in the same cluster, redetermine the center of each cluster according to the principle of the minimum sum of distances. On this basis, obtain the results by clustering again.

Step 5: Repeat Step 4 until clustering result obtained becomes stable. Take this result as the result under *k*.

Step 6: Change the value of *<sup>k</sup>* (*<sup>k</sup>* <sup>≤</sup> <sup>√</sup>*n*) and repeat Steps 3–5 to obtain clustering results under different *k*.

Step 7: Use the elbow method to obtain the optimal *k* and the optimal clustering scheme. Step 8: Output the optimal clustering results.

Figure 2 shows the flow chart of the clustering algorithm.

**Figure 2.** The flow chart of the clustering algorithm.

#### *3.2. The Center-of-Gravity Method*

The center-of-gravity method is used to obtain the optimal distribution center location within each cluster after clustering. The point is to minimize the transportation cost between any demand point and the distribution center.

In the logistics distribution system, the coordinate of each logistics demand point is defined as *i*(*xi*, *yi*)(*i* = 1, 2, . . . ., *n*) according to geographical location, and the coordinate of distribution center P is set as (*x*0, *y*0). Equation (13) shows the objective function:

$$H = \sum\_{i=1}^{n} h\_i w\_i d\_i \tag{13}$$

where *H* is total transportation cost, *hi* is demand of goods at point *i*, *wi* is transportation cost from distribution center to point *i* and *di* is distance of distribution center to point *i*.

Here *hi* is defined as follows:

$$h\_i = \sum\_{s \in S} p\_{is} a\_{is} u\_i = u\_i \sum\_{s \in S} p\_{is} a\_{is} \tag{14}$$

where *S* is the set of scenarios, *pis* is the probability of scenario *s* occurring at point *i*, *ais* is the demand coefficient when scenario *s* occurs at point *i*, *ais* = 1 denotes a stable demand case and *ui* is per capita disposable income of urban residents at point *i*.

In the first stage of this method, the gravity centers of each cluster are computed by the following formulas:

$$\log\_0 = \frac{\sum\_{i=1}^n h\_i w\_i \mathbf{x}\_i / d\_i}{\sum\_{i=1}^n h\_i w\_i / d\_i} \tag{15}$$

$$y\_0 = \frac{\sum\_{i=1}^{n} h\_i w\_i y\_i / d\_i}{\sum\_{i=1}^{n} h\_i w\_i / d\_i} \tag{16}$$

$$d\_i = \sqrt{\left(\mathbf{x}\_0 - \mathbf{x}\_i\right)^2 + \left(y\_0 - y\_i\right)^2} \tag{17}$$

P(*x*0, *y*0) is the extreme point, which is a necessary condition for the optimal solution. If it is not optimal, the iterative method is introduced to find the optimal solution. The solution formulas are as follows:

$$\mathbf{x}\_0^{(q+1)} = \frac{\sum\_{i=1}^n h\_i w\_i \mathbf{x}\_i / d\_i^{(q+1)}}{\sum\_{i=1}^n h\_i w\_i / d\_i^{(q+1)}} \tag{18}$$

$$\mathcal{Y}\_0^{(q+1)} = \frac{\sum\_{i=1}^n h\_i w\_i y\_i / d\_i^{(q+1)}}{\sum\_{i=1}^n h\_i w\_i / d\_i^{(q+1)}} \tag{19}$$

$$d\_i^{(q+1)} = \sqrt{\left(\mathbf{x}\_0^{(q+1)} - \mathbf{x}\_i\right)^2 + \left(y\_0^{(q+1)} - y\_i\right)^2} \tag{20}$$

This computing process ends when the difference between the last two coordinates of distribution center is lower than a specific minimal value.

#### *3.3. Evaluation Function*

In this paper, the evaluation function's object is evaluating the sum of the transportation costs from each distribution center to its responsible demand points. To describe the function accurately, assumptions should be made and given as follows. The capacity of each distribution center can always meet all the demands of the points covered by it. Each demand point is supplied by only one distribution center. Omit the rental, maintenance and other costs related to transporting vehicles during the transportation process. Transportation costs include only freight rates and transportation distance, with the exclusion of labor costs incurred from loading and unloading. Based on the assumptions above, the evaluation function can be calculated as follows:

$$F = \sum\_{i=1}^{k} \sum\_{j \in N\_i} m\_{ij} w\_{ij} f \tag{21}$$

where *i* indicates the distribution center, *j* indicates the demand point, *Ni* is the set of all requirement points covered by *i* distribution center, *mij* is the transportation distance from the distribution center *i* to the demand point *j*, *wij* is the total weight of goods transported from distribution center *i* to the demand point *j*, and *J* is the unit rate from distribution center *i* to the demand point *j*.

#### **4. Case Study**

In this section, taking Inner Mongolia Autonomous Region as an example, the author analyzes the locations of distribution centers in 12 cities to provide a better solution for its logistics transportation. The algorithm is implemented in Python version 3.8.

#### *4.1. Data Collection and Processing*

There are six indicators involved: the longitude, the latitude, the per capita disposable income of urban residents, the permanent urban population, the urban built-up area and the population density.

The longitude and the latitude are obtained from Baidu Maps. The per capita disposable income of urban residents, the permanent population and the urban built-up area of each city in 2019 are downloaded from the official website of the Inner Mongolia Bureau of Statistics. The population density of each city is computed by Equation (2). Table 2 shows the initial data used in this case.

**Table 2.** The initial data.


<sup>1</sup> PCDIUR: Per capita disposable income of urban residents. <sup>2</sup> PUP: Permanent urban population.

In order to eliminate the influence of data dimension on the results, all indicators are normalized in advance.

#### *4.2. The Location Problem's Solution*

According to the methodology proposed in Section 3, the location problem is solved as follows.

Compute the domain radius *ra* and the density factor *Di* of each city using the latitude and the longitude. Sort the density factor in the descend order, as shown in Table 3.

**Table 3.** Domain radius *ra* and density factor *Di* for each city.


Select the top *k* cities in Table 3 as the initial centers for clustering. Use the elbow method to evaluate the clustering results with different *k*.

Since it contains 12 demand points in total, the value range of *k* lies in [2, 4] according to the principle of *<sup>k</sup>* <sup>≤</sup> <sup>√</sup>*n*. Figure <sup>3</sup> shows the sum of the squares of error (*SSE*) when *<sup>k</sup>* takes the values in [2, 6]. The optimal value is *k* = 3.

**Figure 3.** The *SSE* when *k* lies in [2, 6].

After selecting the optimal *k*, clustering results are obtained according to the improved clustering algorithm, which is shown in Table 4.



The result is modified by the center-of-gravity method. In this process, the permanent urban population is taken to represent the demand of each city. The modified results are shown in Table 5 and Figure 4. C0, C1, C2 and Center are in different colors in Figure 4.



**Figure 4.** The location results for the 2SM.

A distribution center cannot be set up outside the city because of its infrastructure requirements. So, the distribution center is sited in a city closest to the location results for the 2SM. The results are shown in Table 6.

**Table 6.** The final location results for the 2SM.


#### **5. Evaluation**

To evaluate the effectiveness of the two-stage model (2SM), the author makes a comparison on the results of it with three models' based on the traditional *K*-means algorithm. For each model, *k* is set as 3. Details are shown in the following.

Model 1: Geographic clustering model (GCM). This model considers only the geographical indicator. It takes the linear distance between cities as the similarity measurement. The distance between cities is computed based on the longitude and the latitude of each city. The location results are obtained by *K*-means clustering. The results are shown in Table 7 and Figure 5. C0, C1, C2 and Center are in different colors in Figure 5.


**Figure 5.** The location results for the GCM.

Model 2: Five-indicator clustering model (5ICM). This model contains five normalized indicators used in this paper: the longitude, the latitude, the per capita disposable income of urban residents, the permanent urban population and the population density. Euclidean distance is used to compute the distance between cities, the results are obtained by *K*-means clustering as shown in Table 8 and Figure 6. C0, C1, C2 and Center are in different colors in Figure 6.


**Table 8.** The location results for the 5ICM.

Model 3: Improved five-indicator clustering model (I5ICM). Based on the result of Model 2, this model incorporates the center-of-gravity method to obtain lower transportation costs. The results are shown in Table 9 and Figure 7. C0, C1, C2 and Center are in different colors in Figure 7.

**Figure 6.** The location results for the 5ICM.


**Table 9.** The location results for the I5ICM.

The total transportation costs of the four models are shown in Table 10.

**Table 10.** The transportation costs of the four models.


As shown in Table 10, the 2SM proposed in this research has the lowest total transportation cost. Furthermore, compared with the three traditional models, the following conclusions can be drawn:


3. If the distance function is modified reasonably and the algorithm is improved into the two-stage model (2SM) proposed in this study, the cost will be further reduced.

**Figure 7.** The location results for the I5ICM.

#### **6. Results and Conclusions**

In this research, a two-stage location selection model based on an improved clustering algorithm and the center-of-gravity method is proposed for an MFLP arising from a real-world problem. This methodology is proved to be effective and contributing to future logistics system planning strategy researches. The following conclusions are drawn as follows.

First, the more indicators introduced into LP, the more realistic the research becomes. In this paper, three socio-economic indicators, namely, economy development, traffic congestion degree and total logistics demand, are introduced in defining a distance function used in clustering, each of which could affect the decision on the site of a distribution center. Different from most of the existing researches, which only used spatial indicators, this research introduces three socio-economic indicators to evaluate a city's potential to accommodate a distribution center, described as the logistics-level score. This method provides a more comprehensive perspective for decision-makers to choose proper sites for distribution.

Second, an improved clustering algorithm is used to divide demand points into different regions. The improved algorithm redefines the traditional distance function, which could not reflect the distance result caused by the socio-economic indicators. The improved distance function takes both the positive effect of the three socio-economic indicators and the passive effect caused by the spatial indicators into consideration, and proved to be more effective than GCM, 5ICM and I5ICM in the case study.

Third, based on the methodology discussed above, this research divides 12 cities in Inner Mongolia into 3 regions and selects 1 city for each region as the regional distribution center. This is consistent with the government's current logistics center planning strategy. Such consistency indicates that this methodology could be used as an optional reference for local governments' logistics planning.

The limitation of this paper may lie in the selection of the clustering indicators. Although three socio-economic indicators are introduced in the 2SM, it is still impossible to describe the complexity of the logistics in real world. This means more socio-economic indicators or even more kinds of indicators should be introduced into LP researches. Besides, during the SFLP process, the author takes only the total transportation cost as the final objective and omits some other influential decision making factors. In subsequent research, other objectives/constraints such as carbon emission reduction and sustainable development can be added to the existing research.

**Author Contributions:** Conceptualization, J.W. and Y.B.; methodology, J.W., X.L. and Y.L.; software, X.L.; validation, J.W., X.L., Y.L. and W.Y.; formal analysis, J.W. and Y.B.; investigation, J.W. and Y.B.; resources, Y.B.; data curation, X.L. and Y.L.; writing—original draft preparation, J.W., X.L. and L.Y.; writing—review and editing, J.W., X.L. and Y.L.; visualization, X.L. and L.Y.; supervision, J.W., W.Y. and Y.B.; project administration, J.W. and Y.B.; funding acquisition, J.W. and Y.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

