Risk assessment based on a new decision-making approach with fermatean fuzzy sets

Plastic injection moulding machine

Plastic injection moulding machines contain many hazards and risks. Plastic injection moulding machines that are potentially hazardous demand careful handling to ensure OHS. Therefore, it is important for individuals working with injection moulding machines to understand these potential risks and adopt OHS measures. In the literature, studies on this machine generally focus on industrial processes such as supplier and raw material selection, factors affecting plastic production, and evaluation of maintenance strategies (Duman, Kongar & Gupta, 2019; Merizalde-Salas, Zumba-Novay & Peralta-Zurita, 2023; Gul, Yucesan & Celik, 2020; Pérez-Domínguez et al., 2018). However, there is a noticeable absence of studies focusing on the assessment of hazards and risks related to OHS through the utilisation of MCDM methods specifically in the context of plastic injection moulding machines.

Risk parameters in risk assessment

In risk assessment methods, depending on the method used, probability, severity, frequency and detectability parameters are used. When the studies in the literature are examined, a number of different risk parameters have been used other than the risk parameters used in classical risk analysis methods. These risk parameters are undetectability (Badida, Janakiraman & Jayaprakash, 2023), cost (Yazdi et al., 2020; Fattahi et al., 2020), time (Yazdi et al., 2020; Fattahi et al., 2020; Bhattacharjee, Dey & Mandal, 2022; Nabizadeh et al., 2021; Ghoushchi et al., 2021), profit (Yazdi et al., 2020; Fattahi et al., 2020), quality (Nabizadeh et al., 2021), operator competency (Bhattacharjee, Dey & Mandal, 2022), pollution ratio (Ghoushchi et al., 2021), sensitivity to personal protective equipment non-utilisation (Gul, Ak & Guneri, 2017; Liu et al., 2021), sensitivity to maintenance non-execution (Gul, Ak & Guneri, 2017), applicability of preventive measures (Gul & Ak, 2022), system reliability (Yener & Can, 2021), prevention action opportunity (Gul & Ak, 2021), effectiveness (Gul & Ak, 2021), maintenance cost (Valipour et al., 2022), maintenance duration (Valipour et al., 2022) and human error probability (La Fata et al., 2021). The use of more risk parameters helps to make more accurate decisions in decision-making processes and to assess risks more effectively. In the literature, it is seen that many different risk parameters are used in addition to the risk parameters used in classical risk assessment methods. This is useful in terms of a better understanding of the different aspects and impacts of risk and the development of an all-round risk assessment strategy.

In this study, the probability, severity, frequency, detectability, human error, machine error and existing safety measure are used as risk parameters. There is not any risk assessment study that includes the parameters of machine error, human error and existing safety measures in the existing literature. The use of these new risk parameters emphasises the unique aspect of this study in the literature.

Fermatean fuzzy sets

In decision-making processes, MCDM methods and fuzzy logic have been combined and used effectively in solving problems in various fields (Wan, Dong & Chen, 2024a, 2024b). In MCDM methods, fuzzy numbers are used to express the ambiguity of abstract information in a more flexible, understandable and realistic way during evaluations. In this context, FFSs that handle fuzzy logic have gained an important place in the literature and offer the opportunity to model uncertainty more flexibly.

There has been some exploration into their properties, building upon the established understanding of FFSs (Liu, Liu & Chen, 2019). Table 1 consolidates information on various studies employing FFSs, along with insights into the application of MCDM methods. The extended literature review is in the File S1. Despite the potential application areas identified for FFSs in existing literature, the utilisation of IVFFSs in comparable domains is yet to be extensively investigated.

When the outputs obtained as a result of the literature review are analysed, it is seen that FFSs have been applied to different MCDM problems. However, the number of studies using IVFF numbers is quite low. According to this literature review, there is no risk assessment study based on the integration of IVFF numbers with the AHP-TOPSIS method. The integration of both methods aims to enhance the approach and provide a more systematic and accurate risk assessment. When the application areas are examined, it is seen that FFSs are used in different areas such as the energy sector, supply chain management, environmental management, construction, logistics and the health sector. For example, numerous studies have been conducted on different topics such as the evaluation of renewable energy sources (Göçer, 2024; Seikh & Chatterjee, 2024), the selection of electric vehicles and charging stations (Golui, Mahapatra & Mahapatra, 2024), and waste management (Debbarma, Chakraborty & Saha, 2024; Seikh & Chatterjee, 2024).

As a result, this study aims to improve the decision-making process using AHP and TOPSIS methods based on the field of plastic injection moulding machine risk assessment.

The AHP method in risk assessment methods

The AHP method is one of the most widely used MCDM methods for prioritising criteria (Saaty, 1994). This method offers the opportunity to evaluate both qualitative and quantitative criteria together by including the priorities of the decision maker in the decision process (Ayyildiz & Taskin Gumus, 2020). It simplifies complex decision-making problems by creating a hierarchical structure (Gul, 2018). In addition, decision-makers make a pairwise comparison of criteria in the AHP method and determine the importance of these criteria (Bakioglu & Atahan, 2021; Biswas et al., 2025). One of the most important advantages of using the AHP method is that it allows the consistency of the decision-maker’s evaluations to be measured, and the group has a structure suitable for decision-making. AHP and AHP’s fuzzy versions are widely used decision-making methods, especially in risk assessment processes (Gul, 2018). They offer an objective and systematic assessment process. The AHP and AHP’s fuzzy extensions used in risk assessment studies in the literature are widely used in the weighting of risk parameters (Gul & Guneri, 2016; Ayvaz et al., 2024) and in the assessment of hazards (Aminbakhsh, Gunduz & Sonmez, 2013; Othman et al., 2016). In this study, the AHP method is used to determine the weight of risk parameters and to obtain hazard scores.

The TOPSIS method in risk assessment methods

TOPSIS is one of the MCDM methods. This method is used effectively, especially in solving complex decision problems (Biscaia, Junior & Colmenero, 2021; Hwang & Yoon, 1981). The TOPSIS method calculates the distances of alternatives from positive ideal and negative ideal solutions and offers the opportunity to rank them according to their proximity to the best solution (Hezer, Gelmez & Özceylan, 2021; Kizielewicz & Sałabun, 2024). Thanks to its distance-based evaluation capability, it provides an effective approach in cases where criteria are not of equal importance and delicate balances need to be observed. It also has the advantage of being able to process a large number of alternatives and criteria without affecting the calculation time (Jati, 2012). The TOPSIS method and its fuzzy extensions are widely used in risk assessment processes. The TOPSIS method has found wide application in the field of risk assessment in the literature (Ali & Maryam, 2014; Jozi & Majd, 2014). In this study, to rank the hazards of plastic injection moulding machines and to determine the risk coefficient value, the TOPSIS method is used.

Research motivation

Although current decision-making methods offer significant advantages in evaluating different criteria and alternatives, there is a need for approaches that can systematically and flexibly integrate the views of decision-makers in complex problems. The main motivation of this study is to develop a new decision-making framework for the separation of an inclusive decision-making group into SDMGs according to evaluation criteria. This new decision-making framework is applied to the risk assessment process, where many different criteria are evaluated. Due to the gap in the literature, the risk assessment of the hazards of the plastic injection moulding machine is carried out. New risk parameters are used in this risk assessment process to go beyond traditional approaches to risk parameters and make the risk assessment more detailed and comprehensive. This allows a multidimensional assessment of the impact of risk parameters on risk value and offers a new perspective. This study presents a method based on the integration of AHP and TOPSIS methods under IVFFS for risk assessment. AHP is preferred in this study because it provides a systematic analysis of the decision-making process, uses pairwise comparison matrices and provides consistency control of these matrices. TOPSIS is preferred because of the practicality and ease of application in transferring the weights and scores obtained by the AHP method directly in the decision matrix. With the integration of these two methods, it is ensured that the analysis process is shortened and has a dynamic structure. This study is extended to the IVFFS environment due to the ability to model the uncertainty and hesitation in decision-makers’ evaluations at a more granular level.

Preliminaries on IVFFSs

Some preliminary information about IVFFs and related notations is provided in this subsection. FFSs are an extension of orthopair fuzzy clusters that can more easily process ambiguous information in the decision-making process (Garg, Shahzadi & Akram, 2020). The concept of FFS was developed by Senapati & Yager (2020, 2019) as a generalisation of intuitive fuzzy sets (IFSs) and Pythagorean fuzzy sets (PFSs). In IFSs, the sum of the membership and non-membership degrees should be at most one, and in PFSs, the sum of their squares should be at most one (Senapati & Yager, 2020). To offer more flexibility to decision-makers, FFS has been developed such that the sum of the cubes of the membership and non-membership degrees is equal to 1 (see Fig. 2). In decision-making processes, decision-makers are confronted with hesitation and partial information that other fuzzy numbers cannot fully express (Simić et al., 2022). FFSs evaluated the uncertainty more comprehensively, taking into account the hesitations and partial information that decision-makers cannot fully articulate. This feature makes FFSs especially effective in situations such as risk assessment, in which there are opinions of decision-makers and uncertain data. FFSs increase credibility in the decision-making process by reflecting expert opinions more accurately. FFSs are more convenient and better suited than other indeterminate cluster extensions to address uncertainty by assigning membership and non-membership rank parameters to a larger domain. This set can be described as an innovative approach to represent unreliable, imprecise, and uncertain information in a fuzzy environment (Garg, Shahzadi & Akram, 2020; Simić et al., 2022). As a generalisation of FFSs, IVFFs are used in this study because they allow flexible comparison of range-valued numerical data. For example, suppose a decision maker defines the membership value as [0.65, 0.7] non-membership value as [0.60, 0.75] of an alternative. Here, the sum of the squares of the upper bound of the membership value (0.7)² and the upper bound of the non-membership value (0.75)² is greater than 1. In this case, it is neither IVIFS nor IVPFS. However, (0.7)³ + (0.75)³ = 0.7648 ≤ 1. In this case, it can be considered as IVFFS.

Linguistic terms are used by decision-makers in the assessment of criteria and hazards. These linguistic terms utilise fuzzy numbers to express the uncertainty of abstract information in a more flexible, understandable, and realistic way in the evaluation processes conducted by decision makers. For this reason, the linguistic terms and their corresponding IVFF number equivalents are given in Table 2 (Jeevaraj, 2021).

Operations for IVFFs are given below (Jeevaraj, 2021):

Definition 1. If we define the universe of discourse as X, the fermatean fuzzy number in this set is presented as in Eq. (1) (Jeevaraj, 2021).

(1) $$\tilde{\mathrm{F}}=\left\{x,{\mu }_{\tilde{F}}\left(x\right),v_{\tilde{F}}\left(x\right);x\in X\right\}$$where ${\mu }_{\tilde{F}}\left(x\right)\mapsto \left[0,1\right]$ and $v_{\tilde{F}}\left(x\right)\mapsto \left[0,1\right]$ shows the closed intervals of the membership and non-membership degrees of the element $x\in X$ to the set $\tilde{\mathrm{F}}$, respectively. Also, the lower and upper limits are denoted by ${\mu }_{\tilde{F}}^L\left(x\right),{\mu }_{\tilde{F}}^U\left(x\right),v_{\tilde{F}}^L\left(x\right),v_{\tilde{F}}^U\left(x\right),$ respectively. Therefore, $\tilde{\mathrm{F}}$ can also be expressed as in Eqs. (2) and (3).

(2) $${\mu }_{\tilde{F}}\left(x\right)=\left[{\mu }_{\tilde{F}}^L\left(x\right),{\mu }_{\tilde{F}}^U\left(x\right)\right]\mapsto \left[0,1\right]$$ (3) $$v_{\tilde{F}}\left(x\right)=\left[v_{\tilde{F}}^L\left(x\right),v_{\tilde{F}}^U\left(x\right)\right]\mapsto \left[0,1\right]$$ (4) $$0\le {\mu }_{\tilde{F}}^U{\left(x\right)}^3+v_{\tilde{F}}^U{\left(x\right)}^3\le 1.$$For each $x\in X,{\pi }_{\tilde{F}}\left(x\right)=\left[{\pi }_{\tilde{\mathcal{F}}}^L\left(x\right),{\pi }_{\tilde{\mathcal{F}}}^U\left(x\right)\right]$ is called as the hesitancy degree in IVFFSs. Lower and upper limits of the hesitancy degree:

(5) $${\pi }_{\tilde{F}}^L\left(x\right)=\sqrt[3]{1-({\mu }_{\tilde{F}}^U\left(x\right){)}^3+(v_{\tilde{F}}^U\left(x\right){)}^3}$$ (6) $${\pi }_{\tilde{F}}^U\left(x\right)=\sqrt[3]{1-({\mu }_{\tilde{F}}^L\left(x\right){)}^3+(v_{\tilde{F}}^L\left(x\right){)}^3}.$$

Definition 2. $\tilde{\mathrm{F}}=\left(\left[{\mu }_{\tilde{F}}^L\left(x\right),{\mu }_{\tilde{F}}^U\left(x\right)\right],\left[v_{\tilde{F}}^L\left(x\right),v_{\tilde{F}}^U\left(x\right)\right]\right)$, $$\widetilde{{\mathrm{F}}_1}=\left(\left[{\mu }_{\tilde{F}1}^L\left(x\right),{\mu }_{\tilde{F}1}^U\left(x\right)\right],\left[v_{\tilde{F}1}^L\left(x\right),v_{\tilde{F}1}^U\left(x\right)\right]\right)$$ and ${\tilde{\mathrm{F}}}_2=\left(\left[{\mu }_{\widetilde{F2}}^L\left(x\right),{\mu }_{\widetilde{F2}}^U\left(x\right)\right],\left[v_{\tilde{F}2}^L\left(x\right),v_{\widetilde{F2}}^U\left(x\right)\right]\right)$ be three IVFFSs $\lambda >0,$ their operations are as follows (Jeevaraj, 2021):

(7) $${\tilde{\mathcal{F}}}_1\oplus {\tilde{\mathcal{F}}}_2=\left(\left[\begin{array}{c}\sqrt[3]{{\left({\mu }_{{\tilde{\mathcal{F}}}_1}^L\right)}^3+{\left({\mu }_{{\tilde{\mathcal{F}}}_2}^L\right)}^3-{\left({\mu }_{{\tilde{\mathcal{F}}}_1}^L\right)}^3{\left({\mu }_{{\tilde{\mathcal{F}}}_2}^L\right)}^3,}\\ \sqrt[3]{{\left({\mu }_{{\tilde{\mathcal{F}}}_1}^U\right)}^3+{\left({\mu }_{{\tilde{\mathcal{F}}}_2}^U\right)}^3-{\left({\mu }_{{\tilde{\mathcal{F}}}_1}^U\right)}^3{\left({\mu }_{{\tilde{\mathcal{F}}}_2}^U\right)}^3}\end{array}\right],\left[v_{{\tilde{\mathcal{F}}}_1}^L{v}_{{\tilde{\mathcal{F}}}_2}^L,v_{{\tilde{\mathcal{F}}}_1}^U{v}_{{\tilde{\mathcal{F}}}_2}^U\right]\right)$$ (8) $${\tilde{\mathcal{F}}}_1\otimes {\tilde{\mathcal{F}}}_2=\left(\left[{\mu }_{{\tilde{\mathcal{F}}}_1}^L{\mu }_{{\tilde{\mathcal{F}}}_2}^L,{\mu }_{{\tilde{\mathcal{F}}}_1}^U{\mu }_{{\tilde{\mathcal{F}}}_2}^U\right],\left[\begin{array}{c}\sqrt[3]{{\left(v_{{\tilde{\mathcal{F}}}_1}^L\right)}^3+{\left(v_{{\tilde{\mathcal{F}}}_2}^L\right)}^3-{\left(v_{{\tilde{\mathcal{F}}}_1}^L\right)}^3{\left(v_{{\tilde{\mathcal{F}}}_2}^L\right)}^3}\\ \sqrt[3]{{\left(v_{{\tilde{\mathcal{F}}}_1}^U\right)}^3+{\left(v_{{\tilde{\mathcal{F}}}_2}^U\right)}^3-{\left(v_{{\tilde{\mathcal{F}}}_1}^U\right)}^3{\left(v_{{\tilde{\mathcal{F}}}_2}^U\right)}^3}\end{array}\right]\right)$$ (9) $$\lambda \tilde{\mathcal{F}}=\left(\left[\sqrt[3]{1-{\left(1-{\left({\mu }_{\tilde{\mathcal{F}}}^L\right)}^3\right)}^{\lambda }},\sqrt[3]{1-{\left(1-{\left({\mu }_{\tilde{\mathcal{F}}}^U\right)}^3\right)}^{\lambda }}\right],\left[{\left(v_{\tilde{\mathcal{F}}}^L\right)}^{\lambda },{\left(v_{\tilde{\mathcal{F}}}^U\right)}^{\lambda }\right]\right)$$ (10) $${\tilde{\mathcal{F}}}^{\lambda }=\left(\left[{\left({\mu }_{\tilde{\mathcal{F}}}^L\right)}^{\lambda },{\left({\mu }_{\tilde{\mathcal{F}}}^U\right)}^{\lambda }\right],\left[\sqrt[3]{1-{\left(1-{\left(v_{\tilde{\mathcal{F}}}^L\right)}^3\right)}^{\lambda }},\sqrt[3]{1-{\left(1-{\left(v_{\tilde{\mathcal{F}}}^U\right)}^3\right)}^{\lambda }}\right]\right).$$Definition 3. Let ${\tilde{\mathrm{F}}}_i=\left(\left[{\mu }_{{\tilde{\mathcal{F}}}_i}^L\left(x\right),{\mu }_{{\tilde{\mathcal{F}}}_i}^U\left(x\right)\right],\left[v_{{\tilde{\mathcal{F}}}_i}^L\left(x\right),v_{{\tilde{\mathcal{F}}}_i}^U\left(x\right)\right]\right)$ ( $i$ = 1, 2, …,n) be a set of IVFFSs and $w={\left(w_1,w_2,\dots ,w_n\right)}^T$ be a weight vector of ${\mathcal{F}}_i$ with $\sum _{i=1}^n{w}_i=1$ (Jeevaraj, 2021),

Interval-Valued Fermatean Fuzzy Weighted Average (IVFFWA) operator is a mapping IVFFWA:

${\tilde{\mathcal{F}}}^n\to \tilde{\mathcal{F}},$ where $\mathrm{I}\mathrm{V}\mathrm{F}\mathrm{F}\mathrm{W}\mathrm{A}({\tilde{\mathcal{F}}}_1,{\tilde{\mathcal{F}}}_2,\dots ,{\tilde{\mathcal{F}}}_n)$

(11) $$\begin{array}{rl}=& ([\sqrt[3]{(1-\prod _{i=1}^n(1-({\mu }_{{\tilde{\mathcal{F}}}_i}^L{)}^3{)}^{w_i})},\sqrt[3]{(1-\prod _{i=1}^n(1-({\mu }_{{\tilde{\mathcal{F}}}_i}^U{)}^3{)}^{w_i})}]\\ & \times [\prod _{i=1}^n{\left(v_{{\tilde{\mathcal{F}}}_i}^L\right)}^{w_i},\prod _{i=1}^n(v_{{\tilde{\mathcal{F}}}_i}^U{)}^{w_i}]).\end{array}$$Interval-valued fermatean fuzzy weighted geometric (IVFFWG) operator is a mapping $\mathrm{I}\mathrm{V}\mathrm{F}\mathrm{F}\mathrm{W}\mathrm{A}:{\tilde{\mathcal{F}}}^n\to \tilde{\mathcal{F}},$ where $IVFFWA({\tilde{\mathcal{F}}}_1,{\tilde{\mathcal{F}}}_2,\dots ,{\tilde{\mathcal{F}}}_n)$

(12) $$\begin{array}{rl}=& ([{\prod }_{i=1}^n{\left({\mu }_i^L\right)}^{w_i},\prod _{i=1}^n{\left({\mu }_i^U\right)}^{w_i}]\\ & \times [3\left(1-\prod _{i=1}^n{\left(1-{\left(v_{{\tilde{\mathcal{F}}}_i}^L\right)}^3\right)}^{w_i}\right),\sqrt[3]{\left(1-\prod _{i=1}^n{\left(1-{\left(v_{{\tilde{\mathcal{F}}}_i}^U\right)}^3\right)}^{w_i}\right)]}]).\end{array}$$

IVFFAHP-TOPSIS method

This section presents the main steps and equations of the IVFFAHP-TOPSIS integrated method. The IVFFAHP method was developed by Alkan & Kahraman (2023), and this study follows the steps suggested for IVFFAHP. The equations of the IVFFAHP method are given from Step 1 to Step 10. Then, in this study, the classical TOPSIS method is integrated into the IVFFAHP method. The equations of the TOPSIS method are given from Step 11 to Step 15.

Step 1. Create the hierarchical framework by determining the criteria and alternatives. The set $A_i=\left(A_1,A_2,\dots ,A_n\right)$, having i $=\left(1,2,\dots ,n\right)$ alternatives, is evaluated by $m$ criteria of set $C_j=\left(C_1,C_2,\dots ,C_m\right)$, with j $=\left(1,2,\dots ,m\right)$. Let $w_j=\left(w_1,w_2,\dots ,w_m\right)$ be the vector set represented the criteria weights, where $w_j>0$ and $\sum _{j=1}^m{w}_j=1$.

Step 2. Decision makers are evaluated and the weight for each decision maker is determined. Table 3 is used to evaluate decision makers. At the same time, the weight of $t$ experts are denoted by ${\psi }^t$, where $k$ is the number of experts. The sum of all expert weights is $\sum _{t=1}^k{\psi }^t=1$. Table 2 presents linguistic terms and their corresponding IVFFSs. Let d be the number of decision makers, and $D{M}_t=\left({\mu }_{\mathrm{t}},v_{\mathrm{t}}\right)$ be the corresponding FFSs to determine weight of decision maker t based on evaluation. Each decision maker is evaluated individually, and the weight of decision maker is determined according to Eq. (13).

(13) $${\psi }_t={\displaystyle \frac{{\mathrm{S}}^{+}\left(D{M}_t\right)}{\sum _{t=1}^d{S}^{+}\left(D{M}_t\right)}={\displaystyle \frac{1+{\mu }_{\mathrm{t}}^3-v_{\mathrm{t}}^3}{\sum _{t=1}^d\left(1+{\mu }_{\mathrm{t}}^3-v_{\mathrm{t}}^3\right)}.}}$$

Step 3. Created the pairwise comparison matrix $Z$ $={\left(z_{ij}\right)}_{m\times m}$ based on the linguistic variables of decision makers given in Table 2.

(14) $$Z=\left[\begin{array}{cccc}1& z_{12}& \cdots & z_{1m}\\ z_{21}& 1& \cdots & z_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ z_{m1}& z_{m2}& \cdots & 1\end{array}\right],\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{z}_{ij}=⟨\left[{\mu }_{ij}^L,{\mu }_{ij}^U\right],\left[v_{ij}^L,v_{ij}^U\right]⟩.$$

Step 4. Calculated for the consistency ratio (CR) of each pairwise comparison matrix ( $Z$). The CR of the matrix is determined according to Satty’s classical consistency procedure. The random index (RI) is given in File S11. λ_max represents the highest eigenvalue, while n represents the number of criteria.

(15) $$CR={\displaystyle \frac{\mathrm{C}\mathrm{I}\left(\mathrm{C}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{y}\mathrm{I}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}\right)}{\mathrm{R}\mathrm{I}}}$$

(16) $$CI={\displaystyle \frac{\lambda \mathrm{\_}max-n}{\mathrm{n}-1}.}$$This step ensures the logical coherence and internal consistency of the decision-maker judgments. By applying Saaty’s CR, the reliability of the pairwise comparison matrices are validated. This step is essential to minimize subjective inconsistency in the decision-making process.

Step 5. Aggregate the opinions of decision makers. The pairwise comparison matrix constituted for each decision maker is aggregated by using IVFFWA aggregation operator given in Eq. (17). Let $A_{ij}^t=$ $\left(\left[{\mu }_{ij}^{Lt},{\mu }_{ij}^{Ut}\right],\left[v_{ij}^{Lt},v_{ij}^{Ut}\right]\right)$ be the pairwise comparison of criteria $i$ and $j$ by decision maker t. $\mathrm{I}\mathrm{V}\mathrm{F}\mathrm{F}\mathrm{W}\mathrm{A}({\tilde{\mathcal{F}}}_1,{\tilde{\mathcal{F}}}_2,\dots ,{\tilde{\mathcal{F}}}_n)$

(17) $$\begin{array}{rl}=& ([\sqrt[3]{\left(1-\prod _{i=1}^n{\left(1-{\left({{\mu }_{{\tilde{\mathcal{F}}}_i}}^L\right)}^3\right)}^{w_i}\right)},\sqrt[3]{\left(1-\prod _{i=1}^n{\left(1-{\left({\mu }_{{{\tilde{\mathcal{F}}}_i}^U}\right)}^3\right)}^{w_i}\right)}]\\ & \times \left[\prod _{i=1}^n{\left(v_{{\tilde{\mathcal{F}}}_i}^L\right)}^{w_i},\prod _{i=1}^n{\left(v_{{\tilde{\mathcal{F}}}_i}^U\right)}^{w_i}\right]).\end{array}$$

Step 6. Find the differences matrix $\mathrm{D}={\left({\mathrm{d}}_{\mathrm{i}\mathrm{j}}\right)}_{\mathrm{m}\times \mathrm{m}}$ between lower and upper points of the membership and non-membership functions using Eqs. (18) and (19).

(18) $$d_{ij}^L={\left({\mu }_{ij}^L\right)}^3-{\left(v_{ij}^U\right)}^3$$

(19) $$d_{ij}^U={\left({\mu }_{ij}^U\right)}^3-{\left(v_{ij}^L\right)}^3.$$In this step, the differences between membership and non-membership degrees are calculated and the superiority and uncertainty levels of the criteria compared to each other are quantitatively revealed. The differences obtained form the basis for a more accurate calculation of the criterion importance in the next steps.

Step 7. Find the interval multiplicative matrix $S={\left(s_{ij}\right)}_{m\times m}$ using Eqs. (20) and (21).

(20) $$s_{ij}^L=\sqrt[3]{{1000}^{d_{ij}^L}}$$

(21) $$s_{ij}^U=\sqrt[3]{{1000}^{d_{ij}^U}}.$$In this step, exponential transformation is applied to increase or decrease the effect of differences. Thanks to this transformation, small differences are suppressed, while the effect of large differences is highlighted. This makes it easier to identify critical criteria in the decision-making process.

Step 8. Obtain the indeterminacy value $T={\left(t_{ij}\right)}_{m\times m}$ of the $z_{ij}$ is Eq. (22).

(22) $$t_{ij}=1-\left({\mu }_{i{j}_U}^3-{\mu }_{i{j}_L}^3\right)-\left(v_{i{j}_U}^3-v_{i{j}_L}^3\right).$$This step reflects the degree of indeterminacy between assessments. Uncertainty in an FF environment refers to the level of indecision in decision-makers’ assessments, and the inclusion of these indeterminacies in weight calculations increases the reliability of evaluations.

Step 9. Multiply the indeterminacy degrees with $S={\left(s_{ij}\right)}_{m\times m}$ matrix to find the matrix of unnormalized weights $R={\left(r_{ij}\right)}_{m\times m}$ using Eq. (23).

(23) $$r_{ij}=\left({\displaystyle \frac{s_{ij}^L+s_{ij}^U}{2}}\right)t_{ij}.$$In this step, criterion weights are determined by considering both the relative differences between the criteria and the indeterminacy levels of the evaluations. Thus, both the importance and the reliability element were included in the evaluation at the same time.

Step 10. Obtain the normalised priority weights w_i by using Eq. (24).

(24) $$w_i={\displaystyle \frac{{\sum }_{j=1}^m{r}_{ij}}{{\sum }_{i=1}^m{\sum }_{j=1}^m{r}_{ij}.}}$$This step allows the sum of the weights of all criteria to be normalized to 1. Thus, decision-makers can clearly and comparably see the relative importance of each criterion in the overall decision process.

Step 11. A weighted decision matrix $V={\left(v_{ij}\right)}_{n\times m}$ is created using the criteria and alternatives. Where $w_j$ is the criterion weight. $a_{ij}$ is the evaluation of the i. alternative according to the j. criterion.

(25) $$v_{ij}=w_j{a}_{ij}.$$

Step 12. Positive and negative ideal solution values are calculated. The maximum value of each column is called the positive ideal solution $\left(v_j^{+}\right)$, and the minimum value is called the negative ideal solution $\left(v_j^{-}\right)$.

$v_j^{+}$ = ( $v_1,v_2,\dots ,v_j)$ maximum values, $v_j^{-}$ = ( $v_1,v_2,\dots ,v_j)$ minimum values.

Step 13. Positive and negative distance values are determined.

(26) $$S^{+}=\sqrt{\sum _{j=1}^m{\left(v_{ij}-v_j^{+}\right)}^2}$$

(27) $$S^{-}=\sqrt{\sum _{j=1}^m{\left(v_{ij}-v_j^{-}\right)}^2}.$$

Step 14. The closeness coefficient is calculated according to the ideal solution. Then the closeness coefficient values are normalized.

(28) $$C_i={\displaystyle \frac{S_i^{-}}{S_i^{-}+S_i^{+}}}$$

(29) $$C_i^{\ast }={\displaystyle \frac{C_i}{{\sum }_{i=1}^n{C}_i}.}$$

Step 15. Ranked from largest to smallest according to the normalised closeness coefficient.

Hierarchy of evaluation criteria, alternatives and decision makers

Risk assessment is carried out by systematically analysing the identified risk parameters and hazards. Risk parameters and hazards in the hierarchical structure in Fig. 3 are determined in line with the literature review, process analysis and expert opinions. Risk parameters are determined to assess the health and safety risks of workers working on plastic injection moulding machines. Probability (P1), severity (P2), frequency (P3), detectability (P4), human error (P5), machine error (P6) and existing safety measures (P7) are used as risk parameters. A total of nine potential hazards that may arise from the use of plastic injection moulding machines have been identified. These hazards and their descriptions are presented in Table 4.

Decision-makers are determined according to their experience in the relevant field (at least 5 years), their level of education (at least high school) and their titles (employer, occupational safety specialist, workplace physician, production manager, employee representative, experienced employee, etc.). Details of the decision-makers are given in Table 5.

Since the decision-making methods used in the risk assessment process take into account subjective judgements, it is aimed to create a decision-making group that can fully comprehend the meaning and importance of the assessment to be made. Accordingly, different SDMGs are determined for each evaluation area. SDMGs are determined in line with the common opinions of the study team, academics and field experts. SDMGs are determined separately according to each assessment area. Details of the SDMGs and evaluation areas are presented in Table 6.

Determining the weights of the risk parameters using IVFF-AHP

In this phase, the weight of the risk parameters is calculated according to the procedures from Step 1 to Step 10 of the IVFF-AHP method. Decision-makers use the linguistic terms and scales presented in Table 2 when assessing risk parameters. An example of calculating the weight of risk parameters is given below.

Step 1: After determining the risk parameters, a hierarchical framework is created as shown in Fig. 3. To evaluate the risk parameters, the occupational safety specialist (DM1), the production manager (DM2), academician (DM3) in the 1st sub-decision maker group (SDMG1) are determined as decision makers.

Step 2: Decision-makers are evaluated according to the linguistic terms given in Table 3. The decision makers determined for the risk parameters are evaluated by the working team, academics, and practitioners as occupational safety specialist (VS), production manager (MS) and academician (VS) according to the linguistic terms in Table 3. Linguistic terms are converted into FFs, and the weight of decision-makers is determined by Eq. (13). For example, the weight of the Occupational safety specialist (DM1) decision-maker is determined as follows:

$D{M}_1=\left({\mu }_1,v_1\right)=\left(0.75,0.30\right)$;

$$\begin{array}{rl}{\psi }_1=& {\displaystyle \frac{1+{\mu }_1^3-v_1^3}{\sum _{t=1}^3\left(1+{\mu }_{\mathrm{t}}^3-v_{\mathrm{t}}^3\right)}}\\ =& {\displaystyle \frac{1+{\left(0.75\right)}^3-{\left(0.30\right)}^3}{\left(1+{\left(0.75\right)}^3-{\left(0.30\right)}^3\right)+\left(1+{\left(0.55\right)}^3-{\left(0.50\right)}^3\right)+\left(1+{\left(0.75\right)}^3-{\left(0.30\right)}^3\right)}}\\ =& \mathrm{0.3641.}\end{array}$$

Step 3: A pairwise comparison matrix of risk parameters is created. Decision-makers assess risk parameters using the linguistic terms in Table 2. The pairwise comparison matrix of risk parameters is given in File S12. These linguistic terms are converted into IVFFSs given in Table 2.

Step 4: The CR of each pairwise comparison matrix is calculated. The CR of the matrix is determined according to Satty’s classical consistency procedure. For example, the CR value of the pairwise comparison matrix of the risk parameters evaluated by the occupational safety specialist (DM1) is calculated as follows. It is obtained as λ_max = 7.042 and n = 7 belonging to the pairwise comparison matrix. According to Eq. (15), CI = (λ_max-n)/(n-1) = (7.042 – 7)/( 7−1) = 0.007. To obtain the CR value, the RI value in the equation is determined as (RI = 1.32) for n = 7. It is then obtained as CR = CI/RI = 0.007/1.32 = 0.005. The CR ratios for each matrix are given in the left column of the tables.

Step 5: The decision-makers’ evaluations are combined using the IVFFWA operator given in Eq. (17). For example, the evaluations Z12 (SLI) of DM1, Z12 (LI) of DM2, Z12 (LI) of DM3 are combined as follows using the IVFFWA operator. These linguistic terms are converted into IVFFSs in Table 2. The linguistic terms in the pairwise comparison matrix of each decision maker are converted into IVFF numbers in the form of $A_{12}^1=\left(\left[0.35,0.4\right],\left[0.6,0.65\right]\right),A_{12}^2=\left(\left[0.2,0.3\right],\left[0.7,0.8\right]\right),A_{12}^3=\left(\left[0.35,0.4\right],\left[0.6,0.65\right]\right)$. Then, the IVFFWA operator operations in Eq. (17) are applied. $\mathrm{I}\mathrm{V}\mathrm{F}\mathrm{F}\mathrm{W}\mathrm{A}({\tilde{\mathcal{F}}}_1,{\tilde{\mathcal{F}}}_2,\dots ,{\tilde{\mathcal{F}}}_n)$

$$\begin{array}{rl}=& (\left[\begin{array}{c}\hfill \sqrt[3]{\left(1-{\left(1-{\left(0.35\right)}^3\right)}^{0.3641}\ast {\left(1-{\left(0.35\right)}^3\right)}^{0.2718}\ast {\left(1-{\left(0.35\right)}^3\right)}^{0.3641}\right)},\\ \hfill \sqrt[3]{\left(1-{\left(1-{\left(0.4\right)}^3\right)}^{0.3641}\ast {\left(1-{\left(0.3\right)}^3\right)}^{0.2718}\ast {\left(1-{\left(0.4\right)}^3\right)}^{0.3641}\right)}\end{array}\right]\\ & \times \left[{\left(0.6\right)}^{0.3641}\ast {\left(0.7\right)}^{0.2718}\ast {\left(0.6\right)}^{0.3641},{\left(0.65\right)}^{0.3641}\ast {\left(0.65\right)}^{0.2718}\ast {\left(0.65\right)}^{0.3641}\right]).\end{array}$$

$\mathrm{I}\mathrm{V}\mathrm{F}\mathrm{F}\mathrm{W}\mathrm{A}({\tilde{\mathcal{F}}}_1,{\tilde{\mathcal{F}}}_2,{\tilde{\mathcal{F}}}_3)$ = $\left(\left[0.35,0.37\right],\left[0.62,0.65\right]\right)$ is obtained.

Step 6: The separation between the upper and lower points in the calculation matrix is calculated according to Eqs. (18) and (19). The separation value between the upper and lower points $\mathrm{D}={\left({\mathrm{d}}_{12}\right)}_{7\times 7}$ evaluated in step 5 is calculated as follows.

$d_{12}^L={\left(0.37\right)}^3-{\left(0.62\right)}^3$ = −0.1877

$d_{12}^U={\left(0.35\right)}^3-{\left(0.65\right)}^3$ = −0.2318

Step 7: The interval generative matrix $S={\left(s_{12}\right)}_{7\times 7}$ is calculated according to Eqs. (20) and (21).

$s_{12}^L=\sqrt[3]{{1000}^{d_{12}^L}}$ = $\sqrt[3]{{1000}^{-0.1877}}$ = 0.6491

$s_{12}^U=\sqrt[3]{{1000}^{d_{12}^U}}$ = $\sqrt[3]{{1000}^{-0.2318}}$ = 0.5864.

Step 8. The uncertainty value $T={\left(t_{12}\right)}_{7\times 7}$ of the $z_{12}$ is obtained.

$t_{12}=1-\left({0.35}^3-{0.37}^3\right)-\left({0.62}^3-{0.65}^3\right)$ = 1.0441.

Step 9. The uncertainty levels are multiplied with $S={\left(s_{12}\right)}_{7\times 7}$ matrix to find the matrix of unnormalized weights $R={\left(r_{12}\right)}_{7\times 7}$ by Eq. (23).

$r_{12}=\left({\displaystyle \frac{s_{12}^L+s_{12}^U}{2}}\right)t_{ij}$ = $\left({\displaystyle \frac{0.6491+0.5864}{2}}\right)\ast 1.0441$ = 0.6450.

Step 10. The computation of the normalised priority weights w_i is as below.

$w_1={\displaystyle \frac{{\sum }_{j=1}^m{r}_{ij}}{{\sum }_{i=1}^m{\sum }_{j=1}^m{r}_{ij}}}$ = 0.1445.

The weights of the risk parameters are determined according to the IVFF-AHP calculation steps in the example above. The weights of the decision makers, linguistic assessments of the risk parameters, and consistency ratios are presented in File S12. The consistency ratio of the pairwise comparison matrices is calculated using Excel during the evaluations. The MATLAB code of the calculations of the IVFF-AHP method is given in the File S2. All risk parameters according to highest weight are the existing safety measures (0.1987), detectability (0.1762), machine error (0.1573), human error (0.1558), probability (0.1445), frequency (0.1195) and severity (0.0481), respectively.

Determination of hazard scores using IVFF-AHP

In this phase, the hazard scores are calculated according to the procedures from Step 1 to Step 10 of the IVFF-AHP method. Decision-makers use the linguistic terms and scales presented in Table 2 when assessing hazards. Pairwise comparison matrices of hazards are created according to each risk parameter. For example, the calculation of the scores of the hazards according to the detectability parameter is given below.

Step 1: After determining the risk parameters and hazards, a hierarchical framework is created as shown in Fig. 3. To evaluate the hazards according to the detectability parameter, the occupational safety specialist (DM1), the production manager (DM2), the academician (DM3) and the employer (DM7), who are in the 5th sub-decision maker group (SDMG5), are determined as decision makers.

Step 2: Decision makers are evaluated according to the linguistic terms given in Table 3. The decision-makers determined for the detectability parameter are evaluated by the working team, academics, and practitioners as occupational safety specialist (AS), production manager (AS), academician (AS) and employer (MS) according to the linguistic terms in Table 3. Linguistic terms are converted into FFs, and the weight of decision-makers is determined by Eq. (13). For example, the weight of the occupational safety specialist decision-maker (DM1) is determined as follows:

$D{M}_1=\left({\mu }_{1,}{v}_1\right)=\left(0.95,0.10\right)$;

$$\begin{array}{rl}{\psi }_1=& {\displaystyle \frac{1+{\mu }_1^3-v_1^3}{\sum _{t=1}^4\left(1+{\mu }_{\mathrm{t}}^3-v_{\mathrm{t}}^3\right)}}\\ =& {\displaystyle \frac{1+{\left(0.95\right)}^3-{\left(0.10\right)}^3}{\left(1+{\left(0.95\right)}^3-{\left(0.10\right)}^3\right)+\left(1+{\left(0.95\right)}^3-{\left(0.10\right)}^3\right)+\left(1+{\left(0.95\right)}^3-{\left(0.10\right)}^3\right)+\left(1+{\left(0.55\right)}^3-{\left(0.50\right)}^3\right)}}\\ =& 0.2808.\end{array}$$

Step 3: A pairwise comparison matrix of hazards is created. Decision-makers assess hazards based on the detectability parameter using the linguistic terms in Table 2. The pairwise comparison matrix of hazards according to the detectability parameter is given in File S13. These linguistic terms are converted into IVFFSs given in Table 2.

Step 4: The CR of each pairwise comparison matrix is calculated. The CR of the matrix is determined according to Satty’s classical consistency procedure. For example, the CR value of the pairwise comparison matrix of hazards assessed by the occupational safety specialist (DM1) is calculated as follows. It is obtained as λ_max = 10.114 and n = 9 belonging to the pairwise comparison matrix. According to Eq. (15), CI = (λ_max-n)/(n-1) = (10.114 – 9)/(9 − 1) = 0.139. To obtain the CR value, the RI value in the equation is determined as (RI = 1.45) for n = 9. It is then obtained as CR = CI/RI = 0.139/1.45 = 0.096. The CR ratios for each matrix are given in the left column of the tables.

Step 5: The decision-makers’ evaluations are combined using the IVFFWA operator given in Eq. (17). For example, the evaluations Z12 (SLI) of DM1, Z12 (SLI) of DM2, Z12 (LI) of DM3 and Z12 (LI) of DM7 are combined as follows using the IVFFWA operator. These linguistic terms are converted into IVFFSs in Table 2. The linguistic terms in the pairwise comparison matrix of each decision maker are converted into IVFF numbers in the form of $A_{12}^1=\left(\left[0.35,0.4\right],\left[0.6,0.65\right]\right)$, $A_{12}^2=\left(\left[0.35,0.4\right],\left[0.6,0.65\right]\right)$, $A_{12}^3=\left(\left[0.2,0.3\right],\left[0.7,0.8\right]\right)$ ve $A_{12}^7=\left(\left[0.2,0.3\right],\left[0.7,0.8\right]\right)$. Then, the IVFFWA operator operations in Eq. (17) are applied. $\mathrm{I}\mathrm{V}\mathrm{F}\mathrm{F}\mathrm{W}\mathrm{A}({\tilde{\mathcal{F}}}_1,{\tilde{\mathcal{F}}}_2,\dots ,{\tilde{\mathcal{F}}}_n)$

$$\begin{array}{rl}=& (\left[\begin{array}{c}\hfill \sqrt[3]{\left(1-{\left(1-{\left(0.35\right)}^3\right)}^{0.2808}\ast {\left(1-{\left(0.35\right)}^3\right)}^{0.2808}\ast {\left(1-{\left(0.2\right)}^3\right)}^{0.2808}\right)\ast {\left(1-{\left(0.2\right)}^3\right)}^{0.1575})},\\ \hfill \sqrt[3]{\left(1-{\left(1-{\left(0.4\right)}^3\right)}^{0.2808}\ast {\left(1-{\left(0.4\right)}^3\right)}^{0.2808}\ast {\left(1-{\left(0.3\right)}^3\right)}^{0.2808}\right)\ast {\left(1-{\left(0.3\right)}^3\right)}^{0.1575}}\end{array}\right]\\ & \times \left[{\left(0.6\right)}^{0.2808}\ast {\left(0.6\right)}^{0.2808}\ast {\left(0.7\right)}^{0.2808}\ast {\left(0.7\right)}^{0.1575},{\left(0.65\right)}^{0.2808}\ast {\left(0.65\right)}^{0.2808}\ast {\left(0.8\right)}^{0.2808}\ast {\left(0.8\right)}^{0.1575}\right])\end{array}$$

$\mathrm{I}\mathrm{V}\mathrm{F}\mathrm{F}\mathrm{W}\mathrm{A}({\tilde{\mathcal{F}}}_1,{\tilde{\mathcal{F}}}_2,{\tilde{\mathcal{F}}}_3)$ = $\left(\left[0.29,0.37\right],\left[0.72,0.70\right]\right).$

Step 6: The difference between the upper and lower points in the calculation matrix is calculated according to Eqs. (18) and (19). The difference value between the upper and lower points $\mathrm{D}={\left({\mathrm{d}}_{12}\right)}_{9\times 9}$ evaluated in step 5 is calculated as follows.

$d_{12}^L={\left(0.37\right)}^3-{\left(0.72\right)}^3$ = −0.3226

$d_{12}^U={\left(0.29\right)}^3-{\left(0.70\right)}^3$ = -0.3186.

Step 7: The interval generative matrix $S={\left(s_{12}\right)}_{9\times 9}$ is calculated according to Eqs. (20) and (21).

$s_{12}^L=\sqrt[3]{{1000}^{d_{12}^L}}$ = $\sqrt[3]{{1000}^{-0.3226}}$ = 0.7777

$s_{12}^U=\sqrt[3]{{1000}^{d_{12}^U}}$ = $\sqrt[3]{{1000}^{-0.3186}}$ = 0.7790

Step 8. The uncertainty value $T={\left(t_{12}\right)}_{9\times 9}$ of the $z_{12}$ is obtained.

$t_{12}=1-\left({\left(0.29\right)}^3-{\left(0.37\right)}^3\right)-\left({\left(0.72\right)}^3-{\left(0.70\right)}^3\right)$ = 0.9960.

Step 9. The uncertainty levels are multiplied with $S={\left(s_{12}\right)}_{9\times 9}$ matrix to find the matrix of unnoted weights $R={\left(r_{12}\right)}_{9\times 9}$ by Eq. (23).

$r_{ij}=\left({\displaystyle \frac{s_{12}^L+s_{12}^U}{2}}\right)t_{12}$ = $\left({\displaystyle \frac{0.7777+0.7790}{2}}\right)\ast 0.9960$ = 0.7764.

Step 10. The computation of the normalised priority weights w_i is as below.

$w_1={\displaystyle \frac{{\sum }_{j=1}^9{r}_{ij}}{{\sum }_{i=1}^9{\sum }_{j=1}^9{r}_{ij}}}$ = 0.0798.

The scores of the hazards are determined according to the IVFFAHP calculation steps in the example given above. The weights of the decision makers, linguistic assessments of hazards assessed according to the parameters of probability, severity, frequency, detectability, human error, machine error, and existing safety precautions, and consistency ratios are presented in the from Files S13 to S19. The consistency ratio of the pairwise comparison matrices is calculated using Excel during the evaluations. The MATLAB codes of the calculations of the IVFF-AHP method are given in the from Files S3 to S9. The hazard scores obtained according to each risk parameter are given in Table 7.

Determination of risk coefficient values using TOPSIS

In this phase, after the hazard scores are obtained with the IVFFAHP method, the procedures from Step 11 to Step 15 of the TOPSIS method are applied to obtain the risk coefficient values. Hazard scores and risk parameter weights obtained by the IVFFAHP method are directly transferred to the weighted decision matrix. In this study, detectability and existing safety measures are defined as cost criteria, while other risk parameters are defined as benefit criteria. The risk coefficient values of the hazards are obtained by the TOPSIS method. According to the risk coefficient values, hazards are ranked starting from the highest risk hazard. An example of the calculation of the ranking of hazards is given below.

Step 11. A weighted decision matrix $V={\left(v_{11}\right)}_{9\times 7}$ is created using the risk parameters and hazard scores. Where $w_1$ is the risk parameters weight. $a_{11}$ is the evaluation of the 1. hazard according to the 1. risk parameter.

$v_{11}=w_1{a}_{11}$ = 0.1445 * 0.0791 = 0.0114

Step 12. Positive and negative ideal solution values are calculated. The maximum value of each column is called the positive ideal solution $\left(v_j^{+}\right)$, and the minimum value is called the negative ideal solution $\left(v_j^{-}\right)$.

$v_j^{+}$ = ( $v_1,v_2,\dots ,v_j)$ = (0.0298, 0.0088, 0.0240, 0.0072, 0.0289, 0.0345, 0.0084)

$v_j^{-}$ = ( $v_1,v_2,\dots ,v_j)$ = (0.0058, 0.0014, 0.0051, 0.0297, 0.0037, 0.0016, 0.0345).

Step 13. Positive and negative distance values are determined.

$$\begin{array}{rl}{S}^{+}& =\sqrt{\sum _{j=1}^7{\left(v_{ij}-v_j^{+}\right)}^2}\\ & =\sqrt{{\left(v_{11}-v_1^{+}\right)}^2+{\left(v_{12}-v_2^{+}\right)}^2+{\left(v_{13}-v_3^{+}\right)}^2+{\left(v_{14}-v_4^{+}\right)}^2+{\left(v_{15}-v_5^{+}\right)}^2+{\left(v_{16}-v_6^{+}\right)}^2+{\left(v_{17}-v_7^{+}\right)}^2}\\ & =\sqrt{\begin{array}{l}{\left(0.0003-0.0298\right)}^2+{\left(0.0000-0.0088\right)}^2+{\left(0.0002-0.0240\right)}^2+{\left(0.0002-0.0297\right)}^2+{\left(0.0002-0.0289\right)}^2\\ +{\left(0.0001-0.0345\right)}^2+{\left(0.0000-0.0345\right)}^2\end{array}}\\ & =0.0366\end{array}$$ $$\begin{array}{rl}{S}^{-}& =\sqrt{\sum _{j=1}^n{\left(v_{ij}-v_j^{-}\right)}^2}\\ & =\sqrt{{\left(v_{11}-v_1^{-}\right)}^2+{\left(v_{12}-v_2^{-}\right)}^2+{\left(v_{13}-v_3^{-}\right)}^2+{\left(v_{14}-v_4^{-}\right)}^2+{\left(v_{15}-v_5^{-}\right)}^2+{\left(v_{16}-v_6^{-}\right)}^2+{\left(v_{17}-v_7^{-}\right)}^2}\\ & =\sqrt{\begin{array}{l}{\left(0.0000-0.0058\right)}^2+{\left(0.0001-0.0014\right)}^2+{\left(0.0000-0.0051\right)}^2+{\left(0.0000-0.0072\right)}^2+{\left(0.0001-0.0037\right)}^2\\ +{\left(0.0006-0.0016\right)}^2+{\left(0.0005-0.0084\right)}^2\end{array}}\\ & =0.0324.\end{array}$$

Step 14. The closeness coefficient is calculated according to the ideal solution. Then the closeness coefficient values are normalized.

$$\begin{array}{rl}{C}_1=& {\displaystyle \frac{S_1^{-}}{S_1^{-}+S_1^{+}}=\frac{0.0324}{0.0324+0.0366}=0.4694}\\ C_1^{\ast }=& {\displaystyle \frac{C_1}{{\sum }_{i=1}^7{C}_i}={\displaystyle \frac{0.4694}{4.3474}=\mathrm{0.1080.}}}\end{array}$$

Step 15. Ranked from largest to smallest according to the normalized closeness coefficient. In this study, normalized closeness is expressed as risk coefficients of hazards. The MATLAB code for the calculations of the TOPSIS method is given in the File S10. Accordingly, the ranking of the hazards with the highest risk coefficient value is as follows: Electricity (0.1559) > Asphyxiating and Toxic Gases (0.1432) > Hot Water Use (0.1379) > Lifting and Transportation Operations (0.1027) > Moving Parts (0.1016) > Hot Surface (0.0959) > Working at Height (0.0941) > Cutting Tool Use (0.0879) > Taking Mould Manually (0.0808).

Sensitivity analysis

A sensitivity analysis is performed to assess how decision-making processes and outcomes react to uncertainties and changes. A sensitivity analysis in this study is performed by changing the risk parameter weights. With sensitivity analysis, the scenarios in the ranking results are evaluated depending on the changes in the weights of the risk parameters. Scenarios for changes in risk parameter weights and their explanations are presented in Table 8. The scenario results and rankings are shown in Table 9. In Table 9, it can be observed that when the risk parameter weights change, the risk coefficient values of the results change. Therefore, the proposed approach is sensitive to the weights of risk parameters. Sensitivity analysis of the proposed methodology shows that changing the risk parameter weights can have an impact on the final results, as expected. However, despite the variation of the risk coefficient values, H3, H4 and H5 are ranked as the most important hazards depending on all situations and are less susceptible to different weightings. This consistency indicates the robustness of the riskiest hazards. In Cases 4 and 5, the order of the nine hazards remains the same, except for minor differences. Differences in rankings for different risk parameter weights indicate the sensitivity of the proposed model to changes. Thus, the applicability and reliability of the proposed model to changing conditions are proven.

Comparative analysis

In this section, the risk coefficient ranking obtained by the IVFFAHP-TOPSIS method is compared with the ranking results of the simple additive weighting (SAW), multi-objective optimization by ratio analysis (MOORA) and evaluation based on distance from average solution (EDAS) methods integrated into the IVFFAHP method. The risk coefficient of the hazards obtained through the software and rankings are given in Table 10. In addition, pearson correlation coefficient values are calculated to examine the relationship between the ranking results of different MCDM methods, and these values are given in Table 10. When the rankings are analysed, it is seen that the highest-risk and lowest-risk hazards are in similar rankings. In TOPSIS, MOORA and EDA S methods, the hazard with the highest risk is determined as electricity (H3). However, the lowest risk hazard differs according to the methods. As can be seen in Table 10, the order of the riskiest hazards is like each other, although there are slight differences in the risk coefficient results between the methods. Similar results obtained because of applying the proposed method to different MCDM methods show the effectiveness and consistency of the proposed method. When the pearson corelation coefficients are examined, it is shown that the highest similarity rate is between the TOPSIS method and the MOORA method (r = 0.8422). There is a strong positive correlation between all MCDM methods, and it is concluded that IVFFAHP-TOPSIS shows a high degree of positive correlation with the results of other methods. The results of the pearson correlation test confirm that the ranking results of the TOPSIS method are consistent with the ranking results of other methods.