In this paper, the forced outage data from the Brazilian National Interconnected System (SIN) is provided and the statistics from the thermal power plants are computed. The SIN is characterized by a marked seasonality in electricity supply. In addition, th e expansion pattern of the Brazilian electric sector shows signs of exhaustion, and the demand for flexible thermal power plants, based on availability, requires outage management as a reliability -centered maintenance (RCM) strategy. In this work, the non -parametric Kruskal -Wallis test and Dunn’s pairwise -comparisons were chosen for evaluating the mean forced outage duration (MFOD) and the forced outage factor (FOF) using the data from the national electricity system operator (ONS) with R Software. The dist ribution fitting was provided using Weibull++ software from Reliasoft. Based on the MFOD and unit failure rate data, the FOF for Brazilian thermal power plants is 3.33% with a 90% probability and 95% confidence level. Finally, Brazilian thermal power plant s were benchmarked against North American power plants. Keywords: Kruskal -Wallis Test, Forced Outage, Thermal Power Plant 1. Introduction 1.1 The Brazilian power grid and its thermal power plants The Brazilian power grid , known as the national interconnected system (SIN ), currently has an installed capacity of 162 GW , of which hydroelectric power plants account for 101 .9 GW ( 62.9%), thermal and nuclear plants account for 22.9 GW (14.1%), and smaller -scale hydropower plant, biomass, wind , and sol ar plants account for the remaining 37.3 GW (23.0%), as shown in Figure 1 (Nacional do Sistema, 2020a). The recent revision of the long-term load forecasting of the SIN due the socioeconomic impacts caused by the new corona virus, indicates an average powe r rating of 65.774 GW in 2020 with a progressive increase to 76.612 GW by 2024 (Nacional do Sistema, 2020a). Figure 2 shows the long -term load forecasting. Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 99 Figure 1: The Brazilian power grid installed capacity by 2019 (Nacional do Sistema, 2020a). Figure 2 : The Brazilian power grid long -term load forecasting (Nacional do Sistema, 2020a). The Brazilian power grid is characterized by a marked seasonality in electricity supply due to the large production capacity during the rainy season and the run-of-river (ROR) operation of the Madeira complex and Belo Monte hydroelectric plants, thus, there is a constant recovery of the water reservoirs and displacement of the depletio n period. However, during the dry season , electricity production is shifted to thermal power plants (Nacional do Sistema, 2020a). According to the national electricity system operator (ONS), the SIN would have gone through a second critical period (90 mon ths) from the summer of 2011/2012 to the end of 2019 , as that of the historic critical period from June 1948 to November 1955 (Nacional do Sistema, 2020a). Due to the predominance of hydraulic -powered electricity generation in the SIN, the contracting of thermal power plants is based on availability. This is because the cost merit , economic benefit , and security risk mitigation needed to supply the national power demand need to the evaluated before considering the use of contracted energy. These factors are evaluated based on the remuneration for the availability of the thermal generator and the reimbursement of operating costs incurred with the actual plant dispatch (Center for Regulatory and Infrastructure Studies , 2017). Figure 3 presents the load in the SIN from 2016 to 2020. Hydro 63%Thermal 14%Wind 9%Biomass 8%Smaller -scale hydropower 4% Solar 2% 162 GW 2019 65,77468,63171,30273,92076,612 2020 2021 2022 2023 2024The Brazilian power grid long -term load forecasting (GW avg) Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 100 Figure 3 : Electrical load in SIN from 2016 to 2020 (MW avg). Adap ted from (Nacional do Sistema, 2017, 2018, 2019a, 2020b). Currently, thermal power generation in Brazil has a high variable unit cost (CVU), with 40% of the installed capacity having a CVU of over 250 BRL/MWh (~47 USD/MWh ). This indicates that thermal pow er plants are used for electrical dispatch only in extremely unfavorable hydrological conditions, thus, increasing the challenges for maintaining the water level in reservoirs (Nacional do Sistema, 2020a). The ONS, a private entity created in 1998, coordi nates and controls the operations of the power generation and transmission facilities integrated in the SIN . In addition, the ONS evaluates the future and short -term conditions of these facilities (Sousa, 2009). The decision -making process of the ONS is regulated by the grid procedures and approved by the Brazilian Electricity Regulatory Agency (ANEEL) through public hearings. In addition, the ONS provides ANEEL with input data and optimization models , to enable t he reproduction of the method used for determining the plants' dispatch (Sousa, 2009). In summary , the ONS functions to minimize the service cost, thus, decreasing the risk of insufficient power supply in each subsystem that composes the SIN. Thermal powe r stations operation provides efficiency gains during the dry and low -pour period, thus, increasing the firm energy of the system (Sousa, 2009). According to the National Confederation of Industry, the national thermoelectric system of Brazil was original ly contracted as a reserve for sporadic electrical performance during unfavorable hydrological periods . However, the data presented by Chamber of Electric Energy Commercialization (CCEE) show that the percentage of thermal power plants in the SIN load is approximately 16.5% (De Comercialização de Energia Elétrica, 2020) at most times in the base operation (National Confederation of Industry, 2018). The base operation is aimed at complementing electrical generation at all load levels. In addition to this, t here is also the peak operation, which occurs in periods of fluctuating demand, and the operation aimed at performing ancillary services such as (Ribeiro, 2019): controlling (i) primary frequency, (ii) secondary frequency, (iii) reactive support, (iv) self -restoration, (v) special protection system (SEP ), and (vi) complementary dispatch for maintaining the operational power reserve (Nacional do Sistema , 2019b). However, the expansion pattern of the Brazilian electric sector shows signs of exhaustion, with the inclusion of the ROR and intermittent renewable sources (wind and solar). Consequently, there is a growing demand for flexible thermal power plants operated in an electric dispatch mode (National Confederation of Industry, 2018). Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 101 The structural trend in the SIN involves the gradual reduction in the capacity of country's reservoirs regularization. Therefore , to recover water storage, thermal power complementation needs to be carried out more frequently and for longer periods , even during hydrological periods close to the long -term average (National Confederation of Industry, 2018). In a study carried out by (Nogueira et al., 2019), they reported that 50% of Brazilian hydroelectric plants are over 20 years old and 32% are over 40 yea rs old. The aging of these plants is the major cause for the increasing frequency and duration of maintenance interventions, whether as scheduled or due to forced shutdowns. ANEEL calculates the unavailability of hydroelectric plants based on two key perf ormance indicators (KPI): scheduled unavailability rate (IP) and equivalent rate of forced unavailability (TEIF). The IP quantifies the percentage of hours that the turbo -generator was turned off to carry out scheduled maintenance , while the TEIF quantifie s the percentage of hours that the equipment remained in forced shutdown (Ministério de Minas e Energias, 2014 ). The reference values of the KPI are shown in Table 1. Table 1: Reference values of the forced outages of hydroelectric plants (Ministério de M inas e Energias, 2014 ). Hydroelectric MW Trb/Gen TEIF (%)¹ 001–029 2.068 030–059 1.982 060–199 1.638 200–699 2.133 700–1300 3.115 ¹ TEIF is harmonized with Forced Outage Factor (FOF) indicator from IEEE Std 762TM-2006 (IEEE, 2007) The Brazilian electrical system (SEB) contracts are based on physical guarantee through energy guaranty contracts certified by the Ministry of Mines and Energy (MME) (Castro et al., 2016). One of the main performance risks for power plants is related to the decreasing physical guarantee due to unavailability, thus, leading to a reduction in the revenue of the commercialization of Electric Energy in Regulated Environments (CCEAR) and an increase in the additional costs for the acquisition of ballast in the S hort-Term Market (MCP) (Martins, 2013). The unavailability of each plant is computed by the ONS using the IP and TEIF . These values are calculated monthly, based on the moving average of the previous 60 months. If the verified unavailability KPIs are high er than those declared, there will be a reduction in the physical guarantee (Martins, 2013). Therefore, to minimize the risk of not supplying the national power demand, the combination of these variables is important : (i) basis operation, (ii) aging of hydroelectric plants, and (iii) flexible dispatch requirements : and demands from thermal power plants , high availability , and reliability standards. 1.2 Maintenance at Thermal Power Plants Maintenance Management at the generating units is important for an economically optimized electrical dispatch in hydrothermal generation systems. However, it is challenging to choose the best schedule for preventive maintenance to minimize the operating co st of the generating agency , maximize the reliability of the system , and extend the units' operational life. In addition, with an increase in the size of the generating unit, the challenges also increase (Balaji et al., 2016). The maintenance costs of ther mal power plants are greatly influenced by: (i) the type of starting, (ii) the frequency of starting, and (iii) the loading pattern (Dipak, 2015a). Figures 4 and 5 show examples of the effect of these factors on a gas turbine. Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 102 123456 1/1 1/10 1/1000Maintenance cost factor Starts per fired hours Fast Start Normal StartEffect of starting on maintenance cost for Thermal Power Plants 1/100 Figure 4: Effect of starting on maintenance cost. Adapted from (Dipak, 2015a). Figure 5: Effect of loading on maintenance cost. Adapted from (Dipak, 2015a). Currently, the steam thermal power plants have larger sizes with a high steam -generating capacit y and highly sophisticated firing system. However, with an increase in the unit size and capacity of thermal plants , forced outages have greater significance not only on the larger loss of revenue , but also greater risk of injury and damage to the plants (Dipak, 2015b). To effectively manage these challenges, the reliability -centered Maintenance (RCM) process is carried out. The RCM process is a common -sense procedure for creating maintenance strategies to preserve assets’ functions. The standard SAE JA101 1, published in 1999, sets the criteria that any process must comply with to be considered an RCM. It establishes that for a process to be acknowledged as RCM, it must follow these seven steps (Sifonte and Reyes -Picknell, 2017): i. Delineate the operational c ontext and the functions , and associated desired standards of assets performance (operational context and functions) ii. Determine how an asset can fail to fulfill its functions (functional failures) iii. Define the causes of each functional failure (failure modes) Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 103 iv. Describe what happens when each failure occurs (failure effects) v. Classify the consequences of failure (failure consequences) vi. Determine what should be performed to predict or prevent each failure (tasks and intervals) vii. Decide whether other failure managemen t strategies may be more effective (one -time changes). Some uncertainties in thermal power plant dispatch can cause deviations from the system balance, which sometimes require inefficient and costly last -minute solutions in the near real -time time frame (Makarov et al., 2017). These uncertainties include : (i) uninstructed deviations of conventional g enerators from their dispatch set points, (ii) generator forced outages, (iii) generator failures to start up, (iv) load drops, (v) losses of major transmission lines, and (vi) frequency variation (Makarov et al., 2017). Outage management organization and administration ensures the safe and effective implementation and control of the maintenance activities during planned and forced outages. Outage planning and performance takes into consideration the safety, quality , and schedule in this order. Thus, the maintenance planning and scheduling process should reflect this (Lipár, 2012). The objective of this study is to perform a multiple pairwise -comparisons of the annual forced shutdown statistics by ONS using the Kruskal -Wallis non -parametric test to verify if the populations’ distributions have changed over the years. This allows for the computing of the reliability analysis for the forced outage data and to analyze the availability of thermal power plants an d comparing it to benchmarks. 2. Materials and Methods 2.1 Thermal power plants: forced shutdown statistics Forced shutdown involves the unscheduled removal of a component out of service due to failure or emergency shutdown. Forced shutdown requires the manual or automatic switching off of an equipment to avoid risks to the physical integrity of the people or the environment, or damage to the equipment, or other consequences to the electrical sy stem. It includes accidental shutdown (without disturbance in the SIN ), incorrect shutdown (with disturbance in the SIN ), and shutdown resulting from SEP or configuration actions (Nacional do Sistema, 2019 c). The ONS provides performance and statistical reports issued from the generating units belonging to the SIN. Therefore, to apply the methodology used in this work, the report ONS DPL -REL -0099/2020: Relatório de Análise Estatística de Desligamentos Forçados de Equipamentos Referente ao ano de 2019 (Statistical Analysis Report on Forced Equipment Shutdowns for 2019) was analyzed (Nacional do Sistema, 2019 c). The ONS reported that the occurrence of the unsatisfactory operational performance of a thermal power plant directly impacts the security and reliability of electricity supply and energy tariffs, given that its unavailability causes the generation of another power unit with a higher CVU (Agência Nacional de Energia Elétrica, 2020). The consolidated data of the forced shutdowns of thermal generators are presented using the indicators : Mean Forced Shutdown Duration for Transmission and Generation Functions (DMDFF ) and Forced Shutdown Rate for Transmission and Generation Functions (TDFF) (Nacional do Sistema, 2019 c). The DMDFF indicator aims to manage the performance of the transmission, transformation, reactive control , and the generation functions, regarding the mean duration of the forced shutdowns during the period considered. The DMDFF is calculated as follow s (Nacional do Sistema, 2019d): 𝐷𝑀𝐷𝐹𝐹 = 𝐹𝑜𝑟𝑐𝑒𝑑 𝑜𝑢𝑡𝑎𝑔𝑒 ℎ𝑜𝑢𝑟𝑠 (𝐹𝑂𝐻 ) 𝑁𝑢𝑚𝑏𝑒𝑟𝑠 𝑜𝑓 𝑢𝑛𝑝𝑙𝑎𝑛𝑛𝑒𝑑 𝑜𝑢𝑡𝑎𝑔𝑒𝑠 𝑓𝑟𝑜𝑚 𝑖𝑛 𝑠𝑒𝑟𝑣𝑖𝑐𝑒 𝑠𝑡𝑎𝑡𝑒 Eq. 1 Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 104 The TDFF indicator aims to manage the rate of forced shutdowns on the transmission and generation functions during service hours at the period considered. The TDFF is calculated is as follows (Nacional do Sistema, 2019d): 𝑇𝐷𝐹𝐹 = 𝑁𝑢𝑚𝑏𝑒𝑟𝑠 𝑜𝑓 𝑢𝑛𝑝𝑙𝑎𝑛𝑛𝑒𝑑 𝑜𝑢𝑡𝑎𝑔𝑒𝑠 𝑓𝑟𝑜𝑚 𝑖𝑛 𝑠𝑒𝑟𝑣𝑖𝑐𝑒 𝑠𝑡𝑎𝑡𝑒 𝑆𝑒𝑟𝑣𝑖𝑐𝑒 𝐻𝑜𝑢𝑟𝑠 (𝑆𝐻)  8760 Eq. 2 Where: The constant value 8760 is the annualization factor – 24 h for 365 d. The TEIF, DMDFF, and TDFF indicators may be harmonized with the IEEE Standard Definitions for use in Reporting Electric Generating Unit Reliability, Availability, and Productivity , IEEE Std 762TM-2006 as follows (IEEE, 2007): 𝑀𝐹𝑂𝐷 = 𝐷𝑀𝐷𝐹𝐹 Eq. 3 𝑈𝑛𝑖𝑡 𝐹𝑎𝑖𝑙𝑢𝑟𝑒 𝑅𝑎𝑡𝑒 = 1 𝑀𝑆𝑇𝐹𝑂=𝑇𝐷𝐹𝐹 Eq. 4 𝐹𝑂𝐹 = (𝐹𝑂𝐻 𝑃𝐻)𝑥100 =𝑇𝐸𝐼𝐹 =𝑈𝑛𝑖𝑡 𝐹𝑎𝑖𝑙𝑢𝑟𝑒 𝑅𝑎𝑡𝑒 𝑥 𝑀𝐹𝑂𝐷 𝑃𝐻𝑥100 Eq. 5 Where MFOD is the mean forced outage duration and MSTFO is the mean service time to forced . PH is the period hours or active hours and it represents the number of hours a unit was in the active state, FOF is the forced outage factor, i.e., the fraction of a given operating period in which a generating unit was not available due to forced outage s. The IEEE nomenclature was adopted to present the consolidated data in this work. 2.2 Kruskal -Wallis test and Multiple pairwise -comparisons The data from the statistic report provided by ONS presents the information grouping over a 5 years interval (2015 to 2019). In this study, the Kruskal -Wallis (KW) test was performed to test the null hypothesis ( H0) that the annual probability distributions of the Brazilian thermal generators are equal . The KW test is a well-known non-parametric test applied in a wide range of disciplines such as engineering and manufacturing applications, medicine, biology, psychology , and education (Ostertagová et al., 2014). The KW test , named in honor of the American statisticians William Krus kal and W. Allen Wallis , was created in 1952 and is a non -parametric test used to compare three or more populations. The test does not make assumptions about normality. However, it assumes that the populations have the same distribution , and that the sampl es are random and independent (Ostertagová et al., 2014; Ronald et al., 2007; Moreno et al., 2019). The KW test requires that the data come from continuous probability distributions. The H0 being tested by the KW statistic assumes that all the population distributions are equal, and the alternative hypothesis ( H1) assumes that at least two population distributions differ from each other (Ostertagová et al., 2014; Ronald et al., 2007). Data are pooled across groups and ranked from the lowest value of the de pendent variable to the highest value. In case of a tie, the average rank is attributed to the tied experimental values (Andrew, 1998). According to McDonald, , the KW test also assumes that the variation within the groups is equal (homoscedasticity), how ever, groups with different standard deviations have different distributions. Thus , if different groups have different shapes, the KW test may give inaccurate results (McDonald, 2007; Fagerland and Sandvik, 2009). Since the KW test assumes that different groups have similar distribution, and groups with different standard deviations have different distributions , if the data are heteroscedastic, KW is no better than the one-way ANOVA and may be worse. For heteroscedastic data , the Welch's ANOVA is preferred (McDonald, 2007). The H0 of KW assumes that the distribution for all k populations are similar, while the H1 assumes that the distribution of at least two population differs. H0 and H1 can be expressed as follows (Andrew, 19 98): 𝐻0:𝐹1(𝑥)=𝐹2(𝑥)=⋯=𝐹𝑘(𝑥) Eq. 6 Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 105 𝐻1:∃1≤𝑖,𝑙≤𝑘: 𝐹𝑖(𝑥)≠𝐹𝑙(𝑥) Eq. 7 The sum of the ranks R i is calculated for each group i (i = 1, 2, . . . , k) of size ni, then the test statistic H is calculated, which represents the variance of the ranks among all groups, with an adjustment for the number of ties (Hecke, 2010). 𝐻=(12 𝑁(𝑁+1)∑𝑅𝑖2 𝑛𝑖𝑘 𝑖=1) −3(𝑁+1), 𝑁=∑𝑛𝑖𝑘 𝑖=1 Eq. 8 Where ni is the sample size for the ith group, Ri is the rank -sum for the ith group, N is the total sample size , and k is the number of groups. If there are many ties in the samples , a correction factor must be applied, and the test statistic corrected for the ties will be (Ostertagová et al., 2014; Hollander et al., 2013): 𝐻∗=𝐻 𝑓∗, 𝑓∗=1−(∑ (𝑡𝑖3−𝑡𝑖)𝑚 𝑖=1 𝑁3−𝑁) Eq. 9 Where ti is the number of ties in the ith group of the m ties groups, f *is the correction factor, H is the KW statistic , and H* is the KW statistic corrected for ties. Whenever H0 is true and either (Ostertagová et al., 2014): {𝑘=3,𝑛𝑖≥6 𝑓𝑜𝑟 𝑖=1,2,3 𝑘>3,𝑛𝑖≥5 𝑓𝑜𝑟 𝑖=1,2,…,𝑘 Eq. 10 The distribution of the test statistic, H, is approximated using the chi -square distribution with ( k-1) degrees of freedom. The H0 is rejected on the right -hand tail of the chi -square distribution (Ostertagová et al., 2014) After the test statistic is calculated, the p -value is then calculated . The p -value is defined as the probability of observing the given value of the test statist ic, or greater, under the H0 (Ferreira and Patino, 2015). That is, the H0 is rejected on a significance level α, when a one-sided p -value is less than the significance level (Ferreira and Patino, 2015): 𝑃[𝜒2(𝑘)>𝐻]< 𝛼 Eq. 11 Where α is the significance level , i.e., the probability of making the wrong decision when the H0 is true, χ² is the quantile of the chi -square distribution , and H is the statistic from the KW test. 2.3 Effect size for the Kruskal –Wallis test The statistics of the effect size for the KW test provide the degree to which the data of one group has higher ranks than that of another group. The effect size is related to the probability that the value from one group will be greater than the value from an other group (Mangiafico, 2016 ). The eta -squared coefficient can be calculated as the measure of the KW test effect size. According to Prajapati et al., eta is a measure of association and is the proportion of the total variance that is attributed to an effect (Prajapati et al., 2010 ). Eta -squared ranges from 0 to 1 , and as a rule , 0.01 is a small effect, 0.06 is a moderate effect , and 0.14 is a large effect (Laerd Statistics, 2020). 𝐸𝑅2=𝐻 (𝑁2−1)/(𝑁+1), {0.01 𝑡𝑜 0.06,𝑠𝑚𝑎𝑙𝑙 𝑒𝑒𝑓𝑓𝑒𝑐𝑡 0.06 𝑡𝑜 0.14,𝑚𝑜𝑑𝑒𝑟𝑎𝑡𝑒 𝑒𝑓𝑓𝑒𝑐𝑡 𝑓𝑟𝑜𝑚 0.14 𝑜𝑛,𝑙𝑎𝑟𝑔𝑒 𝑒𝑓𝑓𝑒𝑐𝑡 Eq. 12 Where H is the statistic from the KW Test, k is the number of groups, and N is the total number of observations. According to Ferreira and Patino, there is a misconception that a very small p -value indicates a highly relevant difference between groups . However, it is necessary to consider the effect size, as it may indicate that the sample size should be increased. The authors recommend preferably reporting the mean values for each group, the difference, and the 95% confidence interval , and then the p -value (Ferreira and Patino, 2015) . Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 106 This significa nt result in a KW test indicates that there are group differences, however it does not indicate which groups. Thus, a post hoc procedure can be used to determine which groups are different from each other (Andrew, 1998). According to Laerd Statistics, if the populations distributions have similar shapes, then, the medians can be compared to evaluate the distribution differences ( Laerd Statistics, 2020) . However, when the distribution shapes are different , the mean rank from the KW test shou ld be considered. This is because with an increase in the group’s mean rank, the observation values in that group increases in comparison to those of the other groups (Minitab 19 Support, 2020). The concept is shown in Figure 6. Figure 6 : Evaluation cri teria for the populations shape in KW test. Adapted from (Laerd Statistics, 2020). The KW test does not assume normality but assumes that the shapes of the distributions in different groups are similar . This suggests that non -parametric tests are not a go od solution for heteroscedastic data (McDonald, 2007). According to Dinno, in the case of a rejected H0, the Dunn’s test should be performed after the KW test (Dinno, 2015) . The Dunn’s test is based on the normal distribution with Bonferroni correction, where nc is the possible number of two -to-two comparisons that will be made between the k groups (Moreno et al., 2019). 𝑛𝑐=𝐶2𝑘=𝑘! 2!(𝑘−2)!=𝑘(𝑘−1) 2 Eq. 13 The test involves the comparison of the module of the differences between the means of the ranks for two groups |𝑅̅𝑖−𝑅̅𝑗| with the least significant difference (LSD) (Moreno et al., 2019). 𝐿𝑆𝐷 =𝑍(𝛼𝑐)√[𝑁(𝑁+1) 12(1 𝑛𝑖+1 𝑛𝑗)], 𝛼𝑐=𝛼 𝑘(𝑘−1)/2 Eq. 14 Where: LSD is the least significant difference , αc is the Bonferroni's correction, α is the significance level, ni and nj are the sample sizes of two compared populations, N is the total sample size, and Zαc comes from negative Z score table for αc (Moreno et al., 2019). When |𝑅̅𝑖−𝑅̅𝑗| ≥𝐿𝑆𝐷, H0 is rejected, and the pairwise comparison shows that two compared populations are significantly different (Moreno et al., 2019). 2.4 Reliability analysis for the for ced outage data The analysis method proposed in this work considers all the thermal generators dispatched centrally and belonging to the selected facilities according to the ONS's criteria (see Table 2) (Nacional do Sistema, 2019c): i. Thermal power plants with an effective power equal to or greater than 300 MW. ii. Thermal power plants with 200 MW ≤ effective power <300 MW, with a transformation equal to or greater than 230 kV. iii. Thermal power plants powered by natural gas, coal, and nuclear sources. Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 107 Table 2 : Centrally dispatched thermal power plants according to the ONS’s criteria for statistical analysis (Nacional do Sistema, 2019c). Number of centrally dispatched Thermal Power Plants 2015 2016 2017 2018 2019 106 108 109 111 112 The MFOD and unit failure rate indicators from the thermal power plants are presented for 2015 to 2019, as shown in Table 3 and Appendix 1. While the causes of the forced outages are shown in Figure 7 (Nacional do Sistema, 2019c): Table 3 : Centrally dispatched Thermal Power Plants KPIs according to ONS’s criteria (Nacional do Sistema, 2019c). Centrally dispatched Thermal Power Plants KPIs KPI/Year 2015 2016 2017 2018 2019 Mean Forced Outage Duration 7.103 7.346 5.383 6.312 11.631 Unit Failure Rate 8.852 4.761 5.452 3.976 5.112 Figure 7 : Causes of forced outages for the thermal power plants. Adapted from (Nacional do Sistema, 2019c). Based on the criteria by ONS , there are four major reasons for the forced power outage of power plants in SIN (Nacional do Sistema, 2019c): i. Internal outage , which is related to the main parts of the power plant such as the insulators, the primary winding from transformers, stator, bearing and generator shaft, circuit breakers , and others. ii. Seco ndary outage , which is related to the complementary or auxiliary equipment of the power plants such as wiring, protection, control, command, auxiliary services, ventilation , and cooling system. In addition, it also includes outages due to the incorrect per formance from the protection system in case of external failures. iii. External outage , which is related to the outage due to the correct performance of the protection system due to failures (acting as a back -up protection) or due to overload caused by outage f rom a third party. In addition, it also includes the outage caused by system configuration. iv. Operational outage , which is related to systemic electrical conditions such as oscillations, overvoltage, over-frequency, overload , and other systemic causes. 3. Results The KW non -parametric test was computed using R software. The R software is a free software used for statistical computing and for the construction of graphics that can be downloaded and distributed free of charge under the GNU license (Landeiro, 2011; R Development Core Team, 2011 ). The prerequisites for computing the KW test using R software require the R packages: (i) tidverse for da ta manipulation and visualization (Wickham et al., 2019), (ii) ggpubr for creating ready to publish plots (Alboukadel, 2020a), and (iii) rstatix , which provides R functions for statistical analyses (Alboukadel, 2020b). 931 813 5710 51508 467 24 1 16586527 45 2 12422 384 18 15 5567504 512 10 Total Secondary External Operational InternalCauses of forced outages for Thermal Power Plants 2015 2016 2017 2018 2019 Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 108 The available data from the ONS were prepared using comma -separated values (CSV) grouped into two columns: weight (KPI value) and group (year). library(tidyverse) library(ggpubr) library(rstatix) mfod <- read.csv2("dmdff.csv") rate <- read.csv2("tdff.csv") Subsequently, the summar y statistics was computed by groups using the kruskal_test function from the rstatix package (Data Novia, 2020 ). Since the number of groups used in this study was greater than 3 and the sample size was greater than 5 for all groups , according to Eq. 10, the distribution of the test statistic , H, was well approximated using the chi -square distribution with (5 - 1) degrees of freedom. The significance level chosen was α = 0.05. Manual calculations were pe rformed to determine the mean ranks ( Ri/ni), and the results are shown in Tables 4 and 5. Table 4: KW test mean rank for the mean forced outage duration . Year Sample size Rank -sum (Rank -sum)²/sample size Mean rank 2015 94 20 699.00 4 557 963.84 220.20 2016 80 16 304.50 3 322 959.00 203.81 2017 86 16 222.50 3 060 110.54 188.63 2018 86 18 561.00 4 005 938.62 215.83 2019 86 21 741.00 5 496 175.36 252.80 Total 432 93 524.00 20 443 147.36 - Table 5: KW test mean rank for the unit failure rate. Year Sample size Rank -sum (Rank -sum)²/sample size Mean rank 2015 97 28 123.00 8 153 640.51 289.93 2016 82 17 365.00 3 677 356.40 211.77 2017 89 20 415.50 4 683 063.37 229.39 2018 88 16 236.00 2 995 542.00 184.50 2019 93 18 885.50 3 835 076.45 203.07 Total 449 101 025.00 23 344 678.74 - Mean forced outage duration statistics : mfod %>% group_by(group) %>% get_summary_stats(weight, type = "common") # A tibble: 5 x 11 group variable n min max median iqr mean sd se ci <int> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 2015 weight 94 0.283 51.8 5.44 5.12 6.91 6.33 0.653 1.30 2 2016 weight 80 0.175 64.7 5.41 5.16 7.14 8.50 0.95 1.89 3 2017 weight 86 0.508 28.2 4.38 4.72 5.74 4.89 0.527 1.05 4 2018 weight 86 0.267 20.2 5.49 5.30 6.29 4.12 0.444 0.883 5 2019 weight 86 0.467 72.0 6.36 7.85 12.1 15.1 1.63 3.24 Unit failure rate statistics : rate %>% group_by(group) %>% get_summary_stats(weight, type = "common") Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 109 # A tibble: 5 x 11 group variable n min m ax median iqr mean sd se ci <int> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 2015 weight 97 1.00 27.4 8.02 9.01 9.10 5.87 0.596 1.18 2 2016 weight 82 1 18.2 5.00 4.02 5.90 4.27 0.472 0.938 3 2017 weight 89 1 20.2 5.02 4.02 6.18 3.79 0.402 0.799 4 2018 weight 88 1 80.3 4.51 4.04 5.64 8.78 0.936 1.86 5 2019 weight 93 1 38.8 4.07 4.04 5.94 5.90 0.612 1.22 Where N is the number of individuals, min is the minimum value, max is maximum value, median is the median, mean is the mean, sd is the standard deviation of the mean, se is the standard error of the mean , and ci is the 95% confidence interval of the mean (Kassambara, 2020a). Mean forced outage duration Kruskal -Wallis test: res.kruskal < - mfod %>% kruskal_test(weight ~ group) res.kruskal # A tibble: 1 x 6 .y. n statistic df p method * <chr> <int> <dbl> <int> <dbl> <chr> 1 weight 432 12.5 4 0.0142 Kruskal -Wallis Unit failure rate Kruskal -Wallis test: res.kruskal < - rate %>% kruskal_test(weight ~ group) res.kruskal # A tibble: 1 x 6 .y. n statistic df p method * <chr> <int> <dbl> <int> <dbl> <chr> 1 weight 449 36.5 4 0.000000231 Kruskal -Wallis Where .y. is the y variable used in the test, n is the sample count, statistic is the Kruskal -Wallis rank -sum statistic used to compute the p -value, df is the degree of freedom, p is the computed p -value , and method is the statistical test used to compare groups (Kassambara, 2020b). The p-values of the computed KW test was less than a 0.05 significance level , thus, rejecting the H0 for both KPIs, indicating that at least two population distributions differ from each other. After computing the KW p -values and plotting the boxplots, the statistics of effect size for the test was verified by computing the eta -squared from the kruskal_effsize function: Mean forced outage duration effect size: mfod %>% kruskal_effsize(weight ~ group) # A tibble: 1 x 5 .y. n effsize method magnitude * <chr> <int> <dbl> <chr> <ord> 1 weight 432 0.0198 eta2[H] small Unit failure rate effect size: rate %>% kruskal_effsize(weight ~ group) # A tibble: 1 x 5 .y. n effsize method magnitude Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 110 * <chr> <int> <dbl> <chr> <ord> 1 weight 449 0.0732 eta2[H] moderate Where .y. is the y variable used in the test, n is the sample counts, effsize is the estimate of the effect size, method is the eta -squared , and magnitude is the magnitude of effect size ( Kassambara, 2020c). This significant results in the KW tests indicate that there were group differences, however it does not indicate which groups. Thus, a Dunn test procedure was used to determine which groups were diffe rent from each other. Mean forced outage duration Dunn test: # Pairwise comparisons pwc <- mfod %>% dunn_test(weight ~ group, p.adjust.method = "bonferroni") pwc # A tibble: 10 x 9 .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr> 1 weight 2015 2016 94 80 -0.863 0.388 1 ns 2 weight 2015 2017 9 4 86 -1.69 0.0902 0.902 ns 3 weight 2015 2018 94 86 -0.235 0.814 1 ns 4 weight 2015 2019 94 86 1.75 0.0801 0.801 ns 5 weight 2016 2017 80 86 -0.782 0.434 1 ns 6 weight 2016 2018 80 86 0.620 0.535 1 ns 7 weight 2016 2019 80 86 2.53 0.0115 0.115 ns 8 weight 2017 2018 86 86 1.43 0.153 1 ns 9 weight 2017 2019 86 86 3.37 0.000751 0.00751 ** 10 weight 2018 2019 86 86 1.94 0.0521 0.521 ns Unit failure rate dunn test: # Pairwise comparisons pwc <- rate %>% dunn_test(weight ~ group, p.adjust.method = "bonferroni") pwc # A tibble: 10 x 9 .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr> 1 weight 2015 2016 97 82 -4.02 0.0000593 0.000593 *** 2 weight 2015 2017 97 89 -3.18 0.00148 0.0148 * 3 weight 2015 2018 97 88 -5.52 0.0000000340 0.000000340 **** 4 weight 2015 2019 97 93 -4.61 0.00000397 0.0000397 **** 5 weight 2016 2017 82 89 0.887 0.375 1 ns 6 weight 2016 2018 82 88 -1.37 0.171 1 ns 7 weight 2016 2019 82 93 -0.443 0.658 1 ns 8 weight 2017 2018 89 88 -2.30 0.0214 0.214 ns 9 weight 2017 2019 89 93 -1.37 0.171 1 ns 10 weight 2018 2019 88 93 0.962 0.336 1 ns Where .y. is the y variable used in the test, group1 and group2 are the compared groups in the pairwise tests, n1 and n2 are the sample counts, statistic is the test statistic (z -value) used to compute the p -value, p is the p -value, p.adj is the adjusted p -value , and p.adj.signif is the significance level of the p- adjusted p -values, respectively (Kassambara, 2020d). Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 111 The following convention for symbols indicates statistical significance (STHDA, 2017): ns means p > 0.05, * means p ≤ 0.05, ** means p ≤ 0.01, *** m eans p ≤ 0.001, and **** means p ≤ 0.0001. To identify the differences between the consolidated annual values, boxplot graphs were plotted, as shown in figures 8 and 9. Mean forced outage duration boxplot: # Visualization: box plots with p-values pwc <- pwc %>% add_xy_position(x = "group") ggboxplot(mfod, x = "group", y = "weight", color = "group", palette = c("#00AFBB", "#E7B800", "#FC4E07", "#3CB371", "#BA55D3"), order = c("2015", "2016", "2017", "2018", "2019"), ylab = "Mean Forced Outage Duration", xlab = "Year" ) + stat_pvalue_manual(pwc, hide.ns = TRUE) + labs( subtitle = get_test_label(res.kruskal, detailed = TRUE), caption = get_pwc_label(pwc) ) Unit failure rate boxplot: # Visualization: box plots with p -values pwc <- pwc %>% add_xy_position(x = "group") ggboxplot(rate, x = "group", y = "weight", color = "group", palette = c("#00AFBB", "#E7B800", "#FC4E07", "#3CB371", "#BA55D3"), order = c("2015", "2016", "2017", "2018", "2019"), ylab = "Unit Failure Rate", xlab = "Year" ) + stat_pvalue_manual(pwc, hide.ns = TRUE) + labs( subtitle = get_test_label(res.kruskal, detailed = TRUE), caption = get_pwc_label(pwc) ) Figure 8 : Mean forced outage duration boxplot. Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 112 Figure 9 : Unit failure rate boxplot. To verify the distribution fitting for the annual data from the MFOD and the unit failure rate, the data was fitted using the distribution wizard from software Weibull++ (version 19.0.2.1075, from Reliasoft ). According to Reliasoft , the distribution wizard performs multiple goodness of fit tests to determine the best distribution for a data set based on the chosen parameter estimation method. In addition, it performs three goodness of fit tests to determine the rank of the distributions (Reliasoft, 2017 ): i. The Kolmogorov -Smirnov test (Goodness of fit – GOF) tests used to determine the statistical difference (the difference between the expected and obtained results). Default weighted in 40% for the Maximum likelihood (MLE) analysis method. ii. The Correlation coefficient test (Plot fit – PLOT) measures how well the plotted points fit a straight line. Default weighted in 10% for the MLE analysis method. iii. The Likelihood value test (Likelihood Ratio – LKV) computes the value of the log -likelihood function, given the parameters of the distribution. Default weighted in 50% for the MLE analysis method. Tables 6 and 7 show the distribution fitting for the MFOD and u nit failure rate. The top ranked distributions implemented for each year and KPI. Table 6 : Mean forced outage duration distribution fitting. Mean forced outage duration distribution fitting 2015 2016 2017 2018 2019 Distribution/ Ranking Distribution/ Ranking Distribution/ Ranking Distribution/ Ranking Distribution/ Ranking Loglogistic 1 Loglogistic 1 G-Gamma 1 G-Gamma 1 Loglogistic 1 G-Gamma 2 G-Gamma 2 3P-Weibull 2 3P-Weibull 2 G-Gamma 2 Gamma 3 Lognormal 3 Gamma 3 2P-Weibull 3 Lognormal 3 Lognormal 3 2P-Exponential 4 Lognormal 4 Gamma 4 3P-Weibull 4 2P-Weibull 4 3P-Weibull 5 2P-Exponential 5 Loglogistic 5 2P-Weibull 5 3P-Weibull 5 Gamma 6 Loglogistic 6 Logistic 6 2P-Exponential 6 2P-Exponential 6 1P-Exponential 7 2P-Weibull 7 Lognormal 7 1P-Exponential 7 Logistic 7 2P-Weibull 8 Logistic 8 2P-Exponential 8 Gamma 7 1P-Exponential 8 Logistic 9 1P-Exponential 9 Normal 9 Logistic 8 Normal 9 Normal 10 Normal 10 Gumbel 10 Normal 9 Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 113 Gumbel 10 Gumbel 11 Gumbel 11 1P-Exponential 11 Gumbel 10 Table 7 : Unit failure rate distribution fitting. Unit failure rate distribution fitting 2015 2016 2017 2018 2019 Distribution/ Ranking Distribution/ Ranking Distribution/ Ranking Distribution/ Ranking Distribution/ Ranking 3P-Weibull 1 3P-Weibull 1 G-Gamma 1 Lognormal 1 G-Gamma 1 G-Gamma 2 Gamma 2 3P-Weibull 2 G-Gamma 2 Lognormal 1 Gamma 3 G-Gamma 3 Gamma 3 3P-Weibull 3 Loglogistic 2 2P-Weibull 4 2P-Exponential 4 2P-Weibull 4 Loglogistic 4 2P-Exponential 3 Lognormal 5 2P-Weibull 5 Loglogistic 5 Gamma 5 3P-Weibull 4 Loglogistic 6 Lognormal 6 Logistic 6 2P-Exponential 6 2P-Weibull 5 Logistic 7 Loglogistic 6 Lognormal 7 2P-Weibull 7 Gamma 6 2P-Exponential 8 Logistic 7 2P-Exponential 8 1P-Exponential 8 1P-Exponential 7 Normal 9 1P-Exponential 8 Normal 9 Logistic 9 Logistic 8 1P-Exponential 10 Normal 9 1P-Exponential 10 Normal 10 Normal 9 Gumbel 10 Gumbel 10 Gumbel 10 Gumbel 11 Gumbel 10 The distribution shapes were plotted to evaluate the scale and skewness within the groups, as an indication of heteroscedasticity, as shown in Figures 10 and 11. Figure 10 : Probability density functions for Mean Forced Outage Duration data. Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 114 Figure 11 : Probability density functions for Unit Failure Rate data. Then, the distribution locations were computed considering the B50% life, assuming the groups had equal variation (homoscedasticity). The median is the value that the variable has a 50% probability of exceeding (Camarillo et al., 2017). Tables 8 and 9 show the B50% life as per computed using the Quick Calculation Pad from Weibull++ Software. The Two-sided Confidence Level of 95% was considered, which is the measure of the imprecision of the true effect size in the population of interest estimated in the study population (Patino and Ferreira, 2015). Table 8 : Median of the mean forced outage duration pdf (hours). Median of mean forced outage duration pdf (hours) Year Upper Bound (0.975) B50% Life Lower Bound (0.025) 2015 6.239337 5.383541 4.645127 2016 5.851687 4.769084 3.886769 2017 5.286304 4.391294 3.647816 2018 6.611890 5.571840 4.695390 2019 8.435479 6.846163 5.556288 Table 9 : Median of unit failure rate pdf (outages/year). Median of unit failure rate pdf (outages/year) Year Upper Bound (0.975) B50% Life Lower Bound (0.025) 2015 9.165611 7.913893 6.848868 2016 5.821421 4.868818 4.089926 2017 6.286079 5.437703 4.703824 2018 4.483247 3.752409 3.140709 2019 4.985610 4.207730 3.551219 Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 115 3.1 Availability analysis of the thermal power plants Using the implemented distribution fitting, the B90% life was calculated for both the MFOD and the unit failure rate with a 95% Confidence Level. The aim of the distribution fitting was to evaluate with a 90% probability , the expected time needed to recover the functions of the thermal power plants and the average number of outages. Data from the KW test and distribution fittin g indicates that the MFOD increased for 2019 and the unit failure rate has been stable since 2016. Accordingly, the data from 2019 were chosen for the analysis. Tables 10 and 11 show the B90% life as computed using the quick calculation pad from Weibull++ Software. The Two-sided confidence level of 95 % was also considered. Figures 12 and 13 show the probability of restoring the power plant function and the failure rate. Table 10 : Mean forced outage duration distribution B90% life. Mean forced outage duration distribution B90% life (hours) Year Upper Bound (0.975) B90% Life Lower Bound (0.025) 2015 16.890513 13.643275 11.020325 2016 20.920577 15.590435 11.618306 2017 14.310876 11.632169 9.454861 2018 13.705005 11.928806 10.382807 2019 32.871988 24.089521 17.65348 Table 11 : Unit failure rate distribution B90% life. Unit failure rate distribution B90% life (outages/year) Year Upper Bound (0.975) B90% Life Lower Bound (0.025) 2015 19.650076 17.116064 14.922217 2016 13.809283 11.638246 9.823966 2017 12.954111 11.233481 9.741394 2018 14.224801 11.177519 8.783036 2019 15.237316 12.101755 9.611436 Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 116 Figure 12 : Probability of repair based on the mean forced outage duration data. Figure 13 : Probability of outages based on the unit failure rate data. Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 117 Using the Eq. 5 and the consolidated data from Tables 10 and 11, the FOF for 2019 was calculated, as shown in Table 12: Table 12 : Forced Outage Factor (FOF) for Brazilian thermal power plants . Forced Outage Factor (%) Year Upper Bound (0.975) FOF Lower Bound (0.025) 2019 5.72 3.33 1.94 The FOF was compared to the data from the North American Electric Reliability Corporation (NERC). NERC is a not -for-profit international regulatory authority: the Electric Reliability Organization (ERO) for North America, subject to oversight by the Federal Energy Regulatory Commission (FERC) and governmental authorities in Canada (NERC, 2020). The consolidated statistics are shown in Table 13, using data from Generating unit statistical brochure (2018), containing data on units reporting events only. Table 1 3: Forced Outage Factor (FOF) for Nort h America thermal power plants (NERC, 2018). Unit Type MW Trb/Gen # of Units FOF Fossil (all fuel types) All Sizes 1018 5.59 Gas turbine All Sizes 660 4.94 Multi -boiler/multi -turbine All Sizes 18 4.31 Jet engine All Sizes 218 2.86 Combined cycle All Sizes 281 2.19 Diesel All Sizes 142 2.13 Nuclear (all types) All Sizes 99 1.25 4. Discussion In this study, the computed KW test p -values was less than a 0.05 significance level , thus, rejecting the H0 for both KPIs, indicating that at least two population distributions differ from each other. In addition, the pairwise comparisons using the Dunn test show that: i. The computed KW for the MFOD indicator rejected the H0, indicating that at least two populati on distributions differ from each other ( P[χ²(4) > 12.5] = 0.0142 < 0.05). In addition, based on the eta- squared computing and the pairwise analysis by the Dunn test for the 0.05 significance level (p ≤ 0.05 for pairwise comparison 2017 -2019 ), there was a difference between the population distributions from 2017 and 2019 with a small size effect. Furthermore, the KW test mean rank was 188.63 and 252.80 for 2017 and 2019, respectively. ii. The computed KW for the unit failure rate indicator rejected the H0, indicating that at least two population distributions differ from each other ( P[χ²(4) > 36.5] = 0.000000231 < 0.05). In addition, based on the eta-squared computing and the pairwise analysis by the Dunn test for the 0.05 significance level (p ≤ 0.001 for pai rwise comparison 2015 -2016, p ≤ 0.05 for pairwise comparison 2015 -2017, and p ≤ 0.0001 for pairwise comparisons 2015 -2018 and 2015 -2019 ), there was a difference between the population distribution from 2015 and that from the other years with a moderate siz e effect. The KW test mean rank was 289.93 for 2015 and varied from 184.50 to 229.39 for 2016 -2020 data. This significant result in the KW test indicated that there are differences between the groups , however it does not indicate which groups (Andrew, 1998). In addition, it does not indicate whether the difference is meaningful, nor does it specify how many of the groups are different from each other (Chan and Walmsley, 1997) . If the result indicated no differences for the 0.05 significance level (p > 0.05), the stability system scenario could be verified by considering the five-years timeframe. Since the result showed significant differences, a multiple comparison Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 118 between treatments was performed to construct pair-wise multiple comparisons to identify the source of the significance (Chan and Walmsley, 1997) . According to McDonald, the KW test does not assume normality , but that the shapes of the distributions in different groups are similar (McDonald and John, 2007) . This indicates that non -parametric tests are not a good solution for heteroscedastic data. According to Laerd Statistics, the mean rank from the KW test should be considered for different shapes. With an increase in the group’s mean rank, the observation values increase in comparison to those of the other groups (Minitab 19 Support, 2020). In addition, as reported by (McDonald, 2007), the standard deviations of the measurements of different groups should always be compared to evaluate their differences . Based on the eta -squared computing, the computed KW for the MFOD indicator indicates a difference between the population distributions from 2017 and 2019 with a small size effect . The KW test mean rank for 2017 and 2019 was 188.63 and 252.80 , respectively, While the standard deviation for 2017 and 2019 was 4.89 and 15. 1 hours , respectively . These results indicate that the observation values from 2019 are higher than those from 2017 . In addition, it also indicates differences in the group distributions . Therefore, it is necessary to investigate for signs of skewness and variance, as proposed by (Fagerland and Sandvik, 2009). The effect size of the KW test was computed to verify the degree to which one group has data with higher ranks than another group. This test is rel ated to the probability that a value from one group will be greater than that from another group (Mangiafico, 2016 ). According to Prajapati et al. , the eta -squared is a measure of association , and the proportion of the total variance attributed to an effec t (Prajapati et al., 2010). Eta -squared ranges from 0 to 1 , and as a rule , 0.01 is a small effect, 0.06 is a moderate effect, and 0.14 is a large effect (Tomczak et al., 2014). Based on the eta-squared ( 𝐸𝑅2 = 0.0198 ≤ 0.06 ), the effect size for the MFOD KW test indicated that there is a difference between the population distributions with a small size effect . In addition, based on the eta-squared (0.06 < 𝐸𝑅2 = 0.0732 ≤ 0.14 ), the effect size for the unit failure rate KW test indicated that there is a di fference between the population distributions with a moderate size effect . A post hoc procedure was performed to determine which groups are different from each other (Ribeiro, 2019.; Andrew, 1998). According to Dinno, when the H0 is rejected, the Dunn’s test should follow the KW test ( Dinno, 2015) . The computed Dunn test for the 0.05 significance level of the MFOD KW test indicated a difference between the population distributions from 2017 and 2019 (p ≤ 0.05 for pairwise comparison 2017 –2019). In addition, the computed Dunn test for the 0.05 significance level for the unit failure rate KW test indicated a difference between the population distributions from 2015 and that of the other years (p ≤ 0.001 for pairwise comparison 2015 -2016, p ≤ 0.05 for pairwise comparison 2015 -2017, and p ≤ 0.0001 for pairwise comparisons 2015 -2018 and 2015 -2019). The computed distribution fitting and the distribution locations analysis show that: i. When the groups had equal variation (homoscedasticity), the distributio n fitting for the MFOD indicator confirmed a small size effect based on the computed eta -squared ( 𝐸𝑅2 = 0.0198 ≤ 0.06) and the pairwise analysis by the Dunn test for the 0.05 significance level. In addition, the distribution location (B50% life) slightly increased for the 2019 data (4.39 in 2017 and 6.84 in 2019), indicating a higher forced outage duration. ii. When the groups had equal variation (homoscedasticity), the distribution fitting for the unit failure rate indicator confirmed the moderate size effect based on the computed eta -squared (0.06 < 𝐸𝑅2 = 0.0732 ≤ 0.14) and the pairwise analysis by the Dunn test for the 0.05 significance level. In addition, the Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 119 distribution location (B50% life) decreased from the 2016 data on, indicating an overall lower unit failure rate (7.91 in 2015, 4.86 in 2016, 5.43 in 2017, 3.75 in 2018, and 4.20 in 2019 ). iii. The distribution fitting for both the MFOD and unit failure rate indicated that the data have a non- normal distributio n, thus, ranking the distributions loglogistic and generalized gama for the MFOD data and 3 -parameters Weibull, and generalized gama and lognormal distributions for the unit failure rate . iv. The graphical analysis of MFOD pdf indicates a right -skewed behavio r with an increased shifting in the 2019’s location parameter , and a slightly higher scale -parameter compared to the other years. v. The graphical analysis of unit failure rate pdf indicates a right -skewed behavior with an increased shifting in the 2015’s location -parameter and a higher scale -parameter compared to the other years. However, there are some limitations with using the KW test. According to McDonald, the KW test cann ot detect the differences between symmetrical distributions with similar location , and very different scale -parameter have different distributions . Sometimes, the KW test is considered a median test for the H0. This is because it assumes that the distribut ions in each group have similar shape, and the KW test can reject the null hypothesis even for same medians (McDonald, 2007). Another situation in which the median test for the H0 is violated is when the distributions have different degrees of skewness, wh ich affects both type I error, rejecting the null hypothesis when it is true (Fagerland and Sandvik, 2009). In addition, according to McDonald, there is no consensus about heteroscedastic data for applying a test that assumes homoscedasticity (McDonald, 2007) . The graphical analysis of the MFOD probability and unit failure rate pdfs indicate a right -skewed behavior for all the available groups, and the medians were considered for the complementary evaluation of the computed data. The Distribution Wizard from Weibull++ Software used in this study is a valuable tool for fitting distribution models. According to Reliasoft , the MLE analysis is an appropriate method for data sets with many observed failures , however, the MLE tends to be statistically distorted for small sample sizes. Therefore, the default weighted composition of three goodness of fit tests were applied to determine the rank of the distributions: (i) 40% for the Kolmogorov -Smirnov test, (ii) 10% for the Correlation coefficient test , and (iii) 50% for the Likelihood value test (Reliasoft, 2017 ). By combining the analytical and graphical tools, it was possible to verify which data pack to choose for applying the B90% life ana lysis: in this case, the 2019 data from the MFOD and unit failure rate reports provided by ONS was used . The BX% life is the lifetime metric , in which X% of the units in a population fail, when considering the reliability analysis (Woo, 2017). Applying B90 % life calculation for the MFOD data indicates a lifetime metric , in which 90% of the power units would recover to available state. Likewise, applying the B90% if calculation for the unit failure rate data indicates the failure/year metric which 90% of the power units would achieve. The choice of the year 2019 to process the availability analysis considers the following scenarios: (i) the atypical rise in the MFOD for 2019 and (ii) the stability of the unit failure rate over the past 4 years. The B90% life and unavailability computing show that: i. The FOF of the Brazilian thermal power plants is 3.33% with a 90% probability and a 95% confidence level, based on the data from the MFOD and unit failure rate data. ii. This value is slightly higher than the upper limit applied by ANEEL for hydroelectric plants (higher value is 3.115%) . iii. The FOF of the Brazilian thermal power plants is also lower than 70% of North America thermal power plants (Fossil, Gas turbine, and Mult i-boiler/multi -turbine types), based solely on the evaluated statistics. Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 120 The analysis of the operating scenario indicates that the combination of: (i) the historical second critical period in SIN (Nacional do Sistema, 2020 a) and (ii) the higher power dem and in the first quarter of 2019 (Nacional do Sistema, 2020 a) (Figure 3), have contributed to the increase in the MFOD. Furthermore, there has been no significant difference in the annual unit failure rate since 2016, assuming that the groups had equal va riation (homoscedasticity) and applying the KW test (p > 0.05). In addition, the FOF of the Brazilian thermal power plants (3.33% at 90% probability and 95% confidence level) is comparable to that from 2436 power plants from the NERC (NERC, 2018), indicati ng that the Brazilian thermal generation is a reliable system regardless of the challenges due to (i) the type of starting, (ii) the frequency of starting, and (iii) the loading pattern (Dipak, 2015a). 5. Conclusion Brazil is a continental -size country and its SIN is characterized by a marked seasonality in its electricity supply. The expansion pattern of the Brazilian electric sector shows signs of exhaustion, with the inclusion of the ROR and intermittent renewable sources . In addition, 50% of the Brazilian hydroelectric plants are over 20 years old and 32% are over 40 years old. Therefore, there is an increase in the demand for flexible thermal power plants based on availability that can be operated in electric dispatch mode. Therefore, th e knowledge acquired from the statistical analysis of forced outages is fundamental to understand the availability of Brazilian thermal power plants, to ensure the security risk mitigation to supplying the national power demand. Assuming the homoscedastic ity, the KW test and Dunn’s pairwise -comparisons are valuable non -parametric approaches for evaluating the differences in the population’s distributions for forced outage data. The tests were applied for available forced outage data provided by ONS and the results indicated that the forced outage statistics of the Brazilian thermal plants are favorable and comparable to the North America's benchmarks. Furthermore, the distribution fitting indicated right -skewed (therefore, non -normal) distributions for the available data from ONS and the KW test fits properly for this type of data , based on the distribution rank from the computed data. However, the limitations to this study are the assumption of homoscedasticity and the effect size for the KW test, which in dicates a small effect (Eta -squared ≤ 0.06) for the MFOD data. This may indicate that the sample size should be increased or separated by outage causes for applying normality tests: internal, secondary, external, and operational causes. The inflexible power generation corresponds to 75% of the Brazilian load in 2020, however, the thermal powe r generation accounts for 7% only. Therefore, strategies for including flexible thermal power units in the SIN, the logistics expansion to explore onshore natural gas in the Northern Brazil and generating power nearby reservoir (the Reservoir -to-wire – R2W ) can increase national energy security, thus, creating a virtual reservoir in the SIN that is safe and reliable for the country's development. Declaration of competing interests The author declares that he has no known competing financial interests or pe rsonal relationships that could have appeared to influence the work reported in this paper. Acknowledgement The author sincerely thanks the Brazilian National Electricity System Operator – ONS – (www.ons.org.br ) for making public valuable information about the performance of the Brazilian National Interconnected System (SIN). This research did not receive any specific grant from funding agencies in the public, commercial, or not - for-profit sectors. The author wo uld like to thank Editage ( www.editage.com ) for English language editing. Asian Institute of Research Engineering and Technology Quarterly Reviews Vol.3, No.2, 2020 121