Introduction

Group decision making (GDM) determines the most ideal alternative according to evaluation information expressed by a group of decision makers (DMs) [1]. Due to the advantages of overcoming the limitations of empirical knowledge and the cognitive abilities of individual DMs, it has been widely used [2]. The uncertainty and fuzziness of the decision making environment make more suitable for DMs to express evaluation information with fuzzy information than with clear numbers. Fuzzy sets (FSs) [3] have become useful tools for expressing uncertainty information in GDM problems [4,5,6], and they have been applied in semiconductor inter-bay handling systems [7], green supply chain management [8, 9], disease diagnosis [10,11,12], and many other fields [13, 14]. Thus far, different extensions of FSs have been proposed. Zhu et al. proposed dual hesitant FSs (DHFSs) in [15]. Using DHFSs, DMs can provide hesitant evaluation information from the perspective of membership degrees (MDs) and non-membership degrees (NMDs). Zhu established probabilistic hesitant FSs (PHFSs) in [16], where different MDs in PHFSs have different occurrence possibilities. However, the occurrence possibilities of membership and non-membership values are not considered in DHFSs, and NMDs and their corresponding importance are not considered in PHFSs. Therefore, DHFSs and PHFSs cannot fully express the evaluation information. Hao et al. [17] established probabilistic DHFSs (PDHFSs), which can comprehensively consider MDs, NMDs and their occurrence possibilities.

However, in GDM problems, it may be more intuitive for DMs to provide evaluation information by using preference relation (PR). PR is a judgment matrix obtained by comparing the paired alternatives. With the development of FS theory, fuzzy PR [18] has developed into other extended forms, e.g., probabilistic hesitant fuzzy PR [19], intuitionistic fuzzy PR (IFPR) [20], and hesitant fuzzy PR [21]. Based on PDHFSs, probabilistic dual hesitant fuzzy PRs (PDHFPRs) is a valuable extension of fuzzy PR. In GDM problem, some DMs think that the possible values of \(x_{i}\) prior to \(x_{j}\) are 0.4 and 0.5, while the possible values of \(x_{i}\) inferior to \(x_{j}\) are 0.2 and 0.3, and the importance of these values is different. However, due to time constraints and limited knowledge of DMs, the importance of these degrees is sometimes difficult to obtain. Therefore, there is an urgent need to study probability calculation method of incomplete PDHFPRs.

The purpose of GDM is to choose the most ideal decision from a group of alternatives. One of the most common methods is to build a model based on the consistency of PRs. Moreover, the consistency of PR is used to measure the degree of agreement among the preference values provided by the individual DMs. The lack of consistency of PR may lead to unreasonable results. To date, the research on the consistency of PR has been fruitful [21,22,23,24,25,26,27,28,29,30,31,32,33,34,35]. However, only Shao and Zhang [37] have focused on the consistency of PDHFPRs. In addition, Shao and Zhang’s method [37] has certain limitations, which is manifested in that the consistency index does not consider the risk attitude of DMs, and the weight of DMs is given subjectively. Thus, it is urgent to study the consensus of PDHFPRs and use it to solve the GDM problem.

The above analysis motivate us to further study some issues related to PDHFPRs, such as how to obtain the probability information of incomplete PDHFPRs, how to measure the consistency level of PDHFPRs, and how to use consistency to solve the GDM problem with PDHFPRs.

The rest of the sections are listed as follows: In “Related work”, the related work is briefly reviewed. In “Preliminaries”, some basic concepts are briefly reviewed. “Probability calculation method based on the multiplicative consistency of PDHFPRs” defines multiplicative consistency of PDHFPRs and proposes a probability calculation method based on the multiplicative consistency of PDHFPRs. In “A GDM approach with incomplete PDHFPRs”, a convergent algorithm to adjust the consistency level for individual PDHFPRs is built, and a GDM approach with incomplete PDHFPRs is constructed. In “Illustrative example”, the proposed approach is used to analyze the impact factors of haze weather. “Conclusion and future work” provides some conclusions.

Related work

In this section, we introduce the GDM method based on PR and review the existing research on probabilistic dual hesitant fuzzy environment that are most related to our work.

As an effective technique in decision making field, GDM based on various types of PR is widely concerned, which are summarized in Table 1. From Table 1, to date, it is rare to discuss the consistency of PDHFPRs [37].

Furthermore, probabilistic dual hesitant fuzzy environment has been widely used in GDM, which are listed in Table 2. Table 2 shows that most of the existing researches concentrated on GDM based on PDHFSs, and only Shao and Zhang [37] discussed GDM based on PDHFPRs.

Table 1 Related work on the consistency of different PRs
Table 2 GDM method in a probabilistic dual hesitant fuzzy environment

Based on the above discussion, it is meaningful to study GDM problems based on PDHFPRs. The main contributions of this paper are summarized as follows:

  1. 1.

    An optimization model based on multiplicative consistency is established to obtain the reasonable associated probabilities of PDHFEs.

  2. 2.

    For an unacceptable consistent PDHFPR, based on a linear optimization model, a convergent algorithm is investigated to improve the consistency level.

  3. 3.

    A GDM approach with incomplete PDHFPRs is constructed, and the effectiveness and applicability of the proposed decision making approach are verified by analyzing the impact factors of haze weather.

Preliminaries

This section reviews some basic concepts of DHFS, PHFS, PDHFS, and PDHFPRs:

Definition 1

[15] A DHFS F on a set of alternatives X is defined by

$$\begin{aligned} F=\{\langle x,h_{F}(x),m_{F}(x)\rangle |x\in X \} \end{aligned}$$
(1)

\(h_{F}(x)\) indicates the possible degrees to the set X of x, \(m_{F}(x)\) indicates the possible NMDs to the set X of x. In addition,

$$\begin{aligned} 0\le \alpha _{i},\beta _{j}\le 0,0\le \alpha ^{+}+\beta ^{+}\le 1, \end{aligned}$$
(2)

for \(i=1,2,\ldots , \# h_{F}(x)\), \(j=1,2,\ldots ,\# m_{F}(x)\), where \(\alpha _{i}\in h_{F}(x)\), \(h_{F}(x)=(\alpha _{1},\alpha _{2}, \ldots , \alpha _{\# h_{F}(x)})\), \(\alpha ^{+}=\bigcup _{\alpha _{i}\in h_{F}(x)}\max \{\alpha _{i}\}\), \(\beta _{j}\in m_{F}(x)\), \(m_{F}(x)=(\beta _{1},\beta _{2}\), \(\ldots , \beta _{\# m_{F}(x)})\), and \(\beta ^{+}=\bigcup _{\beta _{i}\in m_{F}(x)}\max \{\beta _{j}\}\). \(\# h_{F}(x)\), \(\# m_{F}(x)\) are the cardinal numbers of \(h_{F}(x)\) and \(m_{F}(x)\), respectively.

Definition 2

[16] A PHFS H on a set of alternatives X is defined by

$$\begin{aligned} H=\{\langle x,h_{x}(p_{x})\rangle |x\in X \} \end{aligned}$$
(3)

where \(h_{x}\in [0,1]\) indicates the possible values to the set X of x, \(p_{x}\in [0,1]\) is the associated probabilities of \(h_{x}\). In addition, \(h_{x}(p_{x})\) is called a probabilistic hesitant fuzzy element and simplified to be expressed by

$$\begin{aligned} h_{x}(p_{x})=\gamma _{i}(p_{i})|i=1,2,\ldots ,\# h_{x}, \end{aligned}$$
(4)

where \(\gamma _{i}\in h_{x}\), \(p_{i}\in p_{x}\) is the associated probability of \(\gamma _{i}\) and \(\sum _{i=1}^{\# h_x} p_i=1\). \(\#h_x\) is the cardinal number of \(h_{x}(p_{x})\).

Definition 3

[17] A PDHFS B on a set of alternatives X is defined by

$$\begin{aligned} B=\{\langle x, h(x) | p(x) , m(x) | q(x) \rangle |x\in X\} \end{aligned}$$
(5)

h(x) |p(x) and m(x) |q(x) denote some possible elements, where h(x) indicates the possible degrees to the set X of x, p(x) is the associated probabilities of h(x), m(x) indicates the possible NMDs to the set X of x, and q(x) denotes the associated probabilities of m(x). In addition,

$$\begin{aligned} 0&\le \varepsilon ,\delta \le 1, 0\le \varepsilon ^{+}+\delta ^{+}\le 1 \end{aligned}$$
(6)
$$\begin{aligned} 0&\le p_i \le 1,0\le q_i\le 1,\sum _{i=1}^{\# h}p_i=1~and~\sum _{i=1}^{\# m}q_i=1, \end{aligned}$$
(7)

where \(\varepsilon \in h(x)\), \(\varepsilon ^{+}\in h^{+}(x)=\bigcup _{\varepsilon \in h(x)}\max \{\varepsilon \}\), \(\delta \in m(x)\), \(\delta ^+\in m^{+}(x)=\bigcup _{\delta \in m(x)}\max \{\delta \}\), \(p_i\in p(x)\) and \(q_i\in q(x)\). \(\# h\) and \(\# m\) are the cardinal numbers of h(x) |p(x) and m(x) |q(x), respectively.

Remark 1

Hao et al. [17] define \(b=(h(x) | p(x),m(x) | q(x))\) as a PDHFE. For convenience, \(b=(h(x) | p(x),m(x) | q(x))\) is simplified to \(b=(h | p,m | q)\).

Remark 2

Particularly, if \(p=q\) , then a PDHFS is simplified to a DHFS. Similarly, if \(h\ne \phi \), \(m=\phi (q=\phi )\), then a PDHFS reduces to a PHFS.

Definition 4

[37] A PDHFPR D on a set of alternatives \(X=\{x_1,x_2,\ldots ,x_n\}\) is characterized by \(D=(d_{ij})_{n\times n}\), where \(d_{ij}=(h_{ij}|p_{ij},m_{ij}|q_{ij})\) is a PDHFE. \(h_{ij}\) and \(m_{ij}\) are the possible degrees to which \(x_i\) is preferred and non-preferred to \(x_j\), \(p_{ij}\) is the associated probabilities of \(h_{ij}\), and \(q_{ij}\) is the associated probabilities of \(m_{ij}\), \(d_{ij}\) should satisfy the following requirements:

  1. 1.

    For \(i=1,2,\ldots ,n\), \(d_{ij}=(\{0.5|1\},\{0.5|1\})\).

  2. 2.

    For \(i\ne j\), \(i,j=1,2,\ldots ,n\),

    1. (a)

      If \(h_{ij}\ne \phi \), \(m_{ij}\ne \phi \), then

      $$\begin{aligned} \left\{ \begin{array}{l} {h_{ij,r}=m_{ji,r},m_{ij,e}=h_{ji,e},}\\ {p_{ij,r}=q_{ji,r},q_{ij,e}=p_{ji,e},} \end{array}\right. \end{aligned}$$
      (8)

      for all \(r=1,2,\ldots ,\# h_{ij}-1\), \(e=1,2,\ldots ,\# m_{ij}-1\), where \(h_{ij,r}\in h_{ij}\) indicates the rth element in \(h_{ij}\), \(p_{ij,r}\) is the associated probability of \(h_{ij,r}\) and \(\sum _{r=1}^{\# h_{ij}}p_{ij,r}=1\), \(m_{ij,e}\in m_{ij}\) indicates the eth element in \(m_{ij}\), \(q_{ij,e}\) is the associated probability of \(m_{ij,e}\) and \(\sum _{e=1}^{\# m_{ij}}q_{ij,e}=1\), \(h_{ij,r}<h_{ij,r+1}\) and \(m_{ij,e}<m_{ij,e+1}\).

    2. (b)

      If , , then

      $$\begin{aligned}&h_{ij,r}+h_{ji,(\# h_{ij-r+1})}=1, \nonumber \\&p_{ij,r}=p_{ji,(\# h_{ij-r+1})}, \# h_{ij}=\# h_{ji}, \end{aligned}$$
      (9)

      for all \(r=1,2,\ldots ,\# h_{ij}-1,\), where \(h_{ij,r}\in h_{ij}\) indicates the rth element in \(h_{ij}\), \(p_{ij,r}\) is the associated probability of \(h_{ij,r}\) and \(\sum _{r=1}^{\# h_{ij}}p_{ij,r}=1\), and \(h_{ij,r}<h_{ij,r+1}\).

    3. (c)

      If , , then

      $$\begin{aligned}&m_{ij,e}+m_{ji,(\# g_{ij-e+1})}\nonumber \\&\quad =1, q_{ij,e}=q_{ji,(\# g_{ij-e+1})}, \# m_{ij}=\# m_{ji}, \end{aligned}$$
      (10)

      for all \(e=1,2,\ldots ,\# m_{ij}-1\), where \(m_{ij,e}\in m_{ij}\) indicates the eth element in \(m_{ij}\), \(q_{ij,e}\) is the associated probability of \(m_{ij,e}\) and \(\sum _{e=1}^{\# m_{ij}}q_{ij,e}=1\), and \(m_{ij,e}<m_{ij,e+1}\).

Example 1

The decision information of alternatives \(\{x_1,x_2,x_3\}\) is provided by the DM as follows:

$$\begin{aligned} D=\left( \begin{array}{ccc} (\{0.5|1\},\{0.5|1\}) &{} \left. \begin{array}{c} (\{0.5|p_{12,1},0.6|p_{12,2}\}, \\ \{0.2|q_{12,1},0.4|q_{12,2}\}) \\ \end{array}\right. &{} \left. \begin{array}{c} (\{0.4|p_{13,1},0.6|p_{13,2}\}, \\ \{0.2|q_{13,1},0.3|q_{13,2}\}) \\ \end{array}\right. \\ \left. \begin{array}{c} (\{0.2|q_{12,1},0.4|q_{12,2}\}, \\ \{0.5|p_{12,1},0.6|p_{12,2}\}) \\ \end{array}\right. &{} (\{0.5|1\},\{0.5|1\}) &{} \left. \begin{array}{c} (\{0.5|p_{23,1},0.7|p_{23,2}\}, \\ \{0.1|q_{23,1},0.3|q_{23,2}\}) \\ \end{array}\right. \\ \left. \begin{array}{c} (\{0.2|q_{13,1},0.3|q_{13,2}\}, \\ \{0.4|p_{13,1},0.6|p_{13,2}\}) \\ \end{array}\right. &{} \left. \begin{array}{c} (\{0.1|q_{23,1},0.3|q_{23,2}\}, \\ \{0.5|p_{23,1},0.7|p_{23,2}\}) \\ \end{array}\right. &{} (\{0.5|1\},\{0.5|1\}) \\ \end{array}\right) \end{aligned}$$

For \(i\ne j\) and \(i,j=1,2,3\), there are \(0\le p_{ij,r},q_{ij,r}\le 1\), \(\sum _{r=1}^2p_{ij,r}=1\) and \(\sum _{r=1}^2q_{ij,r}=1\).

From Definition 4, D is a PDHFPR.

Probability calculation method based on the multiplicative consistency of PDHFPRs

The limited knowledge of DMs makes it hard for DMs to provide accurate associated probabilities of PDHFEs through subjective evaluation. This section proposes probability calculation method based on the multiplicative consistency of PDHFPRs.

Multiplicative consistency of PDHFPRs

Definition 5

[20] An IFPR on the set \(X=\{x_1,x_2,\ldots ,x_n\}\) is defined by \(A=(a_{ij})_{n\times n}\), where \(a_{ij}=\langle u_{ij},v_{ij}\rangle \) is an intuitionistic fuzzy value (IFV), \(u_{ij}\) indicates the degree that \(x_i\) is prior to \(x_j\), \(v_{ij}\) indicates the degree that \(x_i\) is not prior to \(x_{j}\) and for \(i,j=1,2,\ldots ,n\),

$$\begin{aligned} 0\le u_{ij}+v_{ij}\le 1, u_{ij}=v_{ji}, v_{ij}=u_{ji}, u_{ii}=v_{ii}=0.5. \end{aligned}$$
(11)

Definition 6

[20] Let \(a_{ij}=\langle u_{ij}, v_{ij}\rangle \) be an IFV. The score function of the IFV is defined by

$$\begin{aligned} \triangle (a_{ij})=u_{ij}-v_{ij}. \end{aligned}$$
(12)

Definition 7

[27] An IFPR on the set \(X=\{x_1,x_2,\ldots ,x_n\}\) is defined by \(A=(a_{ij})_{n\times n}\), if there is a normalized intuitionistic fuzzy weight vector \(\phi =(\phi _1,\phi _2,\ldots ,\phi _n)^T\) such that

$$\begin{aligned} a_{ij}=\langle u_{ij},v_{ij}\rangle =\left\{ \begin{array}{ll} \langle 0.5,0.5\rangle , &{} i=j \\ \langle \sqrt{\phi _{i\mu }\cdot \phi _{j\nu }}, \sqrt{\phi _{i\nu }\cdot \phi _{j\mu }}\rangle , &{} i\ne j \end{array}\right. \end{aligned}$$
(13)

then, A is a multiplicatively consistent IFPR, where \(\phi _{i}=\langle \phi _{i\mu }, \phi _{i\nu }\rangle \), \(\phi _{i\mu },\phi _{i\nu }\in (0,1]\), \(\phi _{i\mu }+\phi _{i\nu }\le 1\), \(\sum _{j\ne i}\phi _{j\mu }\le \phi _{i\nu }\), and \(\sum _{j\ne i}\phi _{j\nu }\le \phi _{i\mu }+n-2\), \(i,j=1,2,\ldots ,n\).

The following introduces the order consistency and multiplicative consistency of PDHFPRs. For complete PDHFPRs, order consistency provides a quick and effective way to obtain the order of alternatives.

Definition 8

A PDHFPR \(D=(d_{ij})_{n\times n}\) on the set \(X=\{x_1,x_2,\ldots ,x_n\}\) is an order consistent PDHFPR, if \(d_{sj}\ge d_{tj}\) for all \(j=1,2,\ldots ,n\), where \(s,t=1,2,\ldots ,n\).

Definition 9

A PDHFPR on the set \(X=\{x_1,x_2,\ldots ,x_n\}\) is represented by \(D=(d_{ij})_{n\times n},\) where \(d_{ij}=(h_{ij}|p_{ij},m_{ij}|q_{ij})\). \(U_{d_{ij}}\) and \(V_{d_{ij}}\) are denoted by

$$\begin{aligned} U_{d_{ij}}=\sum _{r=1}^{\# h_{ij}}h_{ij,r}\cdot p_{ij,r},V_{d_{ij}} =\sum _{e=1}^{\# m_{ij}}m_{ij,e}\cdot q_{ij,e}, \end{aligned}$$
(14)

where \(h_{ij,r}\in h_{ij}\) indicates the rth element in \(h_{ij}\), \(p_{ij,r}\) is the associated probability of \(h_{ij,r}\) and \(\sum _{r=1}^{\# h_{ij}}p_{ij,r}=1\), \(m_{ij,e}\in m_{ij}\) indicates the eth element in \(m_{ij}\), \(q_{ij,e}\) is the associated probability of \(m_{ij,e}\) and \(\sum _{e=1}^{\# m_{ij}}q_{ij,e}=1\).

Based on Definition 4, for \(\forall i,j=1,2,\ldots ,n\), then \(\# h_{ij}=\# m_{ji}\), \(\# m_{ij}=\# h_{ji}\) and \(U_{d_{ii}}=V_{d_{ii}}=0.5\). Furthermore,

$$\begin{aligned} U_{d_{ij}}= & {} \sum _{r=1}^{\# h_{ij}}h_{ij,r}\cdot p_{ij,r} =\sum _{e=1}^{\# m_{ji}}m_{ji,r}\cdot q_{ji,r}=V_{d_{ji}},\\ V_{d_{ji}}= & {} \sum _{e=1}^{\# m_{ij}}m_{ij,e}\cdot q_{ij,e} =\sum _{e=1}^{\# h_{ji}}h_{ji,e}\cdot p_{ji,e}=U_{d_{ji}} \end{aligned}$$

Since \(h_{ij,r},m_{ij,e},p_{ij,e},q_{ij,e}\in [0,1]\), then

$$\begin{aligned} 0\le U_{d_{ij}}= & {} \sum _{r=1}^{\# h_{ij}}h_{ij,r}\cdot p_{ij,r} \le 1,0\le V_{d_{ij}}\\= & {} \sum _{e=1}^{\# m_{ij}}h_{ij,e} \cdot p_{ij,e}\le 1.\\ U_{d_{ij}}+V_{d_{ij}}= & {} \sum _{r=1}^{\# h_{ij}}h_{ij,r}\cdot p_{ij,r} +\sum _{e=1}^{\# m_{ij}}h_{ij,e}\cdot p_{ij,e}\\\le & {} \max _{h_{ij,r}\in h_{ij}}\{h_{ij,r}\}\cdot \sum _{r=1}^{\# h_{ij}}p_{ij,r}\\&+\max _{m_{ij,e}\in m_{ij}}\{m_{ij,e}\}\cdot \sum _{e=1}^{\# m_{ij}}q_{ij,e}\le 1. \end{aligned}$$

Thus, \(E_D=\langle U_{d_{ij}},V_{d_{ij}}\rangle \) is an IFPR.

Definition 10

A PDHFPR on the set \(X=\{x_1,x_2,\ldots ,x_n\}\) is denoted as \(D=(d_{ij})n\times n,\) where \(d_{ij}=(h_{ij}|p_{ij},m_{ij}|q_{ij})\). If a normalized weight vector \(W=(w_1,w_2,\ldots ,w_n)^T\) exists such that

$$\begin{aligned} \langle U_{d_{ij}}V_{d_{ij}}\rangle&=\left\langle \sum _{r=1}^{\sharp h_{ij}} h_{ij,r}\cdot p_{ij,r}\sum _{e=1}^{\sharp m_{ij}}h_{ij,e}\cdot p_{ij,e}\right\rangle \nonumber \\&=\left\{ \begin{array}{ll} \langle 0.5,0.5\rangle , &{} i=j, \\ \langle \sqrt{w_{i\mu }\cdot w_{j\nu }}, \sqrt{w_{i\nu }\cdot w_{j\mu }}\rangle , &{} i\ne j, \end{array}\right. \end{aligned}$$
(15)

then PDHFPR D is multiplicative consistent, where \(w_i=\langle w_{i\mu },w_{i\nu }\rangle \) satisfies \(w_{i\mu },w_{i\nu }\in (0,1]\), \(w_{i\mu }+w_{i\nu }\le 1\), \(\sum _{j\ne i}w_{j\mu }\le w_{i\nu }\), and \(\sum _{j\ne i}w_{j\nu }\le w_{i\mu }+n-2\), \(i,j=1,2,\ldots ,n\).

Probability calculation method

Due to the complexity and uncertainty of the decision problem, the DMs cannot provide the accurate associated probabilities of PDHFEs, and it is challenging to give PDHFPRs with multiplicative consistency. At this time, the decision making result is not rational. When a PDHFPR does not satisfy multiplicative consistency, Eq. (15) does not hold, which implies that

$$\begin{aligned}&\sum _{r=1}^{\# h_{ij}}h_{ij,r}\cdot p_{ij,r}\ne \sqrt{w_{i\mu } \cdot w_{j\nu }}~~or\nonumber \\&\sum _{e=1}^{\# m_{ij}}m_{ij,r}\cdot q_{ij,e} \ne \sqrt{w_{i\nu }\cdot w_{j\mu }} \end{aligned}$$
(16)

The non-negative deviation variables \(\tau _{ij}\) and \(\rho _{ij}\) are defined as follows:

$$\begin{aligned} \tau _{ij}= & {} \left| \sum _{r=1}^{\# h_{ij}}h_{ij,r}\cdot p_{ij,r}-\sqrt{w_{i\mu }\cdot w_{j\nu }}\right| \\ \rho _{ij}= & {} \left| \sum _{e=1}^{\# m_{ij}}m_{ij,e}\cdot q_{ij,e}-\sqrt{w_{i\nu }\cdot w_{j\mu }}\right| \end{aligned}$$

The smaller the value of the non-negative deviation variables \(\tau _{ij}\) and \(\rho _{ij}\) , the higher the level of multiplicative consistency of PDHFPR D.

From \(\sum _{r=1}^{\# h_{ij}}h_{ij,r}\cdot p_{ij,r}=\sum _{r=1}^{\# m_{ji}}m_{ji,r}\cdot q_{ji,r}\) and \(\sum _{e=1}^{\# m_{ij}}m_{ij,e}\cdot q_{ij,r}=\sum _{e=1}^{\# h_{ji}}h_{ji,e}\cdot p_{ji,e}\), we can obtain \(\tau _{ij}=|\sum _{r=1}^{\# h_{ij}}h_{ij,r}\cdot p_{ij,r}-\sqrt{w_{i\mu }\cdot w_{j\nu }}|=|\sum _{e=1}^{\# m_{ij}}m_{ji,e}\cdot q_{ji,e}-\sqrt{w_{i\nu }\cdot w_{j\mu }}|=\rho _{ji}\) and \(\tau _{ij}=\rho _{ji}\) Thus, we only need to consider the upper triangular elements of PDHFPR and can simplify the optimization model as follows:

(M-1)

$$\begin{aligned} \min J&= \sum _{i=1}^{n-1}\sum _{j=i+1}^n(\tau _{ij}+\rho _{ij})\nonumber \\&=\sum _{i=1}^{n-1}\sum _{j=i+1}^n\left( \left| \sum _{r=1}^{\# h_{ij}}h_{ij,r} \cdot p_{ij,r}-\sqrt{w_{i\mu }\cdot w_{j\nu }}\right| \right. \nonumber \\&\quad \left. +\left| \sum _{e=1}^{\# m_{ij}}m_{ij,e}\cdot q_{ij,e} -\sqrt{w_{i\nu }\cdot w_{j\mu }}\right| \right) \nonumber \\&\quad \mathrm{{s.t.}}\left\{ \begin{array}{l} p_{ij,r},q_{ij,e}\in [0,1], \displaystyle \sum _{r=1}^{\# h_{ij}}p_{ij,r}=1, \displaystyle \sum _{m=1}^{\# m_{ij}}q_{ij,e}=1 \\ w_{i\mu },w_{j\nu }\in (0,1], \displaystyle \sum _{j=1,j\ne i}^nw_{j\mu }\le w_{i\nu },\\ w_{i\mu }+w_{j\nu }\le 1, \displaystyle \sum _{j=1,j\ne i}^nw_{j\nu }\le w_{i\mu }+n-2\\ r=1,2,\ldots ,\# h_{ij}, e=1,2,\ldots ,\# m_{ij} \end{array}\right. \end{aligned}$$
(17)

Then, there exist four non-negative deviation variables \(o_{ij}^-,o_{ij}^+,\varsigma _{ij}^-,\varsigma _{ij}^+\) satisfying the following conditions \(o_{ij}^-\cdot o_{ij}^+=0\) and \(\varsigma _{ij}^-\cdot \varsigma _{ij}^+=0\) such that model (M-1) can be expressed as follows:

(M-2)

$$\begin{aligned}&\min J =\sum _{i=1}^{n-1}\sum _{j=i+1}^n(g_{ij}^-\cdot o_{ij}^-+g_{ij}^+ \cdot o_{ij}^++n_{ij}^+\cdot \varsigma _{ij}^++n_{ij}^-\cdot \varsigma _{ij}^-)\nonumber \\&\mathrm{{s.t.}}\left\{ \begin{array}{l} \displaystyle \sum _{r=1}^{\# h_{ij}}h_{ij,r}\cdot p_{ij,r}+g_{ij}^-\cdot o_{ij}^{-}-g_{ij}^+ \cdot o_{ij}^{+}=\sqrt{w_{i\mu }\cdot w_{j\nu }}\\ \displaystyle \sum _{e=1}^{\# m_{ij}}m_{ij,e}\cdot q_{ij,e}+n_{ij}^-\cdot \varsigma _{ij}^{-} -n_{ij}^+\cdot \varsigma _{ij}^{+}=\sqrt{w_{i\nu }\cdot w_{j\mu }} \\ o_{ij}^-,o_{ij}^+,\varsigma _{ij}^-,\varsigma _{ij}^+\ge 0,o_{ij}^-\cdot o_{ij}^+=0, \varsigma _{ij}^-\cdot \varsigma _{ij}^+=0 \\ p_{ij,r},q_{ij,e}\in [0,1], \displaystyle \sum _{r=1}^{\# h_{ij}}p_{ij,r}=1, \displaystyle \sum _{e=1}^{\# m_{ij}}q_{ij,e}=1 \\ w_{i\mu },w_{j\nu }\in (0,1], \displaystyle \sum _{j=1,j\ne i}^nw_{j\mu }\le w_{i\nu },\\ w_{i\mu }+w_{j\nu }\le 1, \displaystyle \sum _{j=1,j\ne i}^nw_{j\nu }\le w_{i\mu }+n-2 \\ i<j, i,j=1,2,\ldots ,n, r=1,2,\ldots ,\\ \# h_{ij}, j=1,2,\ldots ,\# m_{ij} \end{array}\right. \end{aligned}$$
(18)

Here, \(g_{ij}^-\), \(g_{ij}^+\) , \(n_{ij}^-\) and \(n_{ij}^+\) are the weights of non-negative deviation variables. In this paper, the weights of these deviation variables were chosen to be equal, namely, \(g_{ij}^-=g_{ij}^+=n_{ij}^-=n_{ij}^+=1\). Furthermore, a simplified model is as follows:

(M-3)

$$\begin{aligned}&\min J=\sum _{i=1}^{n-1}\sum _{j=i+1}^n(o_{ij}^- +o_{ij}^++\varsigma _{ij}^++\varsigma _{ij}^-)\nonumber \\&\mathrm{{s.t.}}\left\{ \begin{array}{l} \displaystyle \sum _{r=1}^{\# h_{ij}}h_{ij,r}\cdot p_{ij,r}-o_{ij}^{+}+o_{ij}^{-}=\sqrt{w_{i\mu }\cdot w_{j\nu }},\\ \displaystyle \sum _{e=1}^{\#m_{ij}}m_{ij,e}\cdot q_{ij,e}-\varsigma _{ij}^{+}+\varsigma _{ij}^{-}=\sqrt{w_{i\nu }\cdot w_{j\mu }}\\ o_{ij}^-,o_{ij}^+,\varsigma _{ij}^-,\varsigma _{ij}^+\ge 0,o_{ij}^-\cdot o_{ij}^+=0, \varsigma _{ij}^-\cdot \varsigma _{ij}^+=0 \\ p_{ij,r},q_{ij,e}\in [0,1], \displaystyle \sum _{r=1}^{\# h_{ij}}p_{ij,r}=1, \displaystyle \sum _{m=1}^{\# m_{ij}}q_{ij,e}=1 \\ w_{i\mu },w_{j\nu }\in (0,1], \displaystyle \sum _{j=1,j\ne i}^nw_{j\mu }\le w_{i\nu },\\ w_{i\mu }+w_{j\nu }\le 1, \displaystyle \sum _{j=1,j\ne i}^nw_{j\nu }\le w_{i\mu }+n-2 \\ i<j, i,j=1,2,\ldots ,n, r=1,2,\ldots ,\# h_{ij}, e=1,2,\ldots ,\# m_{ij} \end{array}\right. \end{aligned}$$
(19)

From the above discussion, solving model (M-3) can obtain the accurate associated probabilities \(p_{ij,r}\) and \(q_{ij,e} (i,j=1,2,\ldots ,n, r=1,2,\ldots ,\#h_{ij}, e=1,2,\ldots ,\# m_{ij})\) of the PDHFE; the optimal variables \(\widetilde{o_{ij}^+}\), \(\widetilde{o_{ij}^-}\), \(\widetilde{\varsigma _{ij}^+}\), \(\widetilde{\varsigma _{ij}^-}\) and the optimal normalized weight vector \(w^{*}\).

If \(\widetilde{J}=0\), since \(\widetilde{o_{ij}^+}, \widetilde{o_{ij}^-}, \widetilde{\varsigma _{ij}^+}, \widetilde{\varsigma _{ij}^-}\ge 0\) and using the constraints of model (M-3), we can get

$$\begin{aligned}&\sum _{r=1}^{\# h_{ij}}h_{ij,r}\cdot p_{ij,r}=\sqrt{\widetilde{w_{i\mu }} \cdot \widetilde{w_{j\nu }}}, \nonumber \\&\sum _{e=1}^{\# m_{ij}}m_{ij,e}\cdot q_{ij,e} =\sqrt{\widetilde{w_{i\nu }}\cdot \widetilde{w_{j\mu }}} \end{aligned}$$
(20)

Thus, \(D=(d_{ij})_{n\times n}\) is a multiplicative consistent PDHFPR.

The rationality of the proposed method is proven by Example 1. First, based on model (M-3), the corresponding model can be obtained. Using the MATLAB Optimization Toolbox to solve the corresponding model, the following expression is obtained:

$$\begin{aligned}&p_{12,1}=0.4778,p_{12,2}=0.5222,p_{13,1}=0,\\&p_{13,2}=1,p_{23,1}=0.5846,\\&p_{23,2}=0.4254, q_{12,1}=0.0909,\\&q_{12,2}=0.9091,q_{13,1}=0.8322,\\&q_{13,2}=0.1678,q_{23,1}=0.1708,\\&q_{23,2}=0.8292,o_{12}^-=o_{12}^+=\varsigma _{12}^-=\varsigma _{12}^+=\\&o_{13}^-=o_{13}^+=\varsigma _{13}^-=\varsigma _{13}^+=o_{23}^+=\varsigma _{23}^-=0,\\&o_{23}^-=0.0714, \varsigma _{23}^+=0.0015,\\&\widetilde{w_1}=\langle 0.4565,0.4391\rangle ,\widetilde{w_2} =\langle 0.3320,0.6680\rangle , \\&\widetilde{w_3}=\langle 0.1070,0.7886\rangle . \end{aligned}$$

Therefore, the complete PDHFPR D in Example 1 is as follows:

$$\begin{aligned} D=\left( \begin{array}{ccc} (\{0.5|1\},\{0.5|1\}) &{} \left. \begin{array}{c} (\{0.5|0.4778,0.6|0.5222\}, \\ \{0.2|0.0909,0.4|0.9091\}) \\ \end{array}\right. &{} \left. \begin{array}{c} (\{0.6|1\}, \\ \{0.2|0.8322,0.3|0.1678\}) \\ \end{array}\right. \\ \left. \begin{array}{c} (\{0.2|0.0909,0.4|0.9091\}, \\ \{0.5|0.4778,0.6|0.5222\}) \\ \end{array}\right. &{} (\{0.5|1\},\{0.5|1\}) &{} \left. \begin{array}{c} (\{0.5|0.5846,0.7|0.4154\}, \\ \{0.1|0.1708,0.3|0.8292\}) \\ \end{array}\right. \\ \left. \begin{array}{c} (\{0.2|8322,0.3|0.1678\}\{0.6|1\}, \\ \{0.6|1\}) \\ \end{array}\right. &{} \left. \begin{array}{c} (\{0.1|0.1708,0.3|0.8292\}, \\ \{0.5|0.5846,0.7|0.4154\}) \\ \end{array}\right. &{} (\{0.5|1\},\{0.5|1\}) \\ \end{array}\right) \end{aligned}$$

In addition, \(\triangle (\widetilde{w_{1}})=0.0174\), \(\triangle (\widetilde{w_{2}})=-0.336\), \(\triangle (\widetilde{w_{3}})=-0.6816\) are obtained, and the ranking of the alternatives is \(x_1\succ x_2\succ x_3\).

According to Definition 8, the complete PDHFPR D satisfies \(d_1\succ d_2\succ d_3\), for all \(j=1,2,3\). Then, the ranking is \(x_1\succ x_2\succ x_3\). Therefore, the feasibility of the probability calculation approach based on the consistency of PDHFPRs is verified.

A GDM approach with incomplete PDHFPRs

A convergent method to improve the consistency of individual PDHFPRs

In the decision making process, PDHFPRs are difficult to achieve full consistency. Considering the risk attitude of DMs, the consistency level for PDHFPRs is tested by a weighted consistency index.

Definition 11

If \(D=(d_{ij})_{n\times n}\) is a PDHFPR, \(\gamma _0\ge 0\) is the given consistency index threshold. A weighted consistency index can be denoted by

$$\begin{aligned} \mathrm{{WCI}}(D)=\lambda \cdot \mathrm{{OCI}}(D)+(1-\lambda )\cdot \mathrm{{PCI}}(D). \end{aligned}$$
(21)

\(\mathrm{{OCI}}(D)=\frac{2}{n\cdot (n-1)} \sum _{i=1}^{n-1} \sum _{j=i+1}^n(o_{ij}^-+o_{ij}^++\varsigma _{ij}^++\varsigma _{ij}^-)\) is the DM’s optimistic consistency index, denoting that attitude of the DM is optimistic; \(\mathrm{{PCI}}(D)=\max \{o_{ij}^-+o_{ij}^+ +\varsigma _{ij}^+ +\varsigma _{ij}^-\}\) is the DM’s pessimistic consistency index, denoting that the DM’s attitude is pessimistic; and \(0\le \lambda \le 1\) s the weight coefficient of \(\mathrm{{OCI}}(D)\), which represents the risk attitude of the DM.

If \(\mathrm{{WCD}}(D)\le \gamma _0\), the consistency of PDHFPR \(D=(d_{ij})_{n\cdot n}\) is acceptable. For an unacceptable consistent PDHFPR, the decision result obtained from it is usually unreliable. Next, a convergent automatic iterative algorithm is used for adjusting the consistency level, and the corresponding ranking weight vector is calculated to make a decision.

Algorithm 1

Input: The initial PDHFPR \(D=(d_{ij})_{n\times n}\), the given consistency index threshold \(\gamma _0\), the weight coefficient \(\lambda \), the adjusted coefficient \(\theta \), and the maximum number of iterations \(t_{max}\). Note that \(d_{ij}=(h_{ij}|p_{ij},m_{ij}|q_{ij})\) is a PDHFE, \(h_{ij,r}\in h_{ij}\) indicates the rth element in \(h_{ij}\), \(p_{ij,r}\) is the associated probability of \(h_{ij,r}\) satisfying \(\sum \limits _{r=1}^{\# h_{ij}}p_{ij,r}=1\), \(m_{ij,e}\in m_{ij}\) indicates the eth element in \(m_{ij}\), and \(q_{ij,e}\) is the associated probability of \(m_{ij,e}\) satisfying \(\sum \limits _{e=1}^{\# m_{ij}}q_{ij,e}=1\). Moreover, the associated probabilities \(p_{ij,r}\) and \(q_{ij,e}\) of PDHFEs are unknown.

Output: The acceptable consistent PDHFPR \(D^{*}=(d_{ij}^{*})_{n\times n}\), the optimal normalized weight vector \(w^{*}\) , and the weighted consistency index \(\mathrm{{WCI}}(D^{*})\).

Step 1: Let \(D^{(t)}=(d_{ij}^{(t)})_{n\times n}\), where \(d_{ij}^{(t)}=(h_{ij}^{(t)}|p_{ij},m_{ij}^{(t)}|q_{ij})\) is a PDHFE, and when \(t=0\), \(D^{(0)}=D\).

Step 2: Model (M-3) is used to obtain the accurate associated probabilities \(p_{ij,r}\) and \(q_{ij,e}\), \((i,j=1,2,\ldots ,n, r=1,2,\ldots ,\# h_{ij}, e=1,2,\ldots ,\# m_{ij})\) of PDHFEs, the optimal variables \(\widetilde{o_{ij}^{(t)+}}, \widetilde{o_{ij}^{(t)-}},\widetilde{\varsigma _{ij}^{(t)+}}, \widetilde{\varsigma _{ij}^{(t)-}}\) and the optimal normalized weight vector \(\widetilde{w^{(t)}}=(\widetilde{w_1^{(t)}},\widetilde{w_2^{(t)}}, \ldots ,\widetilde{w_n^{(t)}})\).

Step 3: The weighted consistency index \(\mathrm{{WCI}}(D^{(t)})\) is calculated according to the risk attitude of DMs:

$$\begin{aligned} \mathrm{{WCI}}(D^{(t)})=\lambda \cdot \mathrm{{OCI}}(D^{(t)})+(1-\lambda )\cdot \mathrm{{PCI}}(D^{(t)}). \end{aligned}$$
(22)

Step 4: If \(\mathrm{{WCI}}(D^{(t)})\le \gamma _0\) or \(t\ge t_{max}\) , turn to Step 6; or else, proceed to the next step.

Step 5: Use Eq. (23) to determine where DMs need to modify preference values.

$$\begin{aligned} \mathrm{{LOC}}&= \bigg \{(z,c)|\widetilde{o_{zc}^{(t)+}}+\widetilde{o_{zc}^{(t)-}} +\widetilde{\varsigma _{zc}^{(t)+}}+\widetilde{\varsigma _{zc}^{(t)-}}\nonumber \\&= \max _{(i,j)}\{\widetilde{o_{ij}^{(t)+}}+\widetilde{o_{ij}^{(t)-}} +\widetilde{\varsigma _{ij}^{(t)+}}+\widetilde{\varsigma _{ij}^{(t)-}}|i,j\in N,i<j\}\bigg \} \end{aligned}$$
(23)
$$\begin{aligned}&\mathrm{{Let}}\left\{ \begin{array}{l} h_{ij,r}^{(t+1)}=\left\{ \begin{array}{l} \min (\theta \cdot \widetilde{h_{zc,r}^{(t)}}+(1-\theta )\cdot h_{zc,r}^{(t)},\\ 1-\max _{g_{zc,r}^{(t+1)}\in g_{zc}^{(t+1)}}\{m_{zc,e}^{(t+1)}\}),i=z,j=c, \\ h_{ij,r}^{(t)},i,j\in N,i<j,i\ne z,j\ne c, \end{array}\right. \\ m_{ij,e}^{(t+1)}=\left\{ \begin{array}{l} \min (\theta \cdot \widetilde{m_{zc,e}^{(t)}}+(1-\theta )\cdot m_{zc,e}^{(t)},\\ 1-\max _{h_{zc,r}^{(t+1)}\in h_{zc}^{(t+1)}}\{h_{zc,r}^{(t+1)}\}),i=z,j=c, \\ m_{ij,e}^{(t)},i,j\in N,i<j,i\ne z,j\ne c, \end{array}\right. \\ \end{array}\right. \end{aligned}$$
(24)

where \(\widetilde{h_{zc,r}^{(t)}}=h_{zc,r}^{(t)} +\widetilde{o_{zc}^{(t)-}}-\widetilde{o_{zc}^{(t)+}}\) and \(\widetilde{m_{zc,e}^{(t)}}=m_{zc,e}^{(t)} +\widetilde{\varsigma _{zc}^{(t)-}}-\widetilde{\varsigma _{zc}^{(t)+}}\).

Furthermore, by Definition 4, the lower triangular elements of PDHFPR \(D^{(t+1)}\) can be obtained, and then the improved PDHFPR \(D^{(t+1)}=(d_{ij}^{(t+1)})_{n\times n}\) can be obtained. Let \(t=t+1\) , and return to Step 2.

Step 6: Let \(D^{*}=D^{(t)}\) and \(w^{*}=\widetilde{w^{(t)}}\), and get the ranking of the alternatives according to \(w^{*}\).

Next, Theorem 1 is used to prove that Algorithm 1 is convergent.

Theorem 1

If \(D=(d_{ij})_{n\times n}\) is a PDHFPR, where \(d_{ij}=(h_{ij}|p_{ij},m_{ij}|q_{ij})\), \(\{D^{(t)}\}\) are PDHFPR sequences in Algorithm 1, \( \mathrm{{WCI}}(D^{(t)})\) is the weighted index of \(D^{(t)}\). Then, for each iteration t,

$$\begin{aligned} \mathrm{{WCI}}(D^{(t+1)})\le \mathrm{{WCI}}(D^{(t)}). \end{aligned}$$

Proof

According to the Eqs. (23) and (24), we have \(i\ne z,j\ne c, \widetilde{o_{ij}^{(t+1)+}}=\widetilde{o_{ij}^{(t)+}}, \widetilde{o_{ij}^{(t+1)-}}=\widetilde{o_{ij}^{(t)-}}, \widetilde{\varsigma _{ij}^{(t+1)+}}=\widetilde{\varsigma _{ij}^{(t)+}}\), and \(\widetilde{\varsigma _{ij}^{(t+1)-}}=\widetilde{\varsigma _{ij}^{(t)-}}\).

If \(\widetilde{o_{ij}^{(t)-}}>0, \widetilde{\varsigma _{ij}^{(t)-}} >0\), then \(\widetilde{o_{ij}^{(t)+}}=0\) and \(\widetilde{\varsigma _{ij}^{(t)+}}=0\). There are two possible conditions as follows:

  1. (a)

    If \(h_{zc,r}^{(t+1)}=\theta \cdot \widetilde{h_{zc,r}^{(t)}} +(1-\theta )\cdot h_{zc,r}^{(t)}\) and \(m_{zc,e}^{(t+1)}=\theta \cdot \widetilde{m_{zc,e}^{(t)}}+(1-\theta )\cdot m_{zc,e}^{(t)}\), then \(h_{zc,r}^{(t+1)}>h_{zc,r}^{(t)}\) and \(m_{zc,e}^{(t+1)}>m_{zc,e}^{(t)}\). Thus, \((\widetilde{o_{zc}^{(t+1)+}} +\widetilde{o_{zc}^{(t+1)-}}+\widetilde{\varsigma _{zc}^{(t+1)+}} +\widetilde{\varsigma _{zc}^{(t+1)-}})< (\widetilde{o_{zc}^{(t)+}} +\widetilde{o_{zc}^{(t)-}}+\widetilde{\varsigma _{zc}^{(t)+}} +\widetilde{\varsigma _{zc}^{(t)-}})\).

  2. (b)

    If \(h_{zc,r}^{(t+1)}=1-\max _{m_{zc,e}^{(t+1)} \in m_{zc}^{(t+1)}}\{m_{zc,e}^{(t+1)}\}\) or

    $$\begin{aligned} m_{zc,e}^{(t+1)} =1-\max _{h_{zc,r}^{(t+1)}\in h_{zc}^{(t+1)}}\{h_{zc,r}^{(t+1)}\}, \end{aligned}$$

    then \(h_{zc,r}^{(t+1)}\ge h_{zc,r}^{(t)}\) and \(m_{zc,e}^{(t+1)} \ge m_{zc,e}^{(t)}\) and

    $$\begin{aligned}&(\widetilde{o_{zc}^{(t+1)+}}+\widetilde{o_{zc}^{(t+1)-}} +\widetilde{\varsigma _{zc}^{(t+1)+}}+\widetilde{\varsigma _{zc}^{(t+1)-}})\\&\quad \le \widetilde{(}{o_{zc}^{(t)+}}+\widetilde{o_{zc}^{(t)-}} +\widetilde{\varsigma _{zc}^{(t)+}}+\widetilde{\varsigma _{zc}^{(t)-}}). \end{aligned}$$

In other cases, the same can be proved. Therefore,

$$\begin{aligned}&(\widetilde{o_{zc}^{(t+1)+}}+\widetilde{o_{zc}^{(t+1)-}} +\widetilde{\varsigma _{zc}^{(t+1)+}}+\widetilde{\varsigma _{zc}^{(t+1)-}})\\&\quad \le (\widetilde{o_{zc}^{(t)+}}+\widetilde{o_{zc}^{(t)-}} +\widetilde{\varsigma _{zc}^{(t)+}}+\widetilde{\varsigma _{zc}^{(t)-}}). \end{aligned}$$

We have \(\mathrm{{OCI}}(D^{(t+1)})\le \mathrm{{OCI}}(D^{(t)})\), \(\mathrm{{PCI}}(D^{(t+1)})\le \mathrm{{PCI}}(D^{(t)})\) and

$$\begin{aligned} \mathrm{{WCI}}(D^{(t+1)})\le \mathrm{{WCI}}(D^{(t)}). \end{aligned}$$

\(\square \)

A GDM approach with incomplete PDHFPRs

Faced with complicated decision making problems, the decision results obtained based on individual PDHFPRs may not be accurate due to the limitations of the expertise and experience of a DM, and a group of DMs need to make decisions at the same time. Because different DMs have various professional knowledge and experience, it is more reasonable for different DMs to give different weights. In other words, the weights of DMs should be proportional to the consistency level for the corresponding PDHFPR.

Suppose that \(X=\{x_1,x_2,\ldots ,x_n\}\) is an alternative set, \(P=\{p_1,p_2,\ldots ,p_m\}\) represents the set of DMs, and \(\psi =\{\psi _1,\psi _2,\ldots ,\psi _m\}\) is the weight vector of DMs. The detailed steps of the GDM approach based on the consistency of PDHFPRs are as follows:

Algorithm 2

Input: A list of PDHFPRs \(D_1,D_2,\ldots ,D_m\), where the associated probabilities of the corresponding PDHFEs are unknown, the given consistency index threshold \(\gamma _0\), the weight coefficient \(\lambda \), the adjusted coefficient \(\theta \) and the maximum number of iterations \(t_\mathrm{{max}}\).

Output: Improved complete PDHFPRs \(D_l^{*}(l=1,2,\ldots ,m)\), the ranking of \(x_i(i=1,2,\ldots ,n)\) and the overall priority weights of alternatives \(\varpi _i\).

Step 1: Construct PDHFPR series \(D_{l}=(d_{ij}^{l})_{n\times n}(1\le l\le m)\) on the condition of the decision information provided by the DMs, where \(d_{ij}^{l}=(h_{ij}^{l}|p_{ij}^l,m_{ij}^{l}|q_{ij}^l)\) is a PDHFE and the associated probabilities of the PDHFEs are unknown.

Step 2: For \(l=1,2,\ldots ,m\), Algorithm 1 is used to obtain the accurate associated probabilities of the PDHFEs and check the consistency of PDHFPRs \(D_{l}=(d_{ij}^{l})_{n\times n}\) . When \(\mathrm{{WCI}}(D_l)>\gamma _0\), the individual consistency can be improved using Algorithm 1. Thus, the improved PDHFPRs \(D_{l}^{*}(1\le l\le m)\) can be constructed, and the associated probabilities \(p_{ij,r}^l,q_{ij,k}^l\), the optimal variables \(\widetilde{o_{ij}^{l+}},\widetilde{o_{ij}^{l-}}, \widetilde{\varsigma _{ij}^{l+}},\widetilde{\varsigma _{ij}^{l-}}\) and the individual priority weights of alternatives \(\widetilde{w^{l}}=(\widetilde{w_1^{l}}, \widetilde{w_2^{l}},\ldots ,\widetilde{w_n^{l}})^T\) can be obtained, where \(\widetilde{w_i^l}=(\widetilde{w_{i\mu }^l}, \widetilde{w_{i\nu }^l})\).

Step 3: The following formula is used to calculate the weights \(\psi _l(l=1,2,\ldots ,m)\) of DMs:

$$\begin{aligned} \psi _l^{'}=1-\mathrm{{WCI}}(D_l),\psi _l=\psi _l^{'}/\sum _{l=1}^m\psi _l^{'}. \end{aligned}$$

Thus, \(0\le \psi _l\le 1\) and \( \sum _{l=1}^m\psi _l=1\).

Step 4: Calculate the overall priority weights: \(\varpi _i=\sum _{k=1}^m\psi _k\cdot \bigtriangleup (\widetilde{w_i^k})\), where \(\bigtriangleup (\widetilde{w_i^k})=\widetilde{w_{i\mu }^k} -\widetilde{w_{i\nu }^k}\).

Step 5: Obtain the ranking of \(x_i(i=1,2,\ldots ,n)\) using the value of \(\varpi _i\).

Assuming that the number of alternatives is n, the number of DMs participating in GDM is m, the maximum number of iterations is \(t_\mathrm{{max}}\), and the number of preference values given by the kth decision expert in the ith row and jth column of the decision-making matrix is \(\#d_{ij}^l\). The time complexity of Algorithm 2 is as follows:

  1. 1.

    The time complexity for the initialization process of the decision-making matrix is \(O(m\times n^2)\);

  2. 2.

    the time complexity for the calculation process of the associated probabilities of the PDHFEs is \(O(m\times n^2)\);

  3. 3.

    the time complexity for the testing process of the consistency is \(O(m\times n^2)\);

  4. 4.

    the time complexity for the improving process of the consistency level in the worst case is \(O(t_\mathrm{{max}}\times m\times n^2)\);

  5. 5.

    the time complexity for the calculation process of the weights of DMs is O(m);

  6. 6.

    the time complexity for the calculation process of the overall priority weights of alternatives is O(n);

    Thus, the total time complexity of of Algorithm 2 is \(O(t_\mathrm{{max}}\times m\times n^2)\).

The space complexity of Algorithm 2 is as follows:

The space required to store the basic parameters of Algorithm 2 is constant, the space required to store the initial preference values given by DMs is \((\sum _{l=1}^{m}\sum _{i}^n\sum _{j}^n(\#d_{ij}^l))+2n\), the space required to store the calculated probabilities is \((\sum _{l=1}^{m}\sum _{i}^n\sum _{j}^n(\#d_{ij}^l))-2n\), the space required to store the optimal variables is \(m\times n^2\), the space required to store the optimal normalized weight vector is 2n, the space required to store the improved preference values is \(2\times \sum _{l=1}^{m}\sum _{i}^n\sum _{j}^n(\#d_{ij}^l)\), the space required to store the weights of DMs is m.

To sum up, the total space complexity of Algorithm 2 is \(O(m\times n^2+\sum _{l=1}^{m}\sum _{i}^n\sum _{j}^n(\#d_{ij}^l))\).

Illustrative example

Application to analyze the impact factors of haze

Haze is produced from the interaction of human activities and specific climatic conditions. In recent years, severe haze weather has occurred in many regions of China. A large number of fine particles will inevitably be discharged because of the social and economic activities of high-density populations. When the emission of fine particulate pollutants is increasing and exceeds than the atmospheric circulation and carrying capacity, its concentration in the air will continue to increase, resulting in haze. The frequent occurrence of haze not only seriously affects air quality but is also the main cause of asthma and other diseases, which seriously affect people’s health and economic development. Therefore, it is of practical significance to study the causes of haze weather.

Shanghai has been harmed by haze weather for a long time. Three related experts \(P=\{p_1,p_2,p_3\}\) use PDHFPRs to express their preference for four evaluation factors of haze weather. The four evaluation factors of haze weather are PM2.5 (\(x_1\)), relative humidity (\(x_2\)),\(NO_2(x_3)\), and PM10 (\(x_4\)). Let the given consistency index threshold \(\gamma _0=0.2\), the weight coefficient \(\lambda =0.2\), and the adjusted coefficient \(\theta =0.5\). We use Algorithm 2 to explore the four evaluation factors of haze weather and analyze the leading causes of haze weather.

Step 1: Based on professional knowledge and experience, PDHFPRs \(D_l=(d_{ij}^l)_{4\times 4}=(h_{ij}^l|p_{ij}^l,m_{ij}^l| q_{ij}^l)_{4\times 4}(l=1,2,3)\) are given by experts and are shown in Tables 3, 4 and 5.

Step 2: Using the PDHFPR \(D_1\) given by expert \(P_1\), the following optimization model is established according to model (M-3):

$$\begin{aligned}&\min J_1 = o_{12}^{1-}+o_{12}^{1+}+\varsigma _{12}^{1-} +\varsigma _{12}^{1+}+o_{13}^{1-}+o_{13}^{1+}+\varsigma _{13}^{1-}+\varsigma _{13}^{1+}\nonumber \\&\qquad \quad +o_{14}^{1-}+o_{14}^{1+}+\varsigma _{14}^{1-} +\varsigma _{14}^{1+}+o_{23}^{1-}+o_{23}^{1+}+\varsigma _{23}^{1-}+\varsigma _{23}^{1+} \nonumber \\&\qquad \quad +o_{24}^{1-}+o_{24}^{1+}+\varsigma _{24}^{1-} +\varsigma _{24}^{1+}+o_{34}^{1-}+o_{34}^{1+}+\varsigma _{34}^{1-}+\varsigma _{34}^{1+}\nonumber \\&\mathrm{{s.t.}}\left\{ \begin{array}{l} (0.2\cdot p_{12,2}^1+0.3\cdot p_{12.2}^1+0.6\cdot p_{12,3}^1)\\ \quad +o_{12}^{1-}-o_{12}^{1+}=\sqrt{w_{1\mu }^1\cdot w_{2\nu }^1} \\ (0.1\cdot q_{12,2}^1+0.3\cdot q_{12.2}^1)+\varsigma _{12}^{1-}-\varsigma _{12}^{1+}=\sqrt{w_{1\nu }^1\cdot w_{2\mu }^1} \\ (0.4\cdot p_{13,1}^1+0.8\cdot p_{13.2}^1) +o_{13}^{1-}-o_{13}^{1+}=\sqrt{w_{1\mu }^1\cdot w_{3\nu }^1},\\ 0.37+\varsigma _{13}^{1-}-\varsigma _{13}^{1+}=\sqrt{w_{1\nu }^1\cdot w_{3\mu }^1} \\ (0.3\cdot p_{14,1}^1+0.5\cdot p_{14.2}^1) +o_{14}^{1-}-o_{14}^{1+}=\sqrt{w_{1\mu }^1\cdot w_{4\nu }^1} \\ (0.14\cdot q_{14,1}^1+0.24\cdot q_{14.1}^1)+\varsigma _{14}^{1-}-\varsigma _{14}^{1+}=\sqrt{w_{1\nu }^1\cdot w_{4\mu }^1} \\ 0.2+o_{23}^{1-}-o_{23}^{1+}=\sqrt{w_{2\mu }^1\cdot w_{3\nu }^1},\\ (0.6\cdot q_{23,1}^1+0.8\cdot q_{23.2}^1)+\varsigma _{23}^{1-}-\varsigma _{23}^{1+}=\sqrt{w_{2\nu }^1\cdot w_{3\mu }^1} \\ 0.25+o_{24}^{1-}-o_{24}^{1+}=\sqrt{w_{2\mu }^1\cdot w_{4\nu }^1},0.2\\ \quad +\varsigma _{24}^{1-}-\varsigma _{24}^{1+}=\sqrt{w_{2\nu }^1\cdot w_{4\mu }^1} \\ (0.2\cdot p_{34,1}^1+0.6\cdot p_{34.2}^1+0.7\cdot p_{34,3}^1)\\ \quad +o_{34}^{1-}-o_{34}^{1+}=\sqrt{w_{3\mu }^1\cdot w_{4\nu }^1} \\ (0.1\cdot q_{34,2}^1+0.2\cdot q_{34.2}^1)+\varsigma _{34}^{1-}-\varsigma _{34}^{1+}=\sqrt{w_{3\nu }^1\cdot w_{4\mu }^1} \\ p_{12,2}^1+p_{12.2}^1+p_{12,3}^1=1,q_{12,2}^1+q_{12,2}^1=1,\\ p_{12,2}^1+p_{12,3}^1=1,q_{14,1}^1+q_{14,2}^1=1, \\ p_{14,1}^1+p_{14.2}^1=1,q_{23,1}^1+q_{23,2}^1=1,\\ p_{34,2}^1+p_{34,3}^1=1,q_{34,1}^1+q_{34,2}^1=1, \\ 0\le p_{12,2}^1,p_{12.2}^1,p_{12,3}^1,q_{12,2}^1,\\ q_{12,2}^1,p_{12,2}^1,p_{12,3}^1,q_{14,1}^1,q_{14,2}^1 \\ p_{14,1}^1,p_{14.2}^1,q_{23,1}^1,q_{23,2}^1,p_{34,2}^1,\\ p_{34,3}^1,q_{34,1}^1,q_{34,2}^1\le 1 \\ o_{ij}^{1+},o_{ij}^{1-},\varsigma _{ij}^{1+},\varsigma _{ij}^{1-}\ge 0,\\ o_{ij}^{1+},o_{ij}^{1-}=0,\varsigma _{ij}^{1+}, \varsigma _{ij}^{1-}=0 \\ 0\le w_{1\mu }^1+w_{1\nu }^1\le 1,\sum _{j=1,j\ne i}^4w_{j\mu }^1\le w_{i\nu }^1,\\ \displaystyle \sum _{j=1,j\ne i}^4w_{j\nu }^1\le w_{i\mu }^1+2\; i<j,i,j=1,2,3,4 \end{array}\right. \end{aligned}$$
(25)

By using the MATLAB Toolbox, we obtain

$$\begin{aligned}&p_{12,1}^1=0.0969,p_{12,2}^1=0.1187,p_{12,3}^1=0.7844, \\&q_{12,1}^1=0.5996,q_{12,2}^1=0.4004,p_{13,1}^1=1,p_{13,2}^1=0,\\&p_{14,1}^1=0,p_{14,2}^1=1,q_{14,1}^1=1,q_{14,2}^1=0,q_{23,1}^1=1,\\&q_{23,2}^1=0,p_{34,1}^1=0.3181,p_{34,2}^1=0.2606,\\&p_{34,3}^1=0.4213,q_{34,1}^1=0.4780,q_{34,2}^1=0.5220,\\&o_{12}^{1-}=o_{12}^{1+}=\varsigma _{12}^{1-} =\varsigma _{12}^{1+}=o_{13}^{1-}=o_{13}^{1+}=\varsigma _{13}^{1+}\\&=o_{14}^{1-}=o_{14}^{1+}=\varsigma _{14}^{1-}=o_{23}^{1-}=o_{23}^{1+}= \varsigma _{23}^{1-}=o_{24}^{1-}\\&=o_{24}^{1+}=\varsigma _{24}^{1-}=\varsigma _{24}^{1+} =o_{34}^{1-}=o_{34}^{1+}=\varsigma _{34}^{1-}=\varsigma _{34}^{1+}=0,\\&\varsigma _{13}^{1-}=0.0009,\varsigma _{14}^{1+}=0.0030,\varsigma _{23}^{1+} =0.0587,\\&\widetilde{w_1^1}=\langle 0.2986,0.4345\rangle ,\\&\widetilde{w_2^1}=\langle 0.0746,0.9254\rangle ,\widetilde{w_3^1} =\langle 0.3166,0.5359\rangle ,\\&\widetilde{w_4^1}=\langle 0.0432,0.8373\rangle . \end{aligned}$$
Table 3 PDHFPR matrix provided by expert
Table 4 PDHFPR matrix provided by expert
Table 5 PDHFPR matrix provided by expert

Therefore, \(\mathrm{{WCI}}(D_1)=0.0490<\gamma _0\) , which implies that \(D_1\) is an acceptable consistent PDHFPR. Similarly, using PDHFPR \(D_2\) given by expert \(p_2\), an optimization model is built in line with model (M-3), and the solution is as follows:

$$\begin{aligned}&p_{13,1}^2=0.8959,p_{13,2}^2=0,p_{13,3}^2=0.1041,q_{23,1}^2=0,q_{23,2}^2=1,\\&p_{24,1}^2=0.9589,p_{24,2}^2=0.0411,p_{34,1}^1=1,p_{34,2}^1=0,\\&\varsigma _{12}^{2+}=0.0786,\varsigma _{13}^{2-}=0.0844,o_{14}^{2+}=0.0665,\\&\varsigma _{24}^{2-}=0.0399,o_{34}^{2+}=0.072,\varsigma _{34}^{2-}=0.0258,\\&\widetilde{w_1^2}=\langle 0.4373,0.4866\rangle , \widetilde{w_2^2}=\langle 0.1008,0.8232\rangle ,\\&\widetilde{w_3^2}=\langle 0.3037,0.6203\rangle ,\widetilde{w_4^2} =\langle 0.0822,0.9178\rangle ,\\&o_{12}^{2-}=o_{12}^{2+}=\varsigma _{12}^{2-} =o_{13}^{2-}=o_{13}^{2+}=\varsigma _{13}^{2+}\\&\quad \quad =o_{14}^{2-}=\varsigma _{14}^{2-}=\varsigma _{14}^{2+} =o_{23}^{2-}=o_{23}^{2+}=\varsigma _{23}^{2+}\\&\quad \quad =\varsigma _{23}^{2-}=o_{24}^{2-}=o_{24}^{2+} =\varsigma _{24}^{2+}=o_{34}^{2-}=\varsigma _{34}^{2+}=0. \end{aligned}$$

Therefore, \(\mathrm{{WCI}}(D_2)=0.0905<\gamma _0\) , which implies that \(D_2\) is an acceptable consistent PDHFPR. Similarly, using PDHFPR \(D_3\) given by expert \(p_3\), an optimization model is built and we can obtain \(\mathrm{{WCI}}(D_3)=0.3623>\gamma _0\), which implies that \(D_3\) is an unacceptable consistent PDHFPR. For this case, Algorithm 1 can be applied to improve the consistency of \(D_3\).

Using Eq. (23) to determine where expert \(e_3\) needs to modify preference values, we can obtain \((u,m)=(2,3)\). The improved PDHFE \(d_{23}^{3*}=\{\{0.46|p_{23,1}^3,0.56|p_{23,2}^3\}, \{0.27|1\}\}\) is obtained according to Eq. (24). Thus, an improved PDHFPR \(D_3^{*}\) can be constructed by introducing \(d_{23}^{3*}\) to \(D_3\). Furthermore, according to model (M-3), the results are obtained as follows:

$$\begin{aligned}&q_{12,1}^3=0,q_{12,2}^3=1,p_{13,1}^3=0.1719,p_{13,2}^3=0.8281,\\&p_{14,1}^3=0,p_{14,2}^3=1,p_{23,1}^3=1,p_{23,2}^3=0,\\&q_{24,1}^3=0.5994,q_{24,2}^3=0.4006,q_{34,1}^3=0,q_{34,2}^3=1,\\&o_{12}^{3-}=\varsigma _{12}^{3+}=o_{13}^{3-}=o_{13}^{3+} =\varsigma _{13}^{3+}=\varsigma _{13}^{3-}\\&\quad =o_{14}^{3-}=o_{14}^{3+}=\varsigma _{14}^{3-} =o_{23}^{3-}=\varsigma _{23}^{3+}=o_{24}^{3-}\\&\quad =o_{24}^{3+}=\varsigma _{24}^{3+}=\varsigma _{24}^{3-} =o_{34}^{3-}=\varsigma _{34}^{3+}=\varsigma _{34}^{3-}=0\\&o_{12}^{3+}=0.0688,\varsigma _{12}^{3-}=0.0137,\varsigma _{14}^{3+}=0.1102,\\&o_{23}^{3+}=0.1381,\varsigma _{23}^{3-}=0.0687,o_{34}^{3+}=0.2131,\\&\widetilde{w_1^3}=\langle 0.3308,0.6692\rangle ,\widetilde{w_2^3} =\langle 0.1470,0.8530\rangle ,\\&\widetilde{w_3^3}=\langle 0.1345,0.7048\rangle ,\widetilde{w_4^3} =\langle 0.2270,0.6122\rangle . \end{aligned}$$

We can obtain \(\mathrm{{WCI}}(D_3^{*})=0.1909<\gamma _0\), which implies that \(D_3^{*}\) is an acceptable consistent PDHFPR.

Step 3: The individual weights of DMs can be calculated as follows: \(\psi _1=0.3562,\psi _2=0.3407,\psi _3=0.3031.\)

Step 4: Compute overall priority weights: \(\varpi _1=-0.1677\), \(\varpi _2=-0.7631\), \(\varpi _3=-0.3588\), \(\varpi _4=-0.6843\).

Step 5: As \(\varpi _1\succ \varpi _3\succ \varpi _4\succ \varpi _2\), the four evaluation factors of haze weather are ranked as

$$\begin{aligned} x_1\succ x_3\succ x_4\succ x_2. \end{aligned}$$

Therefore, PM2.5 is the most critical cause of haze in Shanghai. This result agrees with the calculation results of this method.

Sensitivity analysis of \(\lambda \) and \(\theta \)

For \(0\le \lambda \le 1(\theta =0.5)\), we take 0.1 as the interval and conduct 11 experiments in total. The influence of \(0\le \lambda \le 1(\theta =0.5)\) on DMs’ weights and ranking results of four evaluation factors are shown in Table 6. Figure 1 shows the influence of \(0\le \lambda \le 1(\theta =0.5)\) on the overall priority weights of alternatives \(\varpi _i=(i=1,2,3,4)\).

Table 6 DMs weights and ranking results of four evaluation factors with different \(\lambda \) and \(\theta \)
Fig. 1
figure 1

Results of \(\varpi _i\) of alternatives with different \(\lambda \)

Table 6 shows that DMs weights change with \(0\le \lambda \le 1(\theta =0.5)\). Furthermore, when \(0\le \lambda \le 1\) takes different values, the ranking of haze impact factors changes. When \(\lambda \le 0.1\), the ranking result is \(x_1\succ x_3\succ x_2\succ x_4\). When \(0.1<\lambda \le 1\), the ranking is \(x_1\succ x_3\succ x_4\succ x_2\). Therefore, the weight coefficient of \(\mathrm{{OCI}}(D)\) affects the ranking results, which proves the flexibility of the decision approach in this paper. However, regardless of how \(\lambda \) changes, the most important factor affecting haze is still PM 2.5. Therefore, our proposed approach exhibits good stability.

Moreover, we study the ranking of haze impact factors when \(0\le \theta \le 1(\lambda =0.2)\) changes from 0.1 to 1. The specific ranking results are shown in Table 6, and the values of overall priority weights \(\varpi _i(i=1,2,3,4)\) are shown in Fig 2. The ranking results do not change with \(\theta \) because \(\theta \) is the consistency adjusted coefficient, which only affects the speed of \(D_3\) reaching the consistency threshold. Specifically, the larger the adjusted coefficient \(\theta \) is, the faster \(D_3\) reaches the consistency threshold.

Comparative analysis

To prove the superiority of our approach, we will use the same illustrative example to compare our method with the methods proposed in Hao et al. [17], Garg and Kaur [40, 42] and Zhao et al. [43]. Under PDHFS environment, the method in Hao et al. [17] and Garg and Kaur [40, 42] requires all PDHFEs to be complete. When applying the methods in [17, 40] to address the illustrative example, the missing associated probabilities of PDHFEs are derived from the associated probability calculation method in this paper.

Fig. 2
figure 2

Results of \(\varpi _i\) of alternatives with different \(\theta \)

(1) Hao et al. [17] defined basic operators for PDHFSs and then proposed a visualization GDM method with PDHFSs based on the operator. For the sake of comparison, the individual weights of DMs is set as the calculation result of the method in this paper, that is, \(\psi _1=0.3562,\psi _2=0.3407,\psi _3=0.3031\) Now, we use the method in Hao et al. [17] to analyze the influencing factors of haze weather.

First, the comprehensive PDHFPR matrix \(D^{*}=(d_{ij}^{*})_{n\times n}\) is obtained using the probabilistic dual hesitant fuzzy weighted average (PDHFWA) operator (Eq. (19) in [17]), where the following elements of the comprehensive matrix are shown:

$$\begin{aligned} d_{ii}^{*}= & {} (\{0.5|1\},\{0.5|1\}),\\ d_{12}^{*}= & {} (\{0.488|0.0969,0.5118|0.1187,0.6|0.7844\},\\&\{0.2028|0.5996,0.3|0.4004\}),\\ d_{13}^{*}= & {} (\{0.4361|0.154,0.4664|\\&\quad 0.7419,0.5262|0.0179,0.5517|0.0862\},\{0.3233|1\}), \\ d_{14}^{*}= & {} (\{0.5676|1\},\{0.2325|1\}) \\ d_{21}^{*}= & {} (\{0.2344|0.5996,0.3|0.4004\},\\&\{0.4057|0.0969,0.4687|0.1187,0.6|0.7844\}), \\ d_{23}^{*}= & {} (\{0.3657|1\},\{0.4042|1\}),\\ d_{24}^{*}= & {} (\{0.2826|0.9589,0.3193|0.0411\},\\&\{0.2833|0.5994,0.3031|0.4006\}), \\ d_{31}^{*}= & {} (\{0.3258|1\},\{0.4316|0.154,0.4617|\\&\quad 0.7419,0.484|0.0179,0.5179|0.0862\}), \\ d_{32}^{*}= & {} (\{0.4675|1\},\{0.3031|1\}),\\ d_{34}^{*}= & {} (\{0.4522|0.3181,0.572|0.2606,0.6137|0.4213\},\\&\{0.1928|0.478,0.2468|0.522\}), \\ d_{41}^{*}= & {} (\{0.2881|1\},\{0.5431|1\}),\\ d_{42}^{*}= & {} (\{0.2994|0.5994,0.3371|0.4006\},\\&\{0.2811|0.9589,0.31|0.0411\}), \\ d_{43}^{*}= & {} (\{0.2354|0.478,0.2668|0.522\},\\&\{0.3839|0.3181,0.5677|0.2606,0.5998|0.4213\}). \end{aligned}$$

Second, the preference values \(d_{ij}^{*}\) are aggregated to obtain the aggregated probabilistic dual hesitant fuzzy preference values \(d_{i}^{*}\) of the factor affecting haze \(x_i\). Due to the length limitation of this paper, only some elements of \(d_{i}^{*}\) are shown as follows:

$$\begin{aligned} d_1^{*}= & {} (\{0.5001|0.0149,\ldots ,0.5501|0.014\},\\&\{0.2955|0.5996,\ldots ,0.3259|0.4004\}), \\ d_2^{*}= & {} (\{0.354|0.5749,\ldots ,0.3765|0.0165\},\\&\{0.3904|0.0581,\ldots ,0.4379|0.3142\}), \\ d_3^{*}= & {} (\{0.440|0.3181,\ldots ,0.4868|0.4213\},\\&\{0.3345|0.0736,\ldots ,0.3724|0.045\}), \\ d_4^{*}= & {} (\{0.3391|0.2865,\ldots ,0.3551|0.209\},\\&\{0.4138|0.305,\ldots ,0.4741|0.0173\}). \end{aligned}$$

Next, the score function of \(d_i^{*}\) is calculated as follows:

$$\begin{aligned}&s(d_1)=0.2242,s(d_2)=-0.0664,\\&s(d_3)=0.1169,s(d_4)=-0.0995. \end{aligned}$$

Finally, the ranking of the factors affecting haze is as follows:

$$\begin{aligned} x_1\succ x_3\succ x_2\succ x_4. \end{aligned}$$

(2) Garg and Kaur [40] defined the correlation coefficients on PDHFSs and proposed the multi-criteria decision-making (MCDM) method based on the correlation coefficients for PDHFSs. Now, we use the method in Garg and Kaur [40] to handle the above problem, and the main steps are as follows:

First, by using Algorithm 1 in Garg and Kaur [40], the aggregated decision matrix \(D^{*}=(d_{ij}^{*})_{n\times n}\) is obtained. The MDs and NMDs and their corresponding probabilities are listed as follows:

$$\begin{aligned} d_{11}^{*}= & {} (\{0.5|1\},\{0.5|1\}),\\ d_{12}^{*}= & {} (\{0.2|0.0323,0.3|0.0396,0.6|0.9281\},\\&\{0.1|0.1999,0.3|0.8001\}), \\ d_{13}^{*}= & {} \left( \{0.4|0.3906,0.5|0.5747,0.7|0.0347\},\right. \\&\left. \left\{ 0.3|\frac{2}{3},0.37|\frac{1}{3}\right\} \right) , \\ d_{14}^{*}= & {} \left( \left\{ 0.45|\frac{1}{3},0.5|\frac{1}{3},0.7|\frac{1}{3}\right\} , \left\{ 0.14|\frac{1}{3},0.2|\frac{1}{3},0.5|\frac{1}{3}\right\} \right) , \\ d_{23}^{*}= & {} \left( \left\{ 0.2|\frac{1}{3},0.25|\frac{1}{3},0.6|\frac{1}{3}\right\} , \left\{ 0.2|\frac{1}{3},0.5|\frac{1}{3},0.6|\frac{1}{3}\right\} \right) , \\ d_{24}^{*}= & {} \left( \left\{ 0.25|\frac{1}{3},0.3|0.6530,0.4|0.0136\right\} \right. \\&\left. ,\left\{ 0.2|\frac{1}{3},0.3|\frac{1}{3},0.4|0.1998,0.5|0.1335\right\} \right) , \\ d_{34}^{*}= & {} (\{0.2|0.1060,0.5|0.3334,0.6|0.4202,0.7|0.1404\}, \\&\{0.1|0.1593,0.2|0.5073,0.4|0.3333\}). \end{aligned}$$

Second, the ideal alternative \(d^{*}\) is determined by Eq. (15) in [40]as follows:

$$\begin{aligned} d^{*}=\left. \begin{array}{l} \{< x_1, (\{0.5|1\}, \\ \{0.4|0.3906,0.5|0.5747,0.7|0.0347\})>,\\< x_2, (\{0.2|0.0323,0.3|0.0396,0.6|0.9281\},\\ \{0.1|0.1999,0.3|0.8001\})>\},\\< x_3, (\{0.5|1\}, \{0.3|\frac{2}{3},0.37|\frac{1}{3}\})>\}, \\ < x_4, (\{0.45|\frac{1}{3},0.5|\frac{1}{3},0.7|\frac{1}{3}\},\\ \{0.1|0.1593,0.2|0.5073,0.4|0.3333\})>\}, \end{array}\right. \end{aligned}$$

Then, the correlation coefficient index \(K_1,K_2,K_3,K_4\), (Eq. (10), (11), (12), (13)) in [40]) is used to obtain the measurement between \(d_i^{*}(i=1,2,3,4)\) and \(d^{*}\), and we getthe following:

$$\begin{aligned}&K_1(d_1^{*},d^{*})=0.8528, K_1(d_2,d^{*})=0.7087, \\&K_1(d_3,d^{*})=0.7312, K_1(d_4,d^{*})=0.7838,K_2(d_1^{*},d^{*})=0.8478, \\&K_2(d_2,d^{*})=0.6729, K_2(d_3,d^{*})=0.5961, K_2(d_4,d^{*})=0.6201. \end{aligned}$$

The selection of \(\omega \) in [40] is subjective. When \(\omega =(0.25,0.25,0.25,0.25)^T.\) we get the following:

$$\begin{aligned} K_3(d_i^{*},d^{*})= & {} K_1(d_i,d^{*}),K_4(d_i,d^{*})\\= & {} K_2(d_i,d^{*}).(i=1,2,3,4). \end{aligned}$$

Finally, the ranking of four evaluation factors is \(x_1\succ x_4\succ x_2\succ x_3\) when the correlation coefficient index \(K_1\) and \(K_3\) are used and is \(x_1\succ x_2\succ x_4\succ x_3\) when the correlation coefficient index \(K_2,K_4\) are used.

(3) Garg and Kaur [42] studied some distance measures between PDHFSs and distance measure-based maximum deviation method were set to calculate the weights of criteria. Then, a new MCDM method with unknown attribute weights was proposed. By using the method in Garg and Kaur [42], the following steps are involved:

First, a comprehensive matrix is obtained by aggregating the individual PDHFPRs using Algorithm 1 in Garg and Kaur [42], and the elements of the comprehensive matrix are the same as Garg and Kaur’s [40] method.

Second, because the number of elements in MDs and NMDs of each PDHFE is different, the distance measure \(d_2\) (Eq. (10) in [42]) is used to obtain the following weights:

$$\begin{aligned} \omega =(0.2326,0.262,0.2427,0.2627)^T. \end{aligned}$$

Then, the probabilistic dual hesitant fuzzy weighted Einstein average

(PDHFWEA) operator (Eq. (31) in [42]) is used to obtain the aggregate values of each alternative. Due to the length limitation of this paper, only some elements of the PDHFPR matrix are shown as follows:

$$\begin{aligned} d_1^{*}= & {} (\{0.5054|0.1208,\ldots ,0.4767|0.0004\},\\&\{0.2962|0.0889,\ldots ,0.2995|0.0444\}), \\ d_2^{*}= & {} (\{0.3205|0.0156,\ldots ,0.4192|0.0009\},\\&\{0.4393|0.1031,\ldots ,0.3302|0.0014\}), \\ d_3^{*}= & {} (\{0.5602|0.0156,\ldots ,0.3823|0.0741\},\\&\{0.2981|0.1095,\ldots ,0.3937|0.0018\}), \\ d_4^{*}= & {} (\{0.2724|0.1208,\ldots ,0.4139|0.0072\},\\&\{0.4615|0.0156,\ldots ,0.4698|0.0396\}). \end{aligned}$$

Similarly, the probabilistic dual hesitant fuzzy ordered weighted Einstein average (PDHFOWEA) operator (Eq. (32) in [42]), probabilistic dual hesitant fuzzy weighted Einstein geometric (PDHFWEG) operator (Eq. (33) in [42]), and probabilistic dual hesitant fuzzy ordered weighted Einstein geometric (PDHFOWEG) operator (Eq. (34) in [42]) operators are used to calculate the aggregated values for each alternative, and the score function values and ranking results using the four types of operators are given in Table 7.

Table 7 shows that the most crucial indicator affecting the generation of haze in Shanghai is PM 2.5, which matches the result of this paper.

(4) When the probability values of the MDs and the NMDs are equal, the PDHFPRs degenerates into the dual hesitant fuzzy PRs. With respect to the GDM problem with dual hesitant fuzzy PRs, Zhao et al. [43] developed new GDM methods to obtain the best alternative, which includes a consensus reaching process. The method 1 in Zhao et al. [43] is utilized as the comparison algorithm.

The dual hesitant fuzzy PRs corresponding to PDHFPRs \(D_l=(d_{ij}^l)_{4\times 4}=(h_{ij}^l|p_{ij}^l, m_{ij}^l|q_{ij}^l)_{4\times 4}(l=1,2,3)\) are recorded as \(c_l=(h_{ij}^l,m_{ij}^l)_{4\times 4}\), and the calculation process is as follows:

First, the compatibility degrees \(C_{kl}=CO_1(C_k,C_l),k,l=1,2,3\) between \(C_k\) and \(C_l\) are computed by \(CO_1\) (Eq. (5) in [43]), which are given as follows:

$$\begin{aligned} C_1=(C_{kl})_{3\times 3}=\left( \begin{array}{ccc} 1 &{}\quad 0.9073 &{}\quad 0.8400 \\ 0.9073 &{}\quad 1 &{}\quad 0.9040 \\ 0.8400 &{}\quad 0.9040 &{}\quad 1 \\ \end{array}\right) \end{aligned}$$

Therefore, the weights of experts are determined by Eq. (9) in [43]:\(\psi _1^1=0.3295,\psi _2^1=3416,\psi _3^1=0.3289\).

Table 7 Comparison of different operators
Table 8 Decision-making results by different methods

Second, for all \(l=1,2,3\), let \(C_l^0=C_l\). Then, using the symmetric dual hesitant fuzzy weighted averaging operator (Eq. (10) in [43]), all dual hesitant fuzzy preference relations \(C_1^0\) are fused into \(C^0=(c_{ij0})_{4\star 4}\), where \(c_{ij0}=(h_{ij0},g_{ij0}),(i,j=1,2,3,4)\) . Due to the length limitation of this paper, some elements in \(C^0\) are shown as follows:

$$\begin{aligned} h_{120}= & {} \{0.6000,0.4982,0.4539\},g_{120}\\= & {} \{0.2829,0.3000,0.2019,0.2155\} \\ h_{130}= & {} \left. \begin{array}{c} \{0.4337,0.4667,0.4679,0.5012,0.5056,0.5389, \\ 0.5802,0.6123,0.6135,0.6446,0.6486,0.6784,\},\\ \end{array}\right. \\ g_{130}= & {} \{0.3222\}. \end{aligned}$$

Then, the compatibility degrees \(CO_1(C^0,C_l^0)\) between \(C_l^0, (l=1,2,3)\) and \(C^0\) are calculated using Eq. (5) in [43] as follows:

$$\begin{aligned} CO_1(C_1^0,C^0)= & {} 0.935,CO_1(C_2^0,C^0)\\= & {} 0.9793,CO_1(C_3^0,C^0)=0.9443. \end{aligned}$$

Thus, the group reaches the given threshold.

Using the symmetric dual hesitant fuzzy averaging operator (Eq. (13) in [43]), we integrate the dual hesitant fuzzy preference values of each alternative \(x_i(i=1,2,3,4)\) into the collective dual hesitant fuzzy preferences \(d_i(i=1,2,3,4)\). Then, using Definition 3 in Zhao et al. [43], the score values are as follows:

$$\begin{aligned} s(d_1)= & {} -0.1813,s(d_2)=-0.0938,s(d_3)\\= & {} -0.0568,s(d_4)=-0.1319. \end{aligned}$$

Thus, the ranking of the four evaluation factors is as follows:

$$\begin{aligned} x_1\succ x_3\succ x_2\succ x_4. \end{aligned}$$

In the above process, if we measure the agreement of preferences by compatibility measure (Eq. (6) in [43]) instead of (Eq. (5) in [43]), the ranking of the factors affecting haze is \(x_1\succ x_3\succ x_2\succ x_4\).

Based on the same example of analyzing the impact factors of haze weather, the calculation results of the above four methods are shown in Table 8.

Table 9 Comparison of the proposed approach and other approaches

Table 9 visually demonstrates the comparison of different methods, and the following detailed comparative analysis shows the superiority of the proposed GDM method:

  1. 1.

    By comparing the decision results in Table 8, it is shown that the most crucial factor affecting haze in Shanghai is PM 2.5, which proves the reasonable effectiveness of our method.

  2. 2.

    Our method combines the DMs’ risk attitude when calculating the consistency level for PDHFPR. When the weight coefficient \(\lambda \le 0.1\), the ranking result is \(x_1\succ x_3\succ x_2\succ x_4\), which is the same as the ranking result by Hao et al. [17], Garg and Kaur [42] and Zhao et al. [43]. And when \(0.1<\lambda \le 1\), the ranking result is \(x_1\succ x_3\succ x_4\succ x_2\). Therefore, the DMs’ risk attitude affects the ranking results, which proves the flexibility of our method.

  3. 3.

    Comparison with the methods in Hao et al. [17] and Garg and Kaur [42], all the methods solve GDM problem with probabilistic dual hesitation fuzzy environment. However, Hao et al. [17] and Garg and Kaur [42] used PDHFSs to directly evaluate alternatives, while the PDHFPRs proposed in this paper only require the DMs to give evaluation information through pairwise comparison, so PDHFPRs are more appropriate for the fuzzy decision of DMs in complicated environments. In addition, Hao et al.s [17] and Garg and Kaurs [42] method can only handle GDM problems with known probability information of PDHFEs. In contrast, our method comprehensively considers the situation of incomplete decision making information and unknown weights of DMs, which is more suitable for the increasingly complicated decision environment.

  4. 4.

    From Table 8, although Garg and Kaurs [42] method selects the same most crucial factor affecting haze with our method, there is a different ranking between the four factors affecting haze in the method of Garg and Kaur [42] and our method. This is mainly because Garg and Kaurs [42] method used the correlation coefficients for PDHFSs to get the ranking result while our method uses the optimization model of PDHFPRs to directly obtain the ranking of alternatives and needs to calculate the overall priority weights of the four evaluation factors by improving the consistency level for the individual PDHFPRs. Therefore, our method can obtain more credible decision results.

  5. 5.

    The comparison with the method in Zhao et al. [43] is as follows: First, Zhao et al.’s [43] method neglects the associated probabilities of MDs and NMDs, which may lose and distort the original information. Our method uses the original PDHFPR information for calculation, which can better retain the original information. Second, our method requires higher consistency than Zhao et al. [43]. In our method, the original PDHFPR of decision expert \(p_3\) is unacceptable consistent while the consistency level of decision expert \(p_3\) is not required to be improved in Zhao et al. [43]. Finally, our approach only adjusts one preference value each time, so the adjusted PDHFPR matrix is closer to the original PDHFPRs, which can better retain the wishes of DMs and greatly improve the accuracy of decision-making results.

Based on the above discussion and Table 9, our method compensates for the deficiencies of the existing methods and offer good application value.

Conclusion and future work

In this paper, the concept of PDHFPRs has been first given, and the multiplicative consistency, order consistency and acceptable consistency of PDHFPRs have been proposed. Considering that decision making experts cannot give accurate probabilities of the elements in PDHFPRs, based on multiplicative consistency, an optimization model that minimizes the deviation variable as the objective function has been constructed. Then, the probability calculation method has been given. According to the risk attitude of DMs, the consistency level for PDHFPRs has been tested using the weighted consistency index. For unacceptable consistency PDHFPRs, the obtained optimal deviations have been used to improve the PDHFE with the largest deviation one by one to adjust the consistency level for the original PDHFPRs. Based on the improved acceptable preference relation, the GDM method of PDHFPRs has been set, and the overall priority weights have been obtained. Finally, by analyzing the impact factors of haze weather, the validity of our method has been illustrated.

Although the method proposed in this paper has many advantages, it also has certain disadvantages. On the one hand, this paper presents the weight calculation methods for DMs, but they are not accurate enough, On the other hand, the given consistency threshold in this paper is determined by the DMs, which has a certain degree of subjectivity. Future research could focus on how to give a reasonable calculation method of accurate weights for DMs as well as an objective and appropriate consistency threshold. In addition, due to the different knowledge and experience of DMs, DMs may use heterogeneity preference relations to express their decision information. Therefore, it is interesting to further study the GDM approach with heterogeneity PDHFPRs.