# Developing a Model for Safety Risk Assessment under Uncertainty for the Manufacturing Industry: A Case Study of Pole Factory Hazards in Riyadh, Saudi Arabia

Authors: Aref Ahmed; Mohammed Alkahtani; Abdulaziz M. El-Tamimi; Husam Kaid; Mustufa Haider Abidi
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Many occupational injuries occur in the manufacturing industry due to hazardous events. The available studies and statistics on occupational safety in the Kingdom of Saudi Arabia demonstrate the need for improving the work environment by introducing effective techniques for analyzing and assessing safety risks to control the most hazardous events. This study aims to develop a general model for assessing safety risks by integrating Monte Carlo simulation (MCS) and fuzzy set theory (FST) to overcome the uncertainty and unavailability of data on the severity and likelihood of hazards. MCS uses theModelRisk software for modeling hazards that exhibit randomness and uncertainty and have historical data. In contrast, FST uses a Matlab code to assess expert judgment about hazards featuring epistemic uncertainty or unavailable historical data. The Al-Babtain Pole Factory in Riyadh was selected as a case study in the manufacturing environment to prove the applicability and effectiveness of the developed model. From the 371 hazards identified using the Occupational Health and Safety Assessment Series 18001, only five were analyzed using the two model techniques. The likelihood and severity of these five hazards were collected and analyzed to obtain the risk levels. A list of hazards and their processing priorities were then produced. According to the risk values calculated using both techniques, Hazard5 was found to be the most hazardous event, followed by Hazard1. The results of the proposed model demonstrated the distributions, statistics, percentiles, and risk limits for the selected hazards. These outputs support decision-making and increase the effectiveness and flexibility of safety risk assessments, which means that the proposed model is reliable and applicable for SRA under uncertainty and data unavailability in the manufacturing industry.


---

## Body

## 1. Introduction

The industrial and social development at a global scale has created a vast demand for energy, chemicals, commodities, and food, which has increased the size and complexity of processing plants. This has advanced the development of new hazards and increased risks, which must be confronted with more than just economic benefits. Today, organizations are much concerned about their occupational health and safety (OHS) at work. By assessing their risks and improving their OHS, firms are trying to achieve the best practices in this field. Many companies have realized that a good OHS system must include risk assessment owing to its benefits, such as avoiding penalties and problems with the law, increasing personal productivity, improving the company’s image, and increasing profits.Risk assessment is a structured process that identifies both the probability and extent of adverse consequences arising from a given activity, facility, or system [1]. Glade et al. calculated risk using the following formula: risk = consequences of an undesirable event × the likelihood of that event [2]. Safety risk assessment (SRA) is a concept used in all industries, ranging from information technology, mining, and construction to manufacturing. The risks involved in industries are related to various factors, such as the physical and chemical properties of the input materials; propagation of incidents that have occurred; likelihood of incidents; and activities, locations, and effects of incidents on humans and the environment. Each industry has developed or adopted different health and SRA methodologies, such as hazard and operability studies, failure modes and effects analysis, and fault and event tree techniques, to reduce injuries at the workplace by setting priorities for dealing with high-, medium-, and low-risk incidents [3].So far, various techniques have been used to analyze SRA. Marhavilas et al. reviewed 404 articles on SRA techniques during 2000–2009 and determined that “quantitative” methods were used with the highest relative frequency (quantitative methods were used 65.63% times, qualitative methods were used 27.68% times, and hybrid methods were used 6.70% times) [4]. Quantitative techniques are very popular in reality, most of which have the same value for the probability of a hazard and its corresponding magnitude; thus, the product of these two input values yields the output risk value associated with that hazard. Pasman et al. discussed problems related to the present status of different risk assessment techniques for particular hazards and observed considerable variability in their outputs. They also provided some information about the variability sources that gave rise to large differences in the outputs [5].Due to these variations, most traditional techniques are ineffective in some cases; it is impossible to identify a single approach for an effective risk analysis that is powerful in all situations [6, 7]. Therefore, many authors have recently focused on developing or combining more than one technique to deal with the nature of the hazards that are characteristics of the manufacturing environment and have discovered that the results are more effective than using traditional techniques.Marhavilas studied risk assessment by developing an alternative risk assessment framework that included deterministic and stochastic processes and was applied to a power provider company for its risk assessment [8, 9]. This method is considered an improvement in risk assessment techniques because it can handle both stochastic and deterministic data. Yazdi et al. developed a new approach based on fault trees, fuzzy set theory (FST), and a Bayesian network to perform risk analysis more consistently under uncertain conditions and the statistical dependency of events [10]. The study was performed to overcome the challenges of unavailability of failure data, dependency of failure events, and data uncertainty. Li et al. proposed a model using FST and evidence theory to handle data uncertainty and utilized a Bayesian network to address the model uncertainty [11]. Seker used fuzzy logic as an assistance tool for conducting a bow-tie analysis to process the concept of uncertainty originating from different experts’ evaluations of hazard identification and risk assessment in the manufacturing industry [12]. The related studies conducted over the last 10 years have indicated a need to consider uncertainty in SRA to obtain the actual assessment of the risk of hazards represented by equivalent techniques.In the Kingdom of Saudi Arabia (KSA), health programs are being used to adopt the best practices for employees in organizations. Ibrahim et al. [13] stated that Saudi occupational health systems have significantly developed in recent decades, by focusing on health interventions to ensure that workers have suitable work-related health support and by assisting them in adapting to the physical and mental demands of their workplaces. Although the attention from the government and occupational health systems in the KSA might have improved the workplace safety and decreased the number of injuries over the last year, as reported by the General Organization for Social Insurance (GOSI) (2010–2018), the number of occupational injuries and rate of injury should still not be underestimated.According to GOSI (2018), 36,855 occupational injuries occurred in the KSA in 2018, of which 291 were fatal, which represents approximately 0.79% of all work-related injuries; 3,167 (8.59%) were permanent injuries; and 9,411 (25.53%) were still under medical care when the report was taken. Eighteen percent of all occupational injuries in 2018 occurred in the manufacturing industry, which ranked second after the construction industry. However, according to the annual injury rate indicator, which calculates a frequency measure to compare the hazardous areas, the manufacturing industry ranks the highest (785.4 injuries per 100,000 employees); moreover, the Riyadh region has the highest ratio (20.08%) in terms of occupational injuries (GOSI 2018).Yucesan et al. [14] applied the Pythagorean fuzzy analytical hierarchy process to perform risk assessment for hydroelectric power plants. In that study, experts identified 20 operational hazards, and the linguistic expressions technique was used to determine the weights. Based on the results, three key hazards were identified, for which preventive measures were proposed. Ilbahar et al. [15] proposed a novel framework and named it as Pythagorean fuzzy proportional risk assessment (PFPRA) for occupational safety and risk assessment. It integrated the Fine–Kinney, Pythagorean fuzzy analytic hierarchy process, and a fuzzy inference system. Excavation process in a construction site was presented as a case study on which the developed framework was employed. The results were compared against Pythagorean fuzzy failure modes and effects analysis, and it was reported that the proposed method provided better results in terms of reducing vagueness. Yucesan and Gul [16] proposed a novel technique, called neutrosophic best and worst method, and conducted a case study in an implant manufacturing company. The proposed technique was compared against classical and fuzzy BWM, and the results indicated that the proposed model can reflect uncertainty and ambiguity in real-world problems better than BWM. Lo et al. [17] proposed an integrated risk assessment model for manufacturing organization that combines several techniques (DEMTAL and TOPSIS) to produce a failure mode effect analysis (FMEA) model for comprehensive failure mode ranking generation. Chang et al. [18] presented a hybrid FMEA and multiattribute decision-making (MADM) model to develop an evaluation framework for risk assessment. They combined rough best-worst method (R-BWM) and rough technique for order preference by similarity to an ideal solution technique (R-TOPSIS) to determine the improvement order of failure modes. Moreover, they added the concept of aspiration level to R-TOPSIS technology to make it more effective and reliable. The results revealed that the proposed method can significantly improve the FMEA method for risk and reliability analysis.Tarafdar and Sinha [19] performed a health risk assessment caused by polycyclic aromatic hydrocarbons present in roadside soil in heavy mining area. Monte Carlo simulation technique was applied to study the influence of different parameters, and it was reported that skin adherence factor for soil was the most significant factor followed by exposure duration. Risk assessment is also applied in several fields such as cybersecurity risk assessment [20]. Researchers utilized multicriteria decision analysis method to analyze risk assessment in cybersecurity domain. A case study was used for evaluating the efficacy of the developed framework.Emerging risk management for Industry 4.0 was introduced by a previous study [21]. Various existing methods were discussed to present their applicability in a complex environment, such as Industry 4.0, where considerable attention was paid to human-machine and human-robot interactions. It was concluded that although Industry 4.0 offers numerous benefits, it also involves several risks and hazards, which require further attention. Badri et al. [22] also discussed the impact of Industry 4.0 on addressing OHS issues and concluded that, with the rapid changes in the manufacturing paradigm, the dangers will burgeon and the net impact on OHS will be negative.In addition, a study was conducted to assess OHS in a healthcare environment [23]. A two-stage fuzzy multicriteria approach was implemented, and a case study was presented for a hospital in Turkey. Furthermore, the fuzzy Vikor approach was used for prioritizing the hazard types in each department of the hospital. In another study, the fuzzy-AHP and fuzzy Vikor techniques were applied in the complex environment of a marine port for improvement in OHS [24].However, limited studies have been conducted in the field of SRA in the KSA, especially for the manufacturing industry [25]. The available studies have determined the presence of various hazards in the manufacturing environment that needs to be assessed. In addition, these studies used traditional methods [26]. Balgheeth demonstrated the need for academic and practical assessment of the current environment in the KSA to improve the outcomes [27].Using global research, as well as the safety situation in the KSA, this study examines the need to perform SRA by developing a general model that enhances the development of new techniques that can provide more effective results in assessing the hazards and accounts for data uncertainty in the risk assessment processes. The proposed methodology considerably differs from the existing approaches in that it uses a new uncertainty propagation algorithm that accounts for both data unavailability and type of uncertainty for the hazardous events in SRA. As a case study, the applicability and effectiveness of the proposed model are tested in a pole factory (PF) in the manufacturing environment in Riyadh, Saudi Arabia.The rest of this paper is organized as follows. Section2 presents the concept of uncertainty in safety risk assessment, Section 3 discusses the modeling techniques applied in this study, Section 4 presents details of the proposed methodology, Section 5 provides details of the application of the proposed methodology, and finally, Section 6 concludes this study.


## 2. Uncertainty in SRA

The impact of uncertainty on the results of risk analyses has been extensively studied. In a benchmarking exercise conducted on a major hazard analysis for a chemical plant managed by the joint research center during 1988–1990, 11 teams from different European countries analyzed a reference object—an ammonia storage facility [28]. This study aimed to estimate the degree of uncertainty in risk assessment. The results demonstrated considerable variability in risk estimations among the various analysis teams. A follow-up benchmarking exercise, which included an assessment of uncertainty in risk analyses of chemical establishments, was completed in late 2001. Seven teams from various European countries collaborated on this project. Each joined a team and, based on its input data and expertise, performed a risk analysis on the same risk scenarios in an ammonia storage facility. The results illustrated the considerable deviation between the assessed frequencies and the consequences of hazards. The spread in the results can be attributed to uncertainty. However, this uncertainty will be transferred to the final risk evaluation and, ultimately, affects decision-making. For example, one risk scenario might be acceptable for one team but not for another. Therefore, there is a significant need to analyze uncertainty and include it in SRA. Some critical considerations of the description and representation of uncertainty in risk assessments can be found in the literature [29].In general, uncertainty in risk assessments can be attributed to three main sources—model uncertainty, completeness uncertainty, and parameter uncertainty—which are described as follows [30].Model uncertainty arises because risk models are based on simple mathematical equations that represent reality but cannot completely characterize the complex physical processes or predict the output of a model that considers the risk level in SRA [31].Completeness uncertainty is related to risk contributors that are not addressed during the analytical process. The causes of this source of uncertainty are either known but not included in the analysis or not known and, therefore, not addressed in the analysis, which can significantly impact the analysis output.Parameter or data uncertainty refers to the uncertainty in the values of the parameters or input variables for a model, given that the mathematical form of that model has been agreed to be appropriate. A risk assessment must be capable of building a certain level of confidence in its results to support the decision-making process. This level of confidence, or accuracy, is difficult to achieve when the analytical inputs suffer from imprecision, lack of information, vagueness, and variability. In this case, handling parameters or variable data uncertainty is the most important and expensive part of the analysis if confidence in the results and accurate risk predictions are to be achieved [6, 7].According to Abrahamson [16], the causes of parameter uncertainty can be aleatory or epistemic. For aleatory parameter uncertainty, natural variability in the parameter values can be described using probability distributions. For epistemic parameter uncertainty, the analyst’s knowledge about the parameter values is described using probabilistic and nonprobabilistic distributions [32]. Each cause of parameter uncertainty requires an appropriate technique for presentation, modeling, and treatment. Abdo et al. presented the technique for treating each parameter according to the uncertainty type; its nature and the available information can help in obtaining more realistic results for risk estimation [7]. Other authors (e.g., Blair et al. and Darbra et al.) have proven that uncertainty manipulation using individual approaches is not feasible in real-world situations in which different informational qualities exist for various parameters [33, 34]. This is because some situations require more than one technique from those described previously to be more realistic, such as assigning probability distributions to certain parameters, while fuzzy-logic numbers would represent others. This new approach has led to the use of FST in combination with probabilistic methods to generate a hybrid approach for risk assessment [35]. Arunraj et al. developed a model based on FST and Monte Carlo simulation (MCS) for conducting analytical uncertainty in risk assessments [36]. Their proposed models provided a better measure of uncertainty than traditional risk analysis methods because they incorporate different types of information uncertainty into their risk calculations.


## 3. Modeling Techniques

Our developed model is based on the integration of a probabilistic approach and a fuzzy approach to incorporate the processing of data uncertainty when assessing risk in the manufacturing industry. The probabilistic approach uses MCS to process data randomness when the data are available. A fuzzy approach is used to process imprecision uncertainty for variable information when the data are rare or unavailable.


### 3.1. MCS

MCS is a probabilistic method that uses a computer simulation to combine multiple probability distributions into a risk equation. If the likelihood and severity consequences of an undesired event features uncertain variables and a probabilistic approach is suitable for presenting and propagating these variables, the most commonly used technique for the probabilistic approach is MCS, which assumes that the model variables (e.g., the likelihood and severity) are random and can be represented by a probability function [37].MCS is used to specify probability distributions for inputs and analyze the output. In this study, theModelRisk software was used to perform MCS and obtain and fit probability distributions for the likelihood and severity of hazards using historical data and produce input distributions to create probability distributions for the outputs [38]. After all random variables as inputs are shaped by probability distributions, the propagation of uncertainty for the distributions is calculated using MCS [39]. For example, let us consider a model whose output is a function Z = f(Y) = f (X1, X2, ..., Xn) of n uncertain variables, which are represented by probability distributions. The uncertainty calculated using MCS propagates through the following steps [6]: (a)
eti = 1 and N = number of samples.(b)
Take samplei, yi=x1,x2,x3,…xn, where xn, i=0,…,n are randomly generated from each distribution.(c)
Calculatezi=fyi=fx1,x2,x3,…,xn.(d)
If any sample is left, go back to stepb; otherwise, go to step e.(e)
Generate theprobability or mass distribution function (PDF) of the results (Z) by substituting all values of zi into a histogram and the cumulative distribution function of Z as Fz=1/N∑i=1NHi, where Hi=1 if zi<z and 0 elsewhere.


### 3.2. FST

FST was proposed by Zadeh in 1965 and applied for the assessment of expert systems under uncertainty, which have been widely used to describe uncertainty that cannot be properly presented using the probability theory or to handle epistemic data uncertainty [40]. In this approach, both the likelihood and consequences of an undesired event are assumed to be fuzzily obtained according to the expert judgment [10, 41]. First, the expert opinions regarding the magnitude of the probability and severity for the same hazard are taken qualitatively and transformed into a fuzzy number quantitatively. Then, their opinions on probability and severity are aggregated separately. Finally, the risk values are obtained according to the product of the aggregated probability and severity for each hazard using the following algorithms.


#### 3.2.1. Algorithm to Calculate the Aggregated Fuzzy Probability for Each Hazard

The number of experts or decision-makers (DMs) that assess the fuzzy probability of hazardj during a time period month can be computed as (1)Pˇij=pij−cij,pij+dij,wherei=1,2,…nexperts;For eachj=1,2,3,…,Nhazards.The midvalue of the fuzzy number for the probability of hazardj by expert i is pij, and the minimum value is estimated as (2)pij−cij.The maximum value is estimated as(3)pij+dij.Thus, the aggregated fuzzy probability for hazardj is obtained as follows [42]: (4)P∧Aj=Pj−Cj,Pj,Pj+Dj,where Cj=1/n∑i=1nci and Dj=1/n∑i=1ndi for each j=1,2,…,N.(5)Pj=min1<i<npi+max1<i<npi2.


#### 3.2.2. Algorithm to Calculate the Aggregated Fuzzy Severity for Each Hazard

The number of experts that assess the fuzzy severity of a hazardj can be computed as(6)Sˇij=sij−eij,sij,sij+fij,where i = 1, 2, …, n, for each j = 1, 2, 3, …, N.The midvalue of the fuzzy number required for calculating the severity of hazardj by expert i is sij, and the minimum value is estimated as (7)sij−eij.The maximum value is estimated as(8)sij+fij.Thus, the aggregated fuzzy severity for hazardj is obtained as follows [42]: (9)SˇAj=Sj−Ej,Sj,Sj+Fj,where Ej=1/n∑i=1nei and Fj=1/n∑i=1nfi.(10)Sj=min1<i<nsi+max1<i<nsi2,for eachj=1,2,…,N.


#### 3.2.3. Multiplication of the Probability and Severity for Each Hazard

The multiplication of the two fuzzy numbers for probability and severity, which are the products of the results obtained for the previous algorithms in equations (4) and (9), can be expressed as follows: (11)PˇAj=Pj−Cj,Pj,Pj+Dj,SˇAj=Sj−Ej,Sj,Sj+Fj.For easier processing, the following is performed:(12)PˇAj=Pj−Cj,Pj,Pj+Dj=a1j,a2j,a3j,which noted byA˜x,(13)SˇAj=Sj−Ej,Sj,Sj+Fj=b1j,b2j,b3j,which noted byB˜x.Therefore, the product of the two fuzzy membership functionsA˜x and B˜x can be expressed as(14)μA˜∗B˜x=−D1+D12+x−P/T11/2,P≤x≤Q,−D1+D22+x−R/U11/2,Q≤x≤R,0,otherwise.For each hazardj, the fuzzy probability and severity are expressed as A=a1,a2,a3 and B=b1,b2,b3, respectively, and represented as A∗B, which can be denoted as follows: (15)whereT1=a2−a1b2−b1,T2=a1a2−a1+ab2b2−b1,U1=a2−a1b2−b1,U2=b3a2−a1+a3b2−b1,D1=T1T2,D2=−U2U1,P=a1b1,Q=a2b2,R=a3b3.Evidently, the resulting fuzzy numberA˜∗B˜ is not a triangular fuzzy number. However, in most cases, computation with the resulting fuzzy numbers is highly tedious. Therefore, it is necessary to avoid the second- and higher-degree terms to simplify the computation. Therefore, the product of two fuzzy numbers is reduced to a triangular fuzzy number as P,Q,R=a1b1,a2b2,a3b3 [43].


## 3.1. MCS

MCS is a probabilistic method that uses a computer simulation to combine multiple probability distributions into a risk equation. If the likelihood and severity consequences of an undesired event features uncertain variables and a probabilistic approach is suitable for presenting and propagating these variables, the most commonly used technique for the probabilistic approach is MCS, which assumes that the model variables (e.g., the likelihood and severity) are random and can be represented by a probability function [37].MCS is used to specify probability distributions for inputs and analyze the output. In this study, theModelRisk software was used to perform MCS and obtain and fit probability distributions for the likelihood and severity of hazards using historical data and produce input distributions to create probability distributions for the outputs [38]. After all random variables as inputs are shaped by probability distributions, the propagation of uncertainty for the distributions is calculated using MCS [39]. For example, let us consider a model whose output is a function Z = f(Y) = f (X1, X2, ..., Xn) of n uncertain variables, which are represented by probability distributions. The uncertainty calculated using MCS propagates through the following steps [6]: (a)
eti = 1 and N = number of samples.(b)
Take samplei, yi=x1,x2,x3,…xn, where xn, i=0,…,n are randomly generated from each distribution.(c)
Calculatezi=fyi=fx1,x2,x3,…,xn.(d)
If any sample is left, go back to stepb; otherwise, go to step e.(e)
Generate theprobability or mass distribution function (PDF) of the results (Z) by substituting all values of zi into a histogram and the cumulative distribution function of Z as Fz=1/N∑i=1NHi, where Hi=1 if zi<z and 0 elsewhere.


## 3.2. FST

FST was proposed by Zadeh in 1965 and applied for the assessment of expert systems under uncertainty, which have been widely used to describe uncertainty that cannot be properly presented using the probability theory or to handle epistemic data uncertainty [40]. In this approach, both the likelihood and consequences of an undesired event are assumed to be fuzzily obtained according to the expert judgment [10, 41]. First, the expert opinions regarding the magnitude of the probability and severity for the same hazard are taken qualitatively and transformed into a fuzzy number quantitatively. Then, their opinions on probability and severity are aggregated separately. Finally, the risk values are obtained according to the product of the aggregated probability and severity for each hazard using the following algorithms.


### 3.2.1. Algorithm to Calculate the Aggregated Fuzzy Probability for Each Hazard

The number of experts or decision-makers (DMs) that assess the fuzzy probability of hazardj during a time period month can be computed as (1)Pˇij=pij−cij,pij+dij,wherei=1,2,…nexperts;For eachj=1,2,3,…,Nhazards.The midvalue of the fuzzy number for the probability of hazardj by expert i is pij, and the minimum value is estimated as (2)pij−cij.The maximum value is estimated as(3)pij+dij.Thus, the aggregated fuzzy probability for hazardj is obtained as follows [42]: (4)P∧Aj=Pj−Cj,Pj,Pj+Dj,where Cj=1/n∑i=1nci and Dj=1/n∑i=1ndi for each j=1,2,…,N.(5)Pj=min1<i<npi+max1<i<npi2.


### 3.2.2. Algorithm to Calculate the Aggregated Fuzzy Severity for Each Hazard

The number of experts that assess the fuzzy severity of a hazardj can be computed as(6)Sˇij=sij−eij,sij,sij+fij,where i = 1, 2, …, n, for each j = 1, 2, 3, …, N.The midvalue of the fuzzy number required for calculating the severity of hazardj by expert i is sij, and the minimum value is estimated as (7)sij−eij.The maximum value is estimated as(8)sij+fij.Thus, the aggregated fuzzy severity for hazardj is obtained as follows [42]: (9)SˇAj=Sj−Ej,Sj,Sj+Fj,where Ej=1/n∑i=1nei and Fj=1/n∑i=1nfi.(10)Sj=min1<i<nsi+max1<i<nsi2,for eachj=1,2,…,N.


### 3.2.3. Multiplication of the Probability and Severity for Each Hazard

The multiplication of the two fuzzy numbers for probability and severity, which are the products of the results obtained for the previous algorithms in equations (4) and (9), can be expressed as follows: (11)PˇAj=Pj−Cj,Pj,Pj+Dj,SˇAj=Sj−Ej,Sj,Sj+Fj.For easier processing, the following is performed:(12)PˇAj=Pj−Cj,Pj,Pj+Dj=a1j,a2j,a3j,which noted byA˜x,(13)SˇAj=Sj−Ej,Sj,Sj+Fj=b1j,b2j,b3j,which noted byB˜x.Therefore, the product of the two fuzzy membership functionsA˜x and B˜x can be expressed as(14)μA˜∗B˜x=−D1+D12+x−P/T11/2,P≤x≤Q,−D1+D22+x−R/U11/2,Q≤x≤R,0,otherwise.For each hazardj, the fuzzy probability and severity are expressed as A=a1,a2,a3 and B=b1,b2,b3, respectively, and represented as A∗B, which can be denoted as follows: (15)whereT1=a2−a1b2−b1,T2=a1a2−a1+ab2b2−b1,U1=a2−a1b2−b1,U2=b3a2−a1+a3b2−b1,D1=T1T2,D2=−U2U1,P=a1b1,Q=a2b2,R=a3b3.Evidently, the resulting fuzzy numberA˜∗B˜ is not a triangular fuzzy number. However, in most cases, computation with the resulting fuzzy numbers is highly tedious. Therefore, it is necessary to avoid the second- and higher-degree terms to simplify the computation. Therefore, the product of two fuzzy numbers is reduced to a triangular fuzzy number as P,Q,R=a1b1,a2b2,a3b3 [43].


## 3.2.1. Algorithm to Calculate the Aggregated Fuzzy Probability for Each Hazard

The number of experts or decision-makers (DMs) that assess the fuzzy probability of hazardj during a time period month can be computed as (1)Pˇij=pij−cij,pij+dij,wherei=1,2,…nexperts;For eachj=1,2,3,…,Nhazards.The midvalue of the fuzzy number for the probability of hazardj by expert i is pij, and the minimum value is estimated as (2)pij−cij.The maximum value is estimated as(3)pij+dij.Thus, the aggregated fuzzy probability for hazardj is obtained as follows [42]: (4)P∧Aj=Pj−Cj,Pj,Pj+Dj,where Cj=1/n∑i=1nci and Dj=1/n∑i=1ndi for each j=1,2,…,N.(5)Pj=min1<i<npi+max1<i<npi2.


## 3.2.2. Algorithm to Calculate the Aggregated Fuzzy Severity for Each Hazard

The number of experts that assess the fuzzy severity of a hazardj can be computed as(6)Sˇij=sij−eij,sij,sij+fij,where i = 1, 2, …, n, for each j = 1, 2, 3, …, N.The midvalue of the fuzzy number required for calculating the severity of hazardj by expert i is sij, and the minimum value is estimated as (7)sij−eij.The maximum value is estimated as(8)sij+fij.Thus, the aggregated fuzzy severity for hazardj is obtained as follows [42]: (9)SˇAj=Sj−Ej,Sj,Sj+Fj,where Ej=1/n∑i=1nei and Fj=1/n∑i=1nfi.(10)Sj=min1<i<nsi+max1<i<nsi2,for eachj=1,2,…,N.


## 3.2.3. Multiplication of the Probability and Severity for Each Hazard

The multiplication of the two fuzzy numbers for probability and severity, which are the products of the results obtained for the previous algorithms in equations (4) and (9), can be expressed as follows: (11)PˇAj=Pj−Cj,Pj,Pj+Dj,SˇAj=Sj−Ej,Sj,Sj+Fj.For easier processing, the following is performed:(12)PˇAj=Pj−Cj,Pj,Pj+Dj=a1j,a2j,a3j,which noted byA˜x,(13)SˇAj=Sj−Ej,Sj,Sj+Fj=b1j,b2j,b3j,which noted byB˜x.Therefore, the product of the two fuzzy membership functionsA˜x and B˜x can be expressed as(14)μA˜∗B˜x=−D1+D12+x−P/T11/2,P≤x≤Q,−D1+D22+x−R/U11/2,Q≤x≤R,0,otherwise.For each hazardj, the fuzzy probability and severity are expressed as A=a1,a2,a3 and B=b1,b2,b3, respectively, and represented as A∗B, which can be denoted as follows: (15)whereT1=a2−a1b2−b1,T2=a1a2−a1+ab2b2−b1,U1=a2−a1b2−b1,U2=b3a2−a1+a3b2−b1,D1=T1T2,D2=−U2U1,P=a1b1,Q=a2b2,R=a3b3.Evidently, the resulting fuzzy numberA˜∗B˜ is not a triangular fuzzy number. However, in most cases, computation with the resulting fuzzy numbers is highly tedious. Therefore, it is necessary to avoid the second- and higher-degree terms to simplify the computation. Therefore, the product of two fuzzy numbers is reduced to a triangular fuzzy number as P,Q,R=a1b1,a2b2,a3b3 [43].


## 4. Proposed Methodology

The proposed model for SRA (Figure1) combines MCS and FST for an effective risk assessment in the manufacturing industry. First, the possible hazards in a company are identified as a comprehensive list of possible hazards and classified according to the type of uncertainty and availability of data for the input variables for each hazard, to assign an appropriate approach for presentation and aggregation to produce the output (i.e., the level and distribution of risk). If uncertainty in the input parameters arises from the randomness caused by the inherent variability, these input parameters can be categorized as those with aleatory uncertainty, such as event frequencies estimated using sufficient statistical data, and presented as probability distributions [44]. In case where the available information is very scarce, even if one adopts the elicitation of expert knowledge to incorporate the diffused information, possibility distributions are used to model epistemic uncertainty [45]. The advantage of MCS and FST is that the analysis of uncertainty is not separated from the actual risk calculations but is built into the model calculations and can be turned on or off depending on how the input variables are defined. The proposed model uses software and mathematical equations, and it can be enhanced using software codes to facilitate use for any hazard in various manufacturing industries. To the best of our knowledge, this is the first study to have provided a clear picture for analyzing the risk with uncertainty. The following are the advantages of integrating MCS and FST in the proposed model for SRA:(i)
Perform SRA effectively and comprehensively by considering uncertainty and data availability for hazards, irrespective of whether the industries have historical data or missing, or even, unknown information regarding these hazards. MCS is used to analyze the hazards that have randomized historical data, whereas FST is used to analyze hazards that do not have historical data and require expert opinion as well as to overcome the limitation of the MCS regarding nonrandom data.(ii)
The results of the proposed model indicate the distributions, statistics, percentiles, and limits of risks for the ranking list of hazards and lead to an effective decision-making in SRA.Figure 1
Proposed methodology for SRA under uncertainty.


## 5. Methodology Application

### 5.1. Case Study

To demonstrate and verify the applicability of the developed model for risk assessment, a PF in ABPT Company, which was founded in 1955, was selected as a case study in the manufacturing environment of Riyadh, Saudi Arabia. This company is one of the largest enterprises in the engineering sector and the manufacturing industry, and according to the 2019 GOSI report and based on the information supplied by the company’s safety administration, it contributes to 3% of the total hazards in Riyadh. According to the historical records of the company and the injuries and man-days lost (MDL) reported by the administration, the trend for the severity of PF injuries had a sharp 4.7 times increase in the second quarter of 2019 compared to the first quarter. For these reasons, we selected PF as the case study to apply the developed model.Figure2 shows the total MDL as a measure of the impact or severity of large injuries that occurred in factories of the company during 2013–2019. Evidently, throughout the study period, the MDL in the PF maintained a high value, where the highest value occurred during 2019 and tended to increase beyond this point.Figure 2
Trend of severity for PF relative to others factories in ABPT Company (source: ABPT Company).


### 5.2. Hazard Identification

The PF was certified with the Occupational Health and Safety Assessment Series (OHSAS) 18001 for the occupational health and safety management system, which requires an analysis of the root causes of the hazards and accidents. OHSAS 18001 facilitates a list of hazards and historical data for the PF. Therefore, the information mainly comprised administrative meetings and reviews of related company documents, including the incident statistics and risk assessment histories in the company. The hazard identification stage aims to produce a comprehensive list of all possible hazardous events that have led to accidents or near misses in the company.The hazard records indicate 371 hazards in the production section of PF, and it is difficult to collect the data for frequency and severity of all these hazards. As the purpose here was to apply the model, five common hazards were selected to be analyzed and applied by utilizing the two approaches in the proposed model. The main reason behind selecting only five hazards is that the methodology to be applied for these hazards is almost the same as that for the other types of hazards because the model can handle hazards that may or may not have historical data. In addition, many hazards did not have enough data regarding the last three years because they only occurred for a short period or had already been eliminated or controlled. Therefore, after discussing and meeting with the safety management team in the company, five hazardous events that possessed historical data for the last three years were applied by MCS (other hazards did have not enough data), and expert opinions on these hazards were used to enable the application of the model and comparison of the results of MCS and FST.The specific hazards included the falling of materials due to overstacking and improper arrangement (Hazard1), exposure to grinding flying particles (Hazard2), trip and fall while skimming due to poor housekeeping (Hazard3), falling of materials due to lifting beyond the capacity (Hazard4), and forklift accidents caused by speeding (Hazard5).


### 5.3. Data Collection

After identifying the specific hazards occurring inside the plant, the frequency and consequences of the potential hazardous events were determined (only the consequences and severity on the human body were considered). To perform risk assessment for any hazard, its frequency and severity must be known. In this study, the following two types of data were collected depending on the techniques used in the developed model.


#### 5.3.1. Historical Data from the Company’s Records

The data on the frequencies and magnitudes of severity for the five identified hazardous events over the last three years (2017–2019) were gathered from the company’s records and then transformed in a standard manner, as shown in SupplementaryA.


#### 5.3.2. Qualitative Data according to Experts’ Judgment

For the same studied hazards that were gathered using historical data, it was supposed that the hazards had no historical data and that experts’ knowledge was required to assess their likelihood and severity, which would be processed using FST.Experts prefer to use linguistic terms rather than numerical expressions to justify the severity and likelihood of hazardous events. The linguistic terms for the likelihood are scaled according to the following five levels referring to Aqlan et al. [42] and Seker [12]: expected (E), possible (P), unlikely (UL), very unlikely (VL), and not expected (NE). In addition, the gates for severity are scaled as high (H), medium (M), low (L), very low (VL), and none (N). In practice, experts’ knowledge is a significant and available means for obtaining hazard-related data in the absence of crisp values or accurate PDFs. Therefore, in this study, we supposed that the specified hazards did not have distributions or historical data and required expert judgment. Five experts are employees in safety department in Al-Babtain Company were selected to provide their qualitative opinions regarding the magnitude of likelihood and severity for the selected hazards. The experts have approximately the same educational level with a diploma in safety, and they have more than five years of experience in field of safety; thus, we assumed that all experts had the same weights. A table was designed and given to experts independently to fill their opinions about the assessment of likelihood and severity of selected hazards as presented in Supplementary B.1, which was used by the one expert. Then, the opinions were quantitatively transformed into a fuzzy number, as shown in Table 1.Table 1
Linguistic expressions and their corresponding fuzzy numbers for likelihood and severity.
Linguistic assessment variables for likelihoodLinguistic assessment variables for severityCorresponding fuzzy numbersCharacteristic function of fuzzy numbers (min, mid, max)ExpectedHigh5(4,5,5)PossibleMedium4(3,4,5)UnlikelyLow3(2,3,4)Very unlikelyVery low2(1,2,3)NonexpectedNone1(1,1,2)The data collected for all experts and the corresponding fuzzy numbers are presented in SupplementaryB.2.


### 5.4. Implementation of the Proposed Model

After hazard identification and presentation of the collected data, the proposed methodologies for analyzing the historical data (i.e., MCS) and expert opinions (i.e., FST) were implemented. This is explained in the following sections.


#### 5.4.1. Using MCS

First, randomness tests were conducted for all collected data by running test tools in Minitab 17 to ensure that the data did not exhibit any trend or pattern, which indicates that the observed variations were special-cause variations. Next, the frequency and severity of the hazards were considered as stochastic variables, which were input variables for MCS, and then, the fitting distributions for each input variable were found to produce the output distributions of the risk for specific hazards. The randomness and distribution fittings are presented in Table2. Evidently, the frequency and severity for all hazards were randomly distributed, except for the frequency of Hazard4. Therefore, Hazard4 was excluded from the MCS analysis but was analyzed for the other hazards using FST.Table 2
Randomness test and fitting of the distributions.
HazardsFrequencyf hazardSeverity of hazardRandomness testFitting distributionRandomness testFitting distributionFalling of materials due to overstacking and improper arrangements (H1)P value = 0.499P value > 0.05Binomialn=5,p=0.7616P value = 0.476P value > 0.05Poissonλ=2.626Exposure to grinding flying particles (H2)P value = 0.735P value > 0.05Binomialn=5,p=0.72222P value = 0.423P value > 0.05Geometricp=0.35036Trip and fall while skimming due to poor housekeeping (H3)P value = 0.814P value > 0.05Binomialn=4,p=0.70714P value = 0.542P value > 0.05Binomialn=2,p=0.62201Falling of materials due to lifting over capacity (H4)P value = 0.011P value < 0.05It will exclude from the study for this technique)P value = 0.770P value > 0.05Binomialn=5,p=0.82942Forklift accidents due to speeding (H5)P value = 0.380P value > 0.05Binomialn=4,p=0.65789)P value = 0.122P value > 0.05Binomialn=4,p=0.85566After that, the PDF for the output function expressed the aggregation of the two input variables that were randomly selected from the frequency and severity distributions. This step was repeated 1,000 times, and the obtained 1,000 output values represent the probability distribution of the output function (i.e., risk distribution).The risk distribution forHazard1 will be extensively explained to demonstrate the importance of utilizing this technique for SRA by DMs. Figure 3 illustrates the distribution of the risk of Hazard1, which is dependent on the distribution of the frequency and severity. Figure 4 shows a box plot for the risk values related to frequency and severity within the range of quartile percentiles (between 25% and 75%).Figure 3
Results of aggregation (the risk) of the frequency and severity forHazard1.Figure 4
Box plot for the inputs and outputs forHazard1.The frequency of this hazardous event was higher and more influential than the severity of the risk value. Therefore, if this hazard takes priority for treatment and is difficult to eliminate, at first glance, the focus should be on minimizing its frequency rather than its severity.Furthermore, the sensitivity analysis forHazard1 obtained from the tornado plot shown in Figure 5 affirms that the frequency variable has more influence on the mean of risks for Hazard1. Although the relative importance of both frequency and severity for this hazard is approximate, it is useful to estimate the low, base, and high outcomes for each variable and, therefore, should be focused on by the DMs.Figure 5
Sensitivity analysis forHazard1.In addition, it is important to ensure the reliability and steady state of the simulation results. We ran the simulation process eight times for the sameHazard1 to compare the results of cumulative plots for all runs, as presented in Figure 5. The means and cumulative probability plots for all runs relative to the risk values were found to be identical. Therefore, the simulation conditions with one run of simulations (1,000 samples) are suitable and sufficient to be run and analyze Hazard1. Figure 6 presents considerable useful information about the risk of Hazard1, which includes important statistics about this risk and its possible risk values. This outcome enables us to determine the ranking of this hazard relative to other hazards in the company by comparing them statistically to obtain a list of hazards ranked according to their priority for processing or elimination.Figure 6
Probability plot for the risk ofHazard1.Figure6 shows additional details for Hazard1 that are important for managers or DMs in safety administration, such as the confidence level for their decisions related to determining the processing priority. Regarding the remaining specified hazards, their frequencies and severities are aggregated to produce the values of risk. The risks for the specified hazards are summarized in Table 3, where the risk values for each hazard are presented as percentiles—25%, 50%, 75%, and 95%. In addition, the specific hazards are ranked depending on this summary.Table 3
Comparison of the statistics on the risk of hazards.
HazardParametersRanking for priority to be taken25%Mean75%95%SDFalling of materials due to overstacking and improper arrangements (H1)811.2014173.902Exposure to grinding flying particles (H2)46.69154.223Trip and fall while skimming due to poor housekeeping (H3)45.5792.284Forklift accidents due to speeding (H5)1013.4417214.791Table3 presents a list of the specified hazards that were analyzed, where the priority for taking action to eliminate or minimize them depends on the mean value of the risk. The hazardous event that must be prioritized for elimination or minimization from the studied hazards is Hazard5, followed by Hazard1. The priority sequence for the studied hazards is as follows: H5, H1, H2, and H3. Regarding Hazard4 (H4), because its frequency is not random, it is not suitable to calculate its distribution before analyzing a special cause for the variation. Therefore, it is excluded from this technique.


#### 5.4.2. Using FST

The qualitative inputs gathered by the experts are presented in SupplementaryB.2 and are subjected to the mathematical equations and procedures presented in Section 3.2. To facilitate this process, a Matlab code was developed to read the experts’ judgments and perform calculations to obtain the aggregate of the frequencies and severities and then the fuzzy risks for the hazards. The results obtained for the code are presented in Table 4 and Figure 7, where the mean of the risk for the hazards is critical to rank the priority sequence for eliminating or controlling hazards as H5, H1, H4, H2, and H3. Therefore, mitigation plans related to these hazards should be implemented in this order.Table 4
Results of aggregation of experts’ evaluations using the fuzzy method.
Aggregate fuzzyHazard1Hazard2Hazard3Hazard4Hazard5MinMidMaxMinMidMaxMinMidMaxMinMidMaxMinMidMaxFrequency2.53.54.12.53.54.1344.82342.53.54.5Severity344.42341.52.53.5344.43.54.55.3Risks7.51418510.516.44.51016.8612188.81623.9Figure 7
Results for the risk of hazards calculated using the fuzzy method.Although this technique has been introduced to analyze hazards that have no historical data and require expert judgment, in this study, it was used for hazards that were already analyzed using MCS. Thus, the outcomes became more confident and FST was used to address the shortage with MCS regarding the data that were nonrandom, as in the case ofHazard4.


#### 5.4.3. Comparative Analysis of the Proposed Models

Comparing the priority hazard lists using MCS and FST, both techniques obtained almost similar sequences, where the results of the priority of hazards to be treated by MCS wereH5-H1-H2-H3 (with H4 removed from MCS) and that to be treated by FST were H5-H1-H4-H2-H3. Moreover, the risk values of the hazards obtained using the two techniques were similar. The results showed that the integration of the two techniques for analyzing the different hazards occurs by overcoming the limitation of MCS regarding nonrandom data, as in the case of Hazard4 or when the hazards did not possess the recorded data. However, MCS can better analyze the hazards that have recorded randomized data and when determining expert opinions for the assessment is difficult and costly.


## 5.1. Case Study

To demonstrate and verify the applicability of the developed model for risk assessment, a PF in ABPT Company, which was founded in 1955, was selected as a case study in the manufacturing environment of Riyadh, Saudi Arabia. This company is one of the largest enterprises in the engineering sector and the manufacturing industry, and according to the 2019 GOSI report and based on the information supplied by the company’s safety administration, it contributes to 3% of the total hazards in Riyadh. According to the historical records of the company and the injuries and man-days lost (MDL) reported by the administration, the trend for the severity of PF injuries had a sharp 4.7 times increase in the second quarter of 2019 compared to the first quarter. For these reasons, we selected PF as the case study to apply the developed model.Figure2 shows the total MDL as a measure of the impact or severity of large injuries that occurred in factories of the company during 2013–2019. Evidently, throughout the study period, the MDL in the PF maintained a high value, where the highest value occurred during 2019 and tended to increase beyond this point.Figure 2
Trend of severity for PF relative to others factories in ABPT Company (source: ABPT Company).


## 5.2. Hazard Identification

The PF was certified with the Occupational Health and Safety Assessment Series (OHSAS) 18001 for the occupational health and safety management system, which requires an analysis of the root causes of the hazards and accidents. OHSAS 18001 facilitates a list of hazards and historical data for the PF. Therefore, the information mainly comprised administrative meetings and reviews of related company documents, including the incident statistics and risk assessment histories in the company. The hazard identification stage aims to produce a comprehensive list of all possible hazardous events that have led to accidents or near misses in the company.The hazard records indicate 371 hazards in the production section of PF, and it is difficult to collect the data for frequency and severity of all these hazards. As the purpose here was to apply the model, five common hazards were selected to be analyzed and applied by utilizing the two approaches in the proposed model. The main reason behind selecting only five hazards is that the methodology to be applied for these hazards is almost the same as that for the other types of hazards because the model can handle hazards that may or may not have historical data. In addition, many hazards did not have enough data regarding the last three years because they only occurred for a short period or had already been eliminated or controlled. Therefore, after discussing and meeting with the safety management team in the company, five hazardous events that possessed historical data for the last three years were applied by MCS (other hazards did have not enough data), and expert opinions on these hazards were used to enable the application of the model and comparison of the results of MCS and FST.The specific hazards included the falling of materials due to overstacking and improper arrangement (Hazard1), exposure to grinding flying particles (Hazard2), trip and fall while skimming due to poor housekeeping (Hazard3), falling of materials due to lifting beyond the capacity (Hazard4), and forklift accidents caused by speeding (Hazard5).


## 5.3. Data Collection

After identifying the specific hazards occurring inside the plant, the frequency and consequences of the potential hazardous events were determined (only the consequences and severity on the human body were considered). To perform risk assessment for any hazard, its frequency and severity must be known. In this study, the following two types of data were collected depending on the techniques used in the developed model.


### 5.3.1. Historical Data from the Company’s Records

The data on the frequencies and magnitudes of severity for the five identified hazardous events over the last three years (2017–2019) were gathered from the company’s records and then transformed in a standard manner, as shown in SupplementaryA.


### 5.3.2. Qualitative Data according to Experts’ Judgment

For the same studied hazards that were gathered using historical data, it was supposed that the hazards had no historical data and that experts’ knowledge was required to assess their likelihood and severity, which would be processed using FST.Experts prefer to use linguistic terms rather than numerical expressions to justify the severity and likelihood of hazardous events. The linguistic terms for the likelihood are scaled according to the following five levels referring to Aqlan et al. [42] and Seker [12]: expected (E), possible (P), unlikely (UL), very unlikely (VL), and not expected (NE). In addition, the gates for severity are scaled as high (H), medium (M), low (L), very low (VL), and none (N). In practice, experts’ knowledge is a significant and available means for obtaining hazard-related data in the absence of crisp values or accurate PDFs. Therefore, in this study, we supposed that the specified hazards did not have distributions or historical data and required expert judgment. Five experts are employees in safety department in Al-Babtain Company were selected to provide their qualitative opinions regarding the magnitude of likelihood and severity for the selected hazards. The experts have approximately the same educational level with a diploma in safety, and they have more than five years of experience in field of safety; thus, we assumed that all experts had the same weights. A table was designed and given to experts independently to fill their opinions about the assessment of likelihood and severity of selected hazards as presented in Supplementary B.1, which was used by the one expert. Then, the opinions were quantitatively transformed into a fuzzy number, as shown in Table 1.Table 1
Linguistic expressions and their corresponding fuzzy numbers for likelihood and severity.
Linguistic assessment variables for likelihoodLinguistic assessment variables for severityCorresponding fuzzy numbersCharacteristic function of fuzzy numbers (min, mid, max)ExpectedHigh5(4,5,5)PossibleMedium4(3,4,5)UnlikelyLow3(2,3,4)Very unlikelyVery low2(1,2,3)NonexpectedNone1(1,1,2)The data collected for all experts and the corresponding fuzzy numbers are presented in SupplementaryB.2.


## 5.3.1. Historical Data from the Company’s Records

The data on the frequencies and magnitudes of severity for the five identified hazardous events over the last three years (2017–2019) were gathered from the company’s records and then transformed in a standard manner, as shown in SupplementaryA.


## 5.3.2. Qualitative Data according to Experts’ Judgment

For the same studied hazards that were gathered using historical data, it was supposed that the hazards had no historical data and that experts’ knowledge was required to assess their likelihood and severity, which would be processed using FST.Experts prefer to use linguistic terms rather than numerical expressions to justify the severity and likelihood of hazardous events. The linguistic terms for the likelihood are scaled according to the following five levels referring to Aqlan et al. [42] and Seker [12]: expected (E), possible (P), unlikely (UL), very unlikely (VL), and not expected (NE). In addition, the gates for severity are scaled as high (H), medium (M), low (L), very low (VL), and none (N). In practice, experts’ knowledge is a significant and available means for obtaining hazard-related data in the absence of crisp values or accurate PDFs. Therefore, in this study, we supposed that the specified hazards did not have distributions or historical data and required expert judgment. Five experts are employees in safety department in Al-Babtain Company were selected to provide their qualitative opinions regarding the magnitude of likelihood and severity for the selected hazards. The experts have approximately the same educational level with a diploma in safety, and they have more than five years of experience in field of safety; thus, we assumed that all experts had the same weights. A table was designed and given to experts independently to fill their opinions about the assessment of likelihood and severity of selected hazards as presented in Supplementary B.1, which was used by the one expert. Then, the opinions were quantitatively transformed into a fuzzy number, as shown in Table 1.Table 1
Linguistic expressions and their corresponding fuzzy numbers for likelihood and severity.
Linguistic assessment variables for likelihoodLinguistic assessment variables for severityCorresponding fuzzy numbersCharacteristic function of fuzzy numbers (min, mid, max)ExpectedHigh5(4,5,5)PossibleMedium4(3,4,5)UnlikelyLow3(2,3,4)Very unlikelyVery low2(1,2,3)NonexpectedNone1(1,1,2)The data collected for all experts and the corresponding fuzzy numbers are presented in SupplementaryB.2.


## 5.4. Implementation of the Proposed Model

After hazard identification and presentation of the collected data, the proposed methodologies for analyzing the historical data (i.e., MCS) and expert opinions (i.e., FST) were implemented. This is explained in the following sections.


### 5.4.1. Using MCS

First, randomness tests were conducted for all collected data by running test tools in Minitab 17 to ensure that the data did not exhibit any trend or pattern, which indicates that the observed variations were special-cause variations. Next, the frequency and severity of the hazards were considered as stochastic variables, which were input variables for MCS, and then, the fitting distributions for each input variable were found to produce the output distributions of the risk for specific hazards. The randomness and distribution fittings are presented in Table2. Evidently, the frequency and severity for all hazards were randomly distributed, except for the frequency of Hazard4. Therefore, Hazard4 was excluded from the MCS analysis but was analyzed for the other hazards using FST.Table 2
Randomness test and fitting of the distributions.
HazardsFrequencyf hazardSeverity of hazardRandomness testFitting distributionRandomness testFitting distributionFalling of materials due to overstacking and improper arrangements (H1)P value = 0.499P value > 0.05Binomialn=5,p=0.7616P value = 0.476P value > 0.05Poissonλ=2.626Exposure to grinding flying particles (H2)P value = 0.735P value > 0.05Binomialn=5,p=0.72222P value = 0.423P value > 0.05Geometricp=0.35036Trip and fall while skimming due to poor housekeeping (H3)P value = 0.814P value > 0.05Binomialn=4,p=0.70714P value = 0.542P value > 0.05Binomialn=2,p=0.62201Falling of materials due to lifting over capacity (H4)P value = 0.011P value < 0.05It will exclude from the study for this technique)P value = 0.770P value > 0.05Binomialn=5,p=0.82942Forklift accidents due to speeding (H5)P value = 0.380P value > 0.05Binomialn=4,p=0.65789)P value = 0.122P value > 0.05Binomialn=4,p=0.85566After that, the PDF for the output function expressed the aggregation of the two input variables that were randomly selected from the frequency and severity distributions. This step was repeated 1,000 times, and the obtained 1,000 output values represent the probability distribution of the output function (i.e., risk distribution).The risk distribution forHazard1 will be extensively explained to demonstrate the importance of utilizing this technique for SRA by DMs. Figure 3 illustrates the distribution of the risk of Hazard1, which is dependent on the distribution of the frequency and severity. Figure 4 shows a box plot for the risk values related to frequency and severity within the range of quartile percentiles (between 25% and 75%).Figure 3
Results of aggregation (the risk) of the frequency and severity forHazard1.Figure 4
Box plot for the inputs and outputs forHazard1.The frequency of this hazardous event was higher and more influential than the severity of the risk value. Therefore, if this hazard takes priority for treatment and is difficult to eliminate, at first glance, the focus should be on minimizing its frequency rather than its severity.Furthermore, the sensitivity analysis forHazard1 obtained from the tornado plot shown in Figure 5 affirms that the frequency variable has more influence on the mean of risks for Hazard1. Although the relative importance of both frequency and severity for this hazard is approximate, it is useful to estimate the low, base, and high outcomes for each variable and, therefore, should be focused on by the DMs.Figure 5
Sensitivity analysis forHazard1.In addition, it is important to ensure the reliability and steady state of the simulation results. We ran the simulation process eight times for the sameHazard1 to compare the results of cumulative plots for all runs, as presented in Figure 5. The means and cumulative probability plots for all runs relative to the risk values were found to be identical. Therefore, the simulation conditions with one run of simulations (1,000 samples) are suitable and sufficient to be run and analyze Hazard1. Figure 6 presents considerable useful information about the risk of Hazard1, which includes important statistics about this risk and its possible risk values. This outcome enables us to determine the ranking of this hazard relative to other hazards in the company by comparing them statistically to obtain a list of hazards ranked according to their priority for processing or elimination.Figure 6
Probability plot for the risk ofHazard1.Figure6 shows additional details for Hazard1 that are important for managers or DMs in safety administration, such as the confidence level for their decisions related to determining the processing priority. Regarding the remaining specified hazards, their frequencies and severities are aggregated to produce the values of risk. The risks for the specified hazards are summarized in Table 3, where the risk values for each hazard are presented as percentiles—25%, 50%, 75%, and 95%. In addition, the specific hazards are ranked depending on this summary.Table 3
Comparison of the statistics on the risk of hazards.
HazardParametersRanking for priority to be taken25%Mean75%95%SDFalling of materials due to overstacking and improper arrangements (H1)811.2014173.902Exposure to grinding flying particles (H2)46.69154.223Trip and fall while skimming due to poor housekeeping (H3)45.5792.284Forklift accidents due to speeding (H5)1013.4417214.791Table3 presents a list of the specified hazards that were analyzed, where the priority for taking action to eliminate or minimize them depends on the mean value of the risk. The hazardous event that must be prioritized for elimination or minimization from the studied hazards is Hazard5, followed by Hazard1. The priority sequence for the studied hazards is as follows: H5, H1, H2, and H3. Regarding Hazard4 (H4), because its frequency is not random, it is not suitable to calculate its distribution before analyzing a special cause for the variation. Therefore, it is excluded from this technique.


### 5.4.2. Using FST

The qualitative inputs gathered by the experts are presented in SupplementaryB.2 and are subjected to the mathematical equations and procedures presented in Section 3.2. To facilitate this process, a Matlab code was developed to read the experts’ judgments and perform calculations to obtain the aggregate of the frequencies and severities and then the fuzzy risks for the hazards. The results obtained for the code are presented in Table 4 and Figure 7, where the mean of the risk for the hazards is critical to rank the priority sequence for eliminating or controlling hazards as H5, H1, H4, H2, and H3. Therefore, mitigation plans related to these hazards should be implemented in this order.Table 4
Results of aggregation of experts’ evaluations using the fuzzy method.
Aggregate fuzzyHazard1Hazard2Hazard3Hazard4Hazard5MinMidMaxMinMidMaxMinMidMaxMinMidMaxMinMidMaxFrequency2.53.54.12.53.54.1344.82342.53.54.5Severity344.42341.52.53.5344.43.54.55.3Risks7.51418510.516.44.51016.8612188.81623.9Figure 7
Results for the risk of hazards calculated using the fuzzy method.Although this technique has been introduced to analyze hazards that have no historical data and require expert judgment, in this study, it was used for hazards that were already analyzed using MCS. Thus, the outcomes became more confident and FST was used to address the shortage with MCS regarding the data that were nonrandom, as in the case ofHazard4.


### 5.4.3. Comparative Analysis of the Proposed Models

Comparing the priority hazard lists using MCS and FST, both techniques obtained almost similar sequences, where the results of the priority of hazards to be treated by MCS wereH5-H1-H2-H3 (with H4 removed from MCS) and that to be treated by FST were H5-H1-H4-H2-H3. Moreover, the risk values of the hazards obtained using the two techniques were similar. The results showed that the integration of the two techniques for analyzing the different hazards occurs by overcoming the limitation of MCS regarding nonrandom data, as in the case of Hazard4 or when the hazards did not possess the recorded data. However, MCS can better analyze the hazards that have recorded randomized data and when determining expert opinions for the assessment is difficult and costly.


## 5.4.1. Using MCS

First, randomness tests were conducted for all collected data by running test tools in Minitab 17 to ensure that the data did not exhibit any trend or pattern, which indicates that the observed variations were special-cause variations. Next, the frequency and severity of the hazards were considered as stochastic variables, which were input variables for MCS, and then, the fitting distributions for each input variable were found to produce the output distributions of the risk for specific hazards. The randomness and distribution fittings are presented in Table2. Evidently, the frequency and severity for all hazards were randomly distributed, except for the frequency of Hazard4. Therefore, Hazard4 was excluded from the MCS analysis but was analyzed for the other hazards using FST.Table 2
Randomness test and fitting of the distributions.
HazardsFrequencyf hazardSeverity of hazardRandomness testFitting distributionRandomness testFitting distributionFalling of materials due to overstacking and improper arrangements (H1)P value = 0.499P value > 0.05Binomialn=5,p=0.7616P value = 0.476P value > 0.05Poissonλ=2.626Exposure to grinding flying particles (H2)P value = 0.735P value > 0.05Binomialn=5,p=0.72222P value = 0.423P value > 0.05Geometricp=0.35036Trip and fall while skimming due to poor housekeeping (H3)P value = 0.814P value > 0.05Binomialn=4,p=0.70714P value = 0.542P value > 0.05Binomialn=2,p=0.62201Falling of materials due to lifting over capacity (H4)P value = 0.011P value < 0.05It will exclude from the study for this technique)P value = 0.770P value > 0.05Binomialn=5,p=0.82942Forklift accidents due to speeding (H5)P value = 0.380P value > 0.05Binomialn=4,p=0.65789)P value = 0.122P value > 0.05Binomialn=4,p=0.85566After that, the PDF for the output function expressed the aggregation of the two input variables that were randomly selected from the frequency and severity distributions. This step was repeated 1,000 times, and the obtained 1,000 output values represent the probability distribution of the output function (i.e., risk distribution).The risk distribution forHazard1 will be extensively explained to demonstrate the importance of utilizing this technique for SRA by DMs. Figure 3 illustrates the distribution of the risk of Hazard1, which is dependent on the distribution of the frequency and severity. Figure 4 shows a box plot for the risk values related to frequency and severity within the range of quartile percentiles (between 25% and 75%).Figure 3
Results of aggregation (the risk) of the frequency and severity forHazard1.Figure 4
Box plot for the inputs and outputs forHazard1.The frequency of this hazardous event was higher and more influential than the severity of the risk value. Therefore, if this hazard takes priority for treatment and is difficult to eliminate, at first glance, the focus should be on minimizing its frequency rather than its severity.Furthermore, the sensitivity analysis forHazard1 obtained from the tornado plot shown in Figure 5 affirms that the frequency variable has more influence on the mean of risks for Hazard1. Although the relative importance of both frequency and severity for this hazard is approximate, it is useful to estimate the low, base, and high outcomes for each variable and, therefore, should be focused on by the DMs.Figure 5
Sensitivity analysis forHazard1.In addition, it is important to ensure the reliability and steady state of the simulation results. We ran the simulation process eight times for the sameHazard1 to compare the results of cumulative plots for all runs, as presented in Figure 5. The means and cumulative probability plots for all runs relative to the risk values were found to be identical. Therefore, the simulation conditions with one run of simulations (1,000 samples) are suitable and sufficient to be run and analyze Hazard1. Figure 6 presents considerable useful information about the risk of Hazard1, which includes important statistics about this risk and its possible risk values. This outcome enables us to determine the ranking of this hazard relative to other hazards in the company by comparing them statistically to obtain a list of hazards ranked according to their priority for processing or elimination.Figure 6
Probability plot for the risk ofHazard1.Figure6 shows additional details for Hazard1 that are important for managers or DMs in safety administration, such as the confidence level for their decisions related to determining the processing priority. Regarding the remaining specified hazards, their frequencies and severities are aggregated to produce the values of risk. The risks for the specified hazards are summarized in Table 3, where the risk values for each hazard are presented as percentiles—25%, 50%, 75%, and 95%. In addition, the specific hazards are ranked depending on this summary.Table 3
Comparison of the statistics on the risk of hazards.
HazardParametersRanking for priority to be taken25%Mean75%95%SDFalling of materials due to overstacking and improper arrangements (H1)811.2014173.902Exposure to grinding flying particles (H2)46.69154.223Trip and fall while skimming due to poor housekeeping (H3)45.5792.284Forklift accidents due to speeding (H5)1013.4417214.791Table3 presents a list of the specified hazards that were analyzed, where the priority for taking action to eliminate or minimize them depends on the mean value of the risk. The hazardous event that must be prioritized for elimination or minimization from the studied hazards is Hazard5, followed by Hazard1. The priority sequence for the studied hazards is as follows: H5, H1, H2, and H3. Regarding Hazard4 (H4), because its frequency is not random, it is not suitable to calculate its distribution before analyzing a special cause for the variation. Therefore, it is excluded from this technique.


## 5.4.2. Using FST

The qualitative inputs gathered by the experts are presented in SupplementaryB.2 and are subjected to the mathematical equations and procedures presented in Section 3.2. To facilitate this process, a Matlab code was developed to read the experts’ judgments and perform calculations to obtain the aggregate of the frequencies and severities and then the fuzzy risks for the hazards. The results obtained for the code are presented in Table 4 and Figure 7, where the mean of the risk for the hazards is critical to rank the priority sequence for eliminating or controlling hazards as H5, H1, H4, H2, and H3. Therefore, mitigation plans related to these hazards should be implemented in this order.Table 4
Results of aggregation of experts’ evaluations using the fuzzy method.
Aggregate fuzzyHazard1Hazard2Hazard3Hazard4Hazard5MinMidMaxMinMidMaxMinMidMaxMinMidMaxMinMidMaxFrequency2.53.54.12.53.54.1344.82342.53.54.5Severity344.42341.52.53.5344.43.54.55.3Risks7.51418510.516.44.51016.8612188.81623.9Figure 7
Results for the risk of hazards calculated using the fuzzy method.Although this technique has been introduced to analyze hazards that have no historical data and require expert judgment, in this study, it was used for hazards that were already analyzed using MCS. Thus, the outcomes became more confident and FST was used to address the shortage with MCS regarding the data that were nonrandom, as in the case ofHazard4.


## 5.4.3. Comparative Analysis of the Proposed Models

Comparing the priority hazard lists using MCS and FST, both techniques obtained almost similar sequences, where the results of the priority of hazards to be treated by MCS wereH5-H1-H2-H3 (with H4 removed from MCS) and that to be treated by FST were H5-H1-H4-H2-H3. Moreover, the risk values of the hazards obtained using the two techniques were similar. The results showed that the integration of the two techniques for analyzing the different hazards occurs by overcoming the limitation of MCS regarding nonrandom data, as in the case of Hazard4 or when the hazards did not possess the recorded data. However, MCS can better analyze the hazards that have recorded randomized data and when determining expert opinions for the assessment is difficult and costly.


## 6. Conclusion

This study aims to develop a model for assessing and presenting risks considering uncertainty in data and information. The applicability and workability of this model were tested in the manufacturing environment of Riyadh, Saudi Arabia. The developed model, based on the integration of MCS and FST, accounted for processing the data uncertainty. Therefore, the proposed model was applied using two approaches: MCS and FST. A list of 371 hazards was reported in the hazard records in the production section of the company. However, it is difficult to analyze and apply the proposed model to all hazards, as the purpose here is application of the model. Thus, the factory’s safety administration was discussed and it was determined that only five hazards were to be studied using MCS and FST for performing the analysis and comparison. Using theModelRisk software, MCS was conducted considering that the hazards had 36-month historical data related to their frequency and severity, which were tested to ensure their random distribution. Then, the risk distributions, statistics with confidence levels for the risk values, and priority of the hazards were obtained. FST was used for the same five hazards, assuming that they did not have historical data and required expert judgment. Five experts in the company’s safety administration performed qualitative evaluation processes separately for the frequency and severity of the five hazards, which were transformed into fuzzy risk members using FST. A Matlab code was developed to solve the mathematical calculus for FST to obtain the minimum, medium, and maximum values for the risk of each hazard, and then, the prioritized list of hazards for elimination or treatment was obtained. The implementation of the proposed model demonstrates that integrating MCS and FST for risk assessment is effective and flexible, irrespective of whether the industry has historical data or missing or unknown information regarding the hazards. The risk assessment results obtained using the two techniques show the distributions, statistics, percentiles, and limits of risks for each selected hazard and enable more effective and flexible SRA decision-making. In addition, the results of comparisons of the priority lists for addressing hazards determined using both tools appear approximately the same, where the priority for the MCS-determined list is H5-H1-H2-H3 (with H4 removed from MCS) and that for the FST-determined list is H5-H1-H4-H2-H3. Therefore, mitigation plans related to these hazards should be implemented in these orders. According to these results, the proposed model seems to be reliable and applicable for assessing safety risks under uncertainty in the manufacturing industry. The main contributions of this study are summarized as follows:(1)
A model for SRA is developed for the manufacturing industry by considering the uncertainty and unavailability of data on the severity and likelihood of all hazards in the work environment.(2)
The proposed model considers a new methodology to obtain a clear picture for analyzing the risk with uncertainty through the intergradation of fuzzy sets and probabilistic risk analysis in occupational SRA.(3)
The proposed model helps DMs to obtain more information on the risk and probability of occurrence of hazards with certain risk values, which enable them to make the best decisions on the priority to eliminate or control that hazard [39].In the future, by using the proposed model with bow-tie analysis to determine the values of the root causes of hazards and their contribution to them by calculating the probability and severity using expert judgment to obtain the actual risk, a more accurate risk assessment can be achieved.



---

*Source: 6691124-2021-03-19.xml*