# Risk-Based Clustering for Near Misses Identification in Integrated Deterministic and Probabilistic Safety Analysis

Authors: Francesco Di Maio; Matteo Vagnoli; Enrico Zio
Journal: Science and Technology of Nuclear Installations
 (2015)

---

## Abstract

In integrated deterministic and probabilistic safety analysis (IDPSA),
safe scenarios and prime implicants (PIs) are generated by simulation. In this paper,
we propose a novel postprocessing method, which resorts to a risk-based clustering method
for identifying Near Misses among the safe scenarios. This is important because the possibility
of recovering these combinations of failures within a tolerable grace time allows avoiding
deviations to accident and, thus, reducing the downtime (and the risk) of the system. The
postprocessing risk-significant features for the clustering are extracted from the following: (i)
the probability of a scenario to develop into an accidental scenario, (ii) the severity of the
consequences that the developing scenario would cause to the system, and (iii) the combination of
(i) and (ii) into the overall risk of the developing scenario. The optimal selection of the extracted
features is done by a wrapper approach, whereby a modified binary differential evolution (MBDE) embeds
aK-means clustering algorithm. The characteristics of
the Near Misses scenarios are identified solving a multiobjective optimization problem, using the
Hamming distance as a measure of similarity. The feasibility of the analysis is shown with
respect to fault scenarios in a dynamic steam generator (SG) of a nuclear power plant (NPP).


---

## Body

## 1. Introduction

Integrated deterministic and probabilistic safety analysis (IDPSA) attempts to overcome some limitations of deterministic safety analysis (DSA) and probabilistic safety analysis (PSA). The former is solidly founded on the multibarrier and defense-in-depth concepts and aims at verifying the capability of a nuclear power plant (NPP) to withstand a set of postulated design basis accidents (DBA) [1, 2]. To account for the uncertainties in the model representation of the actual plant behavior, conservatism is introduced in the calculations by thermal-hydraulics (TH) codes under DBA conditions [3]. The latter aims at considering a wider set of possible accidental scenarios and includes the quantification of accident probabilities [4, 5].Both DSA and PSA are scenario-based analyses, where scenario selection and definition are done by expert judgment. State of the art of DSA and PSA approaches can provide relevant and important insights into what is already known to be an “issue,” but they are not capable of revealing what, and to what extent, is not known (i.e., scenarios which are not expert-selected in the DSA and PSA inputs), with the risk of neglecting or underestimating potentially dangerous scenarios [6]. This is due to the difficulties of the static structure of the classic DSA and PSA approaches in treating dynamic variations that usually occur during the operational time of a process [7] due to (i) stochastic disturbances (e.g., equipment failures), (ii) deterministic plant responses (i.e., transients), (iii) controls, and (iv) operator actions [6, 8, 9]. Indeed, the order and timing of the events occurring along a scenario and the values of the process variables at the time of event occurrence are critical in determining the evolution of the scenario itself [10].The development and application of IDPSA in practice must meet the challenge of computational complexity, in both model construction and implementation and in postprocessing for the retrieval of the relevant information from the scenario outcomes. The number of dynamic scenario branches generated in IDPSA increases in power law with the number of occurring events and, thus, is much larger than in classical PSA based on event trees (ET) and fault trees (FT). The a posteriori information retrieval (postprocessing) then becomes quite burdensome and difficult [11, 12]. Continuous event trees (CETs) [13, 14] and dynamic event trees (DETs) [15, 16] provide realistic frameworks for IDPSA. However, their application is limited by their computationally intensive nature, by the need of tailoring the algorithms to the system under consideration and by the need of processing a massive amount of data for any single initiating event considered [17].Postprocessing, in general, consists in classifying the generated dynamic scenarios into safe scenarios and prime implicants (PIs), that is, sequences of events that represent minimal combinations of accident failures necessary for system failure and cannot be covered by more general implicants [18]. Among the safe scenarios, Near Misses are important scenarios to be identified, because they are those sequences of events that reach values of the safety parameters close to, but not exceeding, the corresponding acceptable thresholds [19]. They can, thus, be relevant contributors to the “hidden” risk of the system and should not be neglected, as a small deviation may transform them into accidental scenarios.In the literature, several authors introduce the concept of Near Misses as accident precursors [20, 21]. We here consider Near Misses as sequences of events that incidentally keep the system in a safe state but endangered and insecure. For the purpose of the analysis, they are here defined as sequences of events similar to those leading the system into fault conditions, except for one characteristic which is missing or is slightly different (e.g., sequence time lag, different failure magnitude, and different involved component in an event) [22].The postprocessing analysis entails a “Forward” classification of the dynamic scenarios into classes, that is, safe, PIs, and Near Misses and a “Backward” identification of the similarities of the features of the scenarios (i.e., stochastic event occurrence and deterministic process variables values), which characterize the groups of Near Misses among the whole set of safe scenarios.For the “Forward” classification of the Near Misses sequences, we look at two factors of risk: the probability of occurrence of an undesired event and the severity of the consequence caused by the event [23]. Thus, we describe the sequences of events by (i) the probability (
p
) that the developing scenario is an accidental scenario, (ii) the consequence (
c
) that the developing scenario can cause to the system, and (iii) the overall risk (
r
) of the developing scenario that we compute synthetically as r
=
p
×
c (expected consequence).The optimal features for discerning the Near Misses from the safe scenarios are extracted from the profiles ofp, c, and r of the accidental scenarios and selected by a wrapper algorithm, which takes into account six statistical indicators of p, c, and r, and, through a modified binary differential evolution (MBDE) optimization algorithm, selects the best features, which are fed to a K-means clustering algorithm, which is a simple and well-known clustering algorithm (other classical clustering algorithms, such as mean-shift [24, 25] or fuzzy C-means [19, 26]).The outcomes of this “Forward” classification is, then, interpreted by a “Backward” identification of the similarities of the features of the Near Misses scenarios: the acquired knowledge can be exploited in an online integrated risk monitoring system that can rapidly detect the problem and set up a repair strategy of the occurring failures before the system reaches a fault state.The proposed approach is illustrated with reference to scenarios occurring in the steam generator (SG) of a NPP [27]. We use multiple-valued logic (MVL) theory for modeling the behavior of the system, where timing and sequences of component failure events are determining the system behavior [4]. By using MVL, we increase the limited description capability of binary variables in modeling the different component operational states (e.g., a valve that can be closed, partially closed, or fully open or can fail at different times) and, therefore, perform an IDPSA postprocessing analysis on the whole set of simulated accidental scenarios [17].The paper is organized as follows. In Section2, the SG model used to generate the scenarios for the reliability analysis is presented [27], along with multistate representation of the system dynamics. In Section 3, the PIs are identified and the risk-based “Forward” and “Backward” Near Misses identification method is introduced with reference to the case study considered. In Section 4, conclusions and remarks are drawn.


## 2. Case Study

### 2.1. The U-Tube Steam Generator (UTSG) Model

The U-tube steam generator (UTSG) under consideration is sketched in Figure1. The improper control of the water level, whose difficulties arise from nonminimum phase plant characteristics, that is, plant strong inverse response behavior, particularly at low operating power, due to the so-called “swell and shrink” effects [28], is a major cause of NPP unavailability [28–30].Figure 1
Schematic of the UTSG [31].The reactor coolant enters the UTSG at the bottom and moves upward and then downward in the inverted U-tubes, transferring heat to the secondary fluid before exiting at the bottom. The secondary fluid, the feedwater (Q
e), enters the UTSG at the top of the downcomer, through the space between the tube bundle wrapper and the SG shell. The value of Q
e is regulated by a system of valves: a low flow rate valve, used when the operating power (P
o) is smaller than 15% of nominal power (P
n), and a high flow rate valve when P
o
>
0.15
P
n [27]. In the secondary side of the tube bundle, water heats up, reaches saturation, starts boiling, and turns into a two-phase mixture. The two-phase fluid moves up through the separator/riser section, where steam is separated from liquid water, and through the dryers, which ensure that the exiting steam (Q
v) is essentially dry. The separated water is recirculated back to the downcomer. The balance between the exiting Q
v and the incoming Q
e governs the change in the water level in the SG. Because of the two-phase nature, two types of water level measurements are considered, as shown in Figure 1, each reflecting a different level concept: the narrow range level (N
rl) is calculated by pressure difference between two points close to the water level and indicates the mixture level, whereas the wide range level (W
rl) is calculated by pressure difference between the two extremities of the SG (steam dome and bottom of the downcomer) and indicates the collapsed liquid level that is related to the mass of water in the SG.“Swell and shrink” phenomena are also modeled to reproduce the dynamic behavior of the SG: whenQ
v increases, the steam pressure in the steam dome decreases and the two-phase fluid in the tube bundle expands causing N
rl to initially swell (i.e., rise), instead of decreasing as would have been expected by the mass balance; contrarily, if Q
v decreases or Q
e increases, a shrink effect occurs. A similar model has been presented in [27].TheN
rl is governed by Q
e and Q
v across the tube bundle region of the SG as shown by the following transfer function:(1)
N
rl
s
=
1
T
n
s
Q
e
f
s
-
Q
G
V
s
,where Q
e
f is the flow rate of the incoming water in the tube bundle, (2), Q
G
V is the equivalent steam-water mixture flow rate exiting the tube bundle region, (3), T
n is a time constant that accounts for the N
rl dynamics.The incoming water flow rateQ
e
f is proportional to Q
e:(2)
Q
e
f
s
=
1
1
+
T
h
s
1
+
τ
s
Q
e
s
,where the lag 1
/
(
1
+
τ
s
) accounts for the feed-water valve dynamics and 1
/
(
1
+
T
h
s
) accounts for the water mass transportation dynamics: their values are reported in Table 1.Table 1
Parameters of the UTSG model at different power levels [27].
P
o
0.03
×
P
n
0.04
×
P
n
0.09
×
P
n
0.24
×
P
n
0.30
×
P
n
0.50
×
P
n
P
n
T
n
36
56
63
44
40
40
40
F
g
13
18
10
4
4
4
4
T
h
170
56
30
10
8
5
5
τ
10
10
10
30
30
30
30
T
g
10
10
10
10
10
10
10
T
int
⁡
140
140
140
140
140
140
140The exiting steam-water massQ
G
V is proportional to Q
v:(3)
Q
G
V
s
=
1
-
F
g
T
g
s
1
+
T
g
s
Q
v
s
,where the first-order lag 1
/
(
1
+
T
g
s
) accounts for the elapsed time from the turbine steam demand and the increase of Q
G
V, and the nonminimum phase term (
1
-
F
g
T
g
s
) accounts for the two-phase swell and shrink effects.Combining (1), (2), and (3), N
rl is equal to(4)
N
rl
s
=
1
T
n
s
Q
e
s
1
+
T
h
s
1
+
τ
s
-
1
-
F
g
T
g
s
1
+
T
g
s
Q
v
sand W
rl, that is, the overall water mass in the steam generator, is(5)
W
rl
s
=
1
T
int
s
Q
e
s
-
Q
v
s
,where T
int is a time constant that accounts for the W
rl dynamics.We assumey
1
=
N
rl and y
2
=
W
rl, and u
=
Q
e and d
=
Q
v; the state space representation of the SG model is, thus, (6)
x
˙
t
=
0
0
0
1
T
n
0
-
1
T
h
0
-
1
T
n
0
0
-
1
T
g
0
0
0
0
-
1
τ
x
t
+
0
0
0
1
τ
u
t
+
-
1
T
n
0
1
+
F
g
T
n
0
d
t
,
y
t
=
1
1
1
0
T
n
T
int
0
0
τ
T
int
x
t
.The values of the parameters T
h, T
n, F
g,τ, T
g, and T
int change depending on the power P
o, as shown in Table 1.The goal of the system is to maintain the SG water level at a reference position (N
ref): the SG fails if the N
rl rises (falls) above (below) the threshold N
high (N
low), in which case automatic reactor or turbine trips are triggered. Indeed, if the N
rl exceeds N
high, the steam separator and dryer lose their functionality and excessive moisture is carried in Q
v, degrading the turbine blades profile and the turbine efficiency; if N
rl decreases below N
low, insufficient cooling capability of the primary fluid occurs. Similarly, the W
rl is relevant for the cooling capability of the primary circuit [28]. Prealarms are triggered when N
rl exceeds N
h
l (N
l
l) if a small deviation from N
ref occurs or when N
rl exceeds N
v
h (N
v
l), when the deviation is large. Set points of N
ref and of N
rl depend on P
o, as shown in Figure 2, and, thus, also the alarms thresholds depend on P
o. The N
rl set point is low at low P
o, to partially account for the strong inverse response of N
rl [28]; thus, the low level thresholds are more restrictive than the high level thresholds at low P
o.Figure 2
Set point forN
rl at different power rate P
o values.A dedicated model has been implemented in SIMULINK to simulate the dynamic response of the UTSG at differentP
o values. Both feedforward and feedback digital control schemes have been adopted. The feedback controller is a PID that provides a flow rate Q
pid resulting from the residuals between N
rl and N
ref, whereas the feedforward controller operates a safety relief valve that is opened if and only if N
rl exceeds the N
h
l and removes a constant flow safety flow rate (Q
sf). The block diagram representing the SIMULINK model of the SG is shown in Figure 3: the controlled variable is N
rl, whereas the control variable is Q
e.Figure 3
Block diagram representing the SIMULINK model of the SG.


### 2.2. The Set of Possible Failures

The set of multiple component failures that can occur during the system life are shown in Figure4.(1)
The outlet steam valve can fail stuck at a random time in[
0,4000
] (s) in three different positions: (i) closed; (ii) stuck open at 50% of the nominal Q
v that should be provided at P
o; (iii) stuck open at 150% of the nominal Q
v that should be provided at P
o.
(2)
The safety relief valve can fail stuck at a random time in[
0,4000
] (s) at a uniform random value Q
sf in the range [
0.5,50.5
] (kg/s).
(3)
The communication between the sensor that monitorsN
rl and the PID controller can fail at random times in [
0,4000
] (s), in which case the PID is provided with the same input value of the previous time step.
(4)
The PID controller can fail stuck at random times in[
0,4000
] (s), providing a uniform random flow rate Q
pid belonging to [
-
18,18
]
% of the nominal Q
e that should be provided at P
o.It is worth noticing that in the UTSG there are two PID controllers and, thus, two communications between the sensors measuring N
rl and the PIDs (one for high power feedback control and the other for low power feedback control). The selective action of the PIDs depending on P
o hides some of the failures. For example, if the power profile of the scenario under investigation is a ramp, both PIDs are called in operation: if anyone (or both) failed, their fault state is detectable. On the contrary, if we consider scenarios with constant power profile, for example, low power rate (P
o
<
15
%
P
n), the occurrence of a high power feedback control failure cannot be detected, and, thus, the fault remains hidden.Figure 4
Sketch of the failures that can be injected into the system.Choices and hypotheses for modeling the failures (i.e., the mission time, the number and type of faults, the distributions of failure times, and magnitudes) have been arbitrarily made with the aim of generating multiple failures in the sequences and capturing the dynamic influence of their order, timing, and magnitude. The choice of a mission time (T
miss) equal to 4000 (s) has been made, because it is a long enough interval of time to allow the complete development also of slow dynamic accident scenarios.


### 2.3. The Multistate Representation of System Dynamics

For realistically treating the dynamic behavior of the UTSG when component failures occur, we go beyond the binary state representation and adopt a multiple value logic (MVL) [17, 32] for an approximated description of the continuous time of occurrence of component failures and their magnitude. The MVL allows describing that the components can fail at any (discrete) time (not only the initial time) along the scenario, with different (discrete) magnitudes (not only the most conservative). The discretization of the time and magnitudes values is as follows:(
i
)
time discretization: we use the labelst
mvl
=
1, t
mvl
=
2, t
mvl
=
3, and t
mvl
=
4, for failures occurring in the intervals [
0,1000
] (s), [
1001,2000
] (s), [
2001,3000
] (s), and [
3001,4000
] (s), respectively; if the label t
mvl
=
0, the component does not fail within the time of the whole scenario, T
miss;
(
ii
)
magnitude discretization:
(
a
)
the steam valve magnitude is indicated as 1, 2, or 3 for failure states corresponding to stuck at 0%, stuck at 50%, and stuck at 150% of theQ
e value that should be provided at P
o, respectively; if the steam valve magnitude is indicated as 0, the component does not fail in T
miss;
(
b
)
the safety relief valve fails with magnitude indicated as 1, 2, 3, and 4, if it is stuck between[
0.5,12.6
] (kg/s), (
12.6,25.27
] (kg/s), (
25.27,37.91
] (kg/s), and (
37.91,50.5
] (kg/s), respectively; if the safety relief valve magnitude is indicated as 0, the component does not fail in T
miss;
(
c
)
the communication between the sensor measuringN
rl and the PID controller is labelled 0 if the communication works, 1 otherwise;
(
d
)
the PID controller failure magnitude range is discretized into 8 equally spaced magnitude intervals, labelled from 1 to 8, representative of failure states corresponding to discrete intervals of output value belonging to[
-
18,18
]
% of the Q
e value that should be provided at P
o; if the PID controller magnitude is labelled as 0, the component does not fail in T
miss.The values of time and magnitude and order of failure occurrence for each component are included into a sequence vector that represents a scenario. As an example, the sequence vector of Figure 5 represents a scenario where the steam valve fails stuck at its maximum allowable value at a time in [
3001,4000
] (s) and it is the third event occurring along the sequence; the safety relief valve fails first in [
0,1000
] (s), with a magnitude belonging to [
0.5,12.6
] (kg/s); the communication between the sensor measuring N
rl and the PID controller is the second failure event in the sequence and occurs in [
2001,3000
] (s); finally, the PID controller fails stuck in [
3001,4000
] (s), with a magnitude belonging to [
6,10
]
% of the Q
e value that should be provided at P
o.Figure 5
Sequence vector representing a scenario.The number of possible sequence vectors that arise from the MVL discretization is 100509, each one evolving towards either safe or faulty conditions. To investigate this, a Monte Carlo-driven fault injection engine is used to sample combinations of discrete times and discrete magnitudes of components failures.The (dynamic) analysis has been performed with respect to the two constant power scenarios, 5%P
n (low power level) and 80% P
n (high power level). The system configurations considered are listed in Table 2.Table 2
System configurations.
Systemconfigurations
Failure of theoutlet steam valve
Failure of thesafety relief valve
Level sensor-PID controllercommunication interruption
Failure of thePID controller
1
—
—
—
—
2
X
—
—
—
3
—
X
—
—
4
—
—
X
—
5
—
—
—
X
6
X
X
—
—
7
X
—
X
—
8
X
—
—
X
9
—
X
X
—
10
—
X
—
X
11
—
—
X
X
12
X
X
X
—
13
X
X
—
X
14
X
—
X
X
15
—
X
X
X
16
X
X
X
XThe dynamic analysis shows that the same combination of components failures does not unequivocally lead to only one system end state but rather it depends on when the failures occur and with what magnitude. This is shown in Figure6, where the frequencies of occurrence of the three system end states (“High,” “Safe,” and “Low”) are plotted for the 16 dynamic system configurations of Table 2.Figure 6
Histograms for high power level (a) and low power level (b) of the frequencies of occurrence of the end states for each of the 16 system configurations of Table2, simulated by sampling discrete failure times and magnitudes of components failures.
(a)
(b)Figure7 shows that, at high power operation, the timing of the events is quite important, because with the same system configuration but different times of failure occurrences, the system end state change. Specifically, in Figure 7(a), the safety valve fails stuck at 100% of Q
sf after 1020 seconds and the communication between the sensor measuring N
rl and the PID controller fails at(i)
1052 seconds (solid line),
(ii)
1063 seconds (dashed-dotted line).The two scenarios lead to low and high failure modes, respectively, whereas they would be considered as minimal cuts sets (MCS) in a static reliability analysis presented in Appendix A.Figure 7
Example of dynamic system behavior at 80%P
n.
(a)
(b)Figure7(b) shows the effects of different failures magnitudes on the system end state: the safety relief valve fails stuck in its maximum position at 2000 seconds, the communication fails at 2010 seconds, and the PID controller fails at 2020 seconds with two different magnitudes:(i)
magnitude equal to 13% of the nominalQ
e value that should be provided at 80% P
n (dashed-dotted line),
(ii)
magnitude equal to 12% of the nominalQ
e value that should be provided at 80% P
n (solid line).The low power scenarios also present dynamic effects, as shown in Figure 8. In particular, Figure 8(a) shows the effects of the timing on the system end state: the safety relief valve fails stuck at Q
sf
=
50.5 (kg/s) at 1005 (s) and the steam output valve fails stuck at 150% of the nominal Q
v value that should be provided at 5% P
n at(i)
1046 seconds (dashed-dotted line),
(ii)
1047 seconds (solid line).Figure 8(b) shows the effects of the order of components failure occurrence on the system end state: the safety relief valve fails stuck at Q
sf
=
50.5 (kg/s) and the PID controller fails stuck at its minimum allowable value.(i)
The PID controller failure is the first failure event along the sequence of events (dashed-dotted line).
(ii)
The safety relief valve failure is the first failure event along the sequence of events (solid line).Hereafter, without loss of generality, among the system configurations of Table 2, we focus only on the classification of the PIs and Near Misses of the high level failure mode at high power level (P
o
=
80
%
P
n).Figure 8
Example of dynamic system behavior at 5%P
n.
(a)
(b)


## 2.1. The U-Tube Steam Generator (UTSG) Model

The U-tube steam generator (UTSG) under consideration is sketched in Figure1. The improper control of the water level, whose difficulties arise from nonminimum phase plant characteristics, that is, plant strong inverse response behavior, particularly at low operating power, due to the so-called “swell and shrink” effects [28], is a major cause of NPP unavailability [28–30].Figure 1
Schematic of the UTSG [31].The reactor coolant enters the UTSG at the bottom and moves upward and then downward in the inverted U-tubes, transferring heat to the secondary fluid before exiting at the bottom. The secondary fluid, the feedwater (Q
e), enters the UTSG at the top of the downcomer, through the space between the tube bundle wrapper and the SG shell. The value of Q
e is regulated by a system of valves: a low flow rate valve, used when the operating power (P
o) is smaller than 15% of nominal power (P
n), and a high flow rate valve when P
o
>
0.15
P
n [27]. In the secondary side of the tube bundle, water heats up, reaches saturation, starts boiling, and turns into a two-phase mixture. The two-phase fluid moves up through the separator/riser section, where steam is separated from liquid water, and through the dryers, which ensure that the exiting steam (Q
v) is essentially dry. The separated water is recirculated back to the downcomer. The balance between the exiting Q
v and the incoming Q
e governs the change in the water level in the SG. Because of the two-phase nature, two types of water level measurements are considered, as shown in Figure 1, each reflecting a different level concept: the narrow range level (N
rl) is calculated by pressure difference between two points close to the water level and indicates the mixture level, whereas the wide range level (W
rl) is calculated by pressure difference between the two extremities of the SG (steam dome and bottom of the downcomer) and indicates the collapsed liquid level that is related to the mass of water in the SG.“Swell and shrink” phenomena are also modeled to reproduce the dynamic behavior of the SG: whenQ
v increases, the steam pressure in the steam dome decreases and the two-phase fluid in the tube bundle expands causing N
rl to initially swell (i.e., rise), instead of decreasing as would have been expected by the mass balance; contrarily, if Q
v decreases or Q
e increases, a shrink effect occurs. A similar model has been presented in [27].TheN
rl is governed by Q
e and Q
v across the tube bundle region of the SG as shown by the following transfer function:(1)
N
rl
s
=
1
T
n
s
Q
e
f
s
-
Q
G
V
s
,where Q
e
f is the flow rate of the incoming water in the tube bundle, (2), Q
G
V is the equivalent steam-water mixture flow rate exiting the tube bundle region, (3), T
n is a time constant that accounts for the N
rl dynamics.The incoming water flow rateQ
e
f is proportional to Q
e:(2)
Q
e
f
s
=
1
1
+
T
h
s
1
+
τ
s
Q
e
s
,where the lag 1
/
(
1
+
τ
s
) accounts for the feed-water valve dynamics and 1
/
(
1
+
T
h
s
) accounts for the water mass transportation dynamics: their values are reported in Table 1.Table 1
Parameters of the UTSG model at different power levels [27].
P
o
0.03
×
P
n
0.04
×
P
n
0.09
×
P
n
0.24
×
P
n
0.30
×
P
n
0.50
×
P
n
P
n
T
n
36
56
63
44
40
40
40
F
g
13
18
10
4
4
4
4
T
h
170
56
30
10
8
5
5
τ
10
10
10
30
30
30
30
T
g
10
10
10
10
10
10
10
T
int
⁡
140
140
140
140
140
140
140The exiting steam-water massQ
G
V is proportional to Q
v:(3)
Q
G
V
s
=
1
-
F
g
T
g
s
1
+
T
g
s
Q
v
s
,where the first-order lag 1
/
(
1
+
T
g
s
) accounts for the elapsed time from the turbine steam demand and the increase of Q
G
V, and the nonminimum phase term (
1
-
F
g
T
g
s
) accounts for the two-phase swell and shrink effects.Combining (1), (2), and (3), N
rl is equal to(4)
N
rl
s
=
1
T
n
s
Q
e
s
1
+
T
h
s
1
+
τ
s
-
1
-
F
g
T
g
s
1
+
T
g
s
Q
v
sand W
rl, that is, the overall water mass in the steam generator, is(5)
W
rl
s
=
1
T
int
s
Q
e
s
-
Q
v
s
,where T
int is a time constant that accounts for the W
rl dynamics.We assumey
1
=
N
rl and y
2
=
W
rl, and u
=
Q
e and d
=
Q
v; the state space representation of the SG model is, thus, (6)
x
˙
t
=
0
0
0
1
T
n
0
-
1
T
h
0
-
1
T
n
0
0
-
1
T
g
0
0
0
0
-
1
τ
x
t
+
0
0
0
1
τ
u
t
+
-
1
T
n
0
1
+
F
g
T
n
0
d
t
,
y
t
=
1
1
1
0
T
n
T
int
0
0
τ
T
int
x
t
.The values of the parameters T
h, T
n, F
g,τ, T
g, and T
int change depending on the power P
o, as shown in Table 1.The goal of the system is to maintain the SG water level at a reference position (N
ref): the SG fails if the N
rl rises (falls) above (below) the threshold N
high (N
low), in which case automatic reactor or turbine trips are triggered. Indeed, if the N
rl exceeds N
high, the steam separator and dryer lose their functionality and excessive moisture is carried in Q
v, degrading the turbine blades profile and the turbine efficiency; if N
rl decreases below N
low, insufficient cooling capability of the primary fluid occurs. Similarly, the W
rl is relevant for the cooling capability of the primary circuit [28]. Prealarms are triggered when N
rl exceeds N
h
l (N
l
l) if a small deviation from N
ref occurs or when N
rl exceeds N
v
h (N
v
l), when the deviation is large. Set points of N
ref and of N
rl depend on P
o, as shown in Figure 2, and, thus, also the alarms thresholds depend on P
o. The N
rl set point is low at low P
o, to partially account for the strong inverse response of N
rl [28]; thus, the low level thresholds are more restrictive than the high level thresholds at low P
o.Figure 2
Set point forN
rl at different power rate P
o values.A dedicated model has been implemented in SIMULINK to simulate the dynamic response of the UTSG at differentP
o values. Both feedforward and feedback digital control schemes have been adopted. The feedback controller is a PID that provides a flow rate Q
pid resulting from the residuals between N
rl and N
ref, whereas the feedforward controller operates a safety relief valve that is opened if and only if N
rl exceeds the N
h
l and removes a constant flow safety flow rate (Q
sf). The block diagram representing the SIMULINK model of the SG is shown in Figure 3: the controlled variable is N
rl, whereas the control variable is Q
e.Figure 3
Block diagram representing the SIMULINK model of the SG.


## 2.2. The Set of Possible Failures

The set of multiple component failures that can occur during the system life are shown in Figure4.(1)
The outlet steam valve can fail stuck at a random time in[
0,4000
] (s) in three different positions: (i) closed; (ii) stuck open at 50% of the nominal Q
v that should be provided at P
o; (iii) stuck open at 150% of the nominal Q
v that should be provided at P
o.
(2)
The safety relief valve can fail stuck at a random time in[
0,4000
] (s) at a uniform random value Q
sf in the range [
0.5,50.5
] (kg/s).
(3)
The communication between the sensor that monitorsN
rl and the PID controller can fail at random times in [
0,4000
] (s), in which case the PID is provided with the same input value of the previous time step.
(4)
The PID controller can fail stuck at random times in[
0,4000
] (s), providing a uniform random flow rate Q
pid belonging to [
-
18,18
]
% of the nominal Q
e that should be provided at P
o.It is worth noticing that in the UTSG there are two PID controllers and, thus, two communications between the sensors measuring N
rl and the PIDs (one for high power feedback control and the other for low power feedback control). The selective action of the PIDs depending on P
o hides some of the failures. For example, if the power profile of the scenario under investigation is a ramp, both PIDs are called in operation: if anyone (or both) failed, their fault state is detectable. On the contrary, if we consider scenarios with constant power profile, for example, low power rate (P
o
<
15
%
P
n), the occurrence of a high power feedback control failure cannot be detected, and, thus, the fault remains hidden.Figure 4
Sketch of the failures that can be injected into the system.Choices and hypotheses for modeling the failures (i.e., the mission time, the number and type of faults, the distributions of failure times, and magnitudes) have been arbitrarily made with the aim of generating multiple failures in the sequences and capturing the dynamic influence of their order, timing, and magnitude. The choice of a mission time (T
miss) equal to 4000 (s) has been made, because it is a long enough interval of time to allow the complete development also of slow dynamic accident scenarios.


## 2.3. The Multistate Representation of System Dynamics

For realistically treating the dynamic behavior of the UTSG when component failures occur, we go beyond the binary state representation and adopt a multiple value logic (MVL) [17, 32] for an approximated description of the continuous time of occurrence of component failures and their magnitude. The MVL allows describing that the components can fail at any (discrete) time (not only the initial time) along the scenario, with different (discrete) magnitudes (not only the most conservative). The discretization of the time and magnitudes values is as follows:(
i
)
time discretization: we use the labelst
mvl
=
1, t
mvl
=
2, t
mvl
=
3, and t
mvl
=
4, for failures occurring in the intervals [
0,1000
] (s), [
1001,2000
] (s), [
2001,3000
] (s), and [
3001,4000
] (s), respectively; if the label t
mvl
=
0, the component does not fail within the time of the whole scenario, T
miss;
(
ii
)
magnitude discretization:
(
a
)
the steam valve magnitude is indicated as 1, 2, or 3 for failure states corresponding to stuck at 0%, stuck at 50%, and stuck at 150% of theQ
e value that should be provided at P
o, respectively; if the steam valve magnitude is indicated as 0, the component does not fail in T
miss;
(
b
)
the safety relief valve fails with magnitude indicated as 1, 2, 3, and 4, if it is stuck between[
0.5,12.6
] (kg/s), (
12.6,25.27
] (kg/s), (
25.27,37.91
] (kg/s), and (
37.91,50.5
] (kg/s), respectively; if the safety relief valve magnitude is indicated as 0, the component does not fail in T
miss;
(
c
)
the communication between the sensor measuringN
rl and the PID controller is labelled 0 if the communication works, 1 otherwise;
(
d
)
the PID controller failure magnitude range is discretized into 8 equally spaced magnitude intervals, labelled from 1 to 8, representative of failure states corresponding to discrete intervals of output value belonging to[
-
18,18
]
% of the Q
e value that should be provided at P
o; if the PID controller magnitude is labelled as 0, the component does not fail in T
miss.The values of time and magnitude and order of failure occurrence for each component are included into a sequence vector that represents a scenario. As an example, the sequence vector of Figure 5 represents a scenario where the steam valve fails stuck at its maximum allowable value at a time in [
3001,4000
] (s) and it is the third event occurring along the sequence; the safety relief valve fails first in [
0,1000
] (s), with a magnitude belonging to [
0.5,12.6
] (kg/s); the communication between the sensor measuring N
rl and the PID controller is the second failure event in the sequence and occurs in [
2001,3000
] (s); finally, the PID controller fails stuck in [
3001,4000
] (s), with a magnitude belonging to [
6,10
]
% of the Q
e value that should be provided at P
o.Figure 5
Sequence vector representing a scenario.The number of possible sequence vectors that arise from the MVL discretization is 100509, each one evolving towards either safe or faulty conditions. To investigate this, a Monte Carlo-driven fault injection engine is used to sample combinations of discrete times and discrete magnitudes of components failures.The (dynamic) analysis has been performed with respect to the two constant power scenarios, 5%P
n (low power level) and 80% P
n (high power level). The system configurations considered are listed in Table 2.Table 2
System configurations.
Systemconfigurations
Failure of theoutlet steam valve
Failure of thesafety relief valve
Level sensor-PID controllercommunication interruption
Failure of thePID controller
1
—
—
—
—
2
X
—
—
—
3
—
X
—
—
4
—
—
X
—
5
—
—
—
X
6
X
X
—
—
7
X
—
X
—
8
X
—
—
X
9
—
X
X
—
10
—
X
—
X
11
—
—
X
X
12
X
X
X
—
13
X
X
—
X
14
X
—
X
X
15
—
X
X
X
16
X
X
X
XThe dynamic analysis shows that the same combination of components failures does not unequivocally lead to only one system end state but rather it depends on when the failures occur and with what magnitude. This is shown in Figure6, where the frequencies of occurrence of the three system end states (“High,” “Safe,” and “Low”) are plotted for the 16 dynamic system configurations of Table 2.Figure 6
Histograms for high power level (a) and low power level (b) of the frequencies of occurrence of the end states for each of the 16 system configurations of Table2, simulated by sampling discrete failure times and magnitudes of components failures.
(a)
(b)Figure7 shows that, at high power operation, the timing of the events is quite important, because with the same system configuration but different times of failure occurrences, the system end state change. Specifically, in Figure 7(a), the safety valve fails stuck at 100% of Q
sf after 1020 seconds and the communication between the sensor measuring N
rl and the PID controller fails at(i)
1052 seconds (solid line),
(ii)
1063 seconds (dashed-dotted line).The two scenarios lead to low and high failure modes, respectively, whereas they would be considered as minimal cuts sets (MCS) in a static reliability analysis presented in Appendix A.Figure 7
Example of dynamic system behavior at 80%P
n.
(a)
(b)Figure7(b) shows the effects of different failures magnitudes on the system end state: the safety relief valve fails stuck in its maximum position at 2000 seconds, the communication fails at 2010 seconds, and the PID controller fails at 2020 seconds with two different magnitudes:(i)
magnitude equal to 13% of the nominalQ
e value that should be provided at 80% P
n (dashed-dotted line),
(ii)
magnitude equal to 12% of the nominalQ
e value that should be provided at 80% P
n (solid line).The low power scenarios also present dynamic effects, as shown in Figure 8. In particular, Figure 8(a) shows the effects of the timing on the system end state: the safety relief valve fails stuck at Q
sf
=
50.5 (kg/s) at 1005 (s) and the steam output valve fails stuck at 150% of the nominal Q
v value that should be provided at 5% P
n at(i)
1046 seconds (dashed-dotted line),
(ii)
1047 seconds (solid line).Figure 8(b) shows the effects of the order of components failure occurrence on the system end state: the safety relief valve fails stuck at Q
sf
=
50.5 (kg/s) and the PID controller fails stuck at its minimum allowable value.(i)
The PID controller failure is the first failure event along the sequence of events (dashed-dotted line).
(ii)
The safety relief valve failure is the first failure event along the sequence of events (solid line).Hereafter, without loss of generality, among the system configurations of Table 2, we focus only on the classification of the PIs and Near Misses of the high level failure mode at high power level (P
o
=
80
%
P
n).Figure 8
Example of dynamic system behavior at 5%P
n.
(a)
(b)


## 3. Near Misses Identification

The Near Misses identification is here treated as a classification problem, in which Near Misses are sorted out from the safe scenarios, among the whole set of accidental transients simulated. In practice, the PIs are first identified among the whole set of 100509 possible scenarios and, then, the Near Misses are separated out among the remaining safe scenarios.


### 3.1. Prime Implicants Identification

A PI is a set of variables that represents a minimal combination of accident component failures necessary for system failure and cannot be covered by a more reduced implicant [17, 18]. Note that in our case the “PIs” identification task may consider noncoherent structure functions, for which both failed and working states of the same components can lead the system to failure. In such circumstances, traditional methods, for example, based on minimal cut sets analysis, cannot be applied, whereas dynamic reliability methods need to be applied for the identification of the PIs [33, 34].The PIs identification among the whole set of 100509 possible scenarios is performed by means of the visual interactive method presented in [34]. The basic idea it relies on is that PIs are those scenarios with as few as possible events that are capable of leading the system into a failure state [35]; then, we first select as most important feature for the PIs identification the literal cost of the sequence vector (i.e., the number of components whose behavior is specified in the accident sequence) and then the accident sequences associated with the lowest literal cost are selected and stored as PIs. In fact, these are the most reduced sequences (i.e., with least number of events) that cannot be covered by any other implicant, and, thus, these are PIs by definition. The selected PIs and the implicants covered by them are deleted from the set of implicants and the procedure is repeated for the remaining implicants until all are covered. By so doing, 1255 PIs are identified for the high level failure mode, covering 36128 minterms. The total computational time approximately required for the identification of the PIs is 780 (s) on an Intel Core 2 Duo T9300 CPU @2.50 GHz.


### 3.2. The “Forward” Classification

Once the (1255) PIs for the SG high level failure mode have been identified, they are removed from the set of all possible scenarios, which is left with 64381 safe scenarios. For the identification of Near Misses among these, we resort to their definition as sequences of failure events that indeed keep the system in a safe condition but endangered (i.e., a quasifault system state). To this aim, we introduce a risk-based characterization of these remaining scenarios, calculating their associated risk, at each time instantt, as [23](7)
Risk
t
=
p
t
×
c
t
,where p
t is the probability that at time t the scenario can lead the system into an accidental scenario and c
(
t
) is the consequence that the developing scenario would cause to the system.In this view, we build a functional relationship such thatp increases as N
rl moves further away from the reference level N
ref, in a way that p
=
0 if N
rl is equal to N
ref and p
=
1 if N
rl reaches N
high. Such relationship is given in (8) below, assuming that scenarios whose N
rl
(
t
) approaches N
high are more prone to failure than those with N
rl
(
t
) close to N
ref; that is, (8) “filters out” (i.e., neglects) scenarios whose N
rl
t is close to N
ref and “mines” (i.e., weighs more) scenarios whose N
rl
t is close to N
high:(8)
p
t
=
φ
N
rl
t
-
μ
+
5
σ
σ
=
∫
N
ref
N
rl
t
1
2
π
σ
e
N
rl
t
-
μ
+
5
σ
2
/
2
σ
2
d
N
rl
,where φ is the cumulative probability function of the Gaussian distribution with mean μ
=
N
ref, and standard deviation σ
=
(
N
high
-
N
ref
)
/
5. Figure 9 shows the trend of p
(
t
).Figure 9
Probability functionp
(
t
) for the definition of risk.The consequencec
(
t
) of a scenario increases as N
rl approaches the failure threshold N
high, and c
(
t
) can be calculated at time t as [23] (9)
c
t
=
A
NRL
t
-
μ
+
3
σ
/
NRL
t
-
μ
,where A is the intensity coefficient that accounts for the closeness of N
rl to the thresholds N
h
l, N
v
h, and N
high, and for the exceedance time between the first event of the failure sequence (hereafter called initiating event (IE)) and the time of exceeding the threshold: the shorter this time, the more critical the scenario. Thus, the larger A is, the faster and closer N
rl approaches a threshold; we assume A
=
100 (no consequences) if no threshold is exceeded; A
=
200 (low consequences) if N
rl exceeds N
h
l after at least 2001 (s) from IE or if N
rl exceeds N
v
h after at least 3001 (s) from IE; A
=
300 (medium consequences) if N
rl exceeds N
h
l within 2000 (s) from IE, if N
rl exceeds N
v
h and the elapsed time is in [
1001,3000
] (s), and if N
rl exceeds N
high after at least 2001 (s) from IE; A
=
400 (catastrophic consequences) if N
rl exceeds N
v
h within 1000 (s) from IE or if N
rl exceeds N
high and the elapsed time from IE is in [
1,2000
] (s). A matrix representation of the intensity coefficient is shown in Figure 10.Figure 10
Matrix representation of the intensity coefficientA.By so doing, the available 64381 remaining safe scenarios are fully described at each time instantt
=
1,2
,
…
,
4000 [s] by their values of probability p
(
t
), consequence c
t, and overall risk r
(
t
). An example of the p
(
t
), c
(
t
), and r
(
t
) evolutions for two generic trends of N
rl
(
t
) is shown in Figure 11. More specifically, the N
rl
(
t
) behaviors represented in Figure 11(a) are due to(i)
solid line: the PID controller fails at 100 (s) with magnitude 4 and the safety relief valve fails at 190 (s) with magnitude 2;
(ii)
dashed-dotted line: the safety relief valve fails at 100 (s) with magnitude 1, the communication between the sensor measuringN
rl and the PID controller is interrupted at 136 (s), and the PID controller fails at 3917 (s) with magnitude 5.It is worth analysing the behavior of p
(
t
), c
(
t
), and, thus, r
(
t
) (Figures 11(b), 11(c), and 11(d), resp.): all three abovementioned functions increase as N
rl
(
t
) moves further away from N
ref and decrease as N
rl
(
t
) approaches N
ref. The steps shown in the consequences and risk plots (around 800 [s] for the solid line scenario and around 3500 [s] for the dashed-dotted line scenario) are due to the change of the discrete consequence intensity coefficient A along the scenarios. The solid line scenario is faster than the dashed-dotted line scenario (upper plot) and, thus, the value of the parameters A for the former scenario is 400 (catastrophic consequences, see Figure 10), due to the fact that N
rl
(
t
) exceeds N
v
h within 1000 (s), whereas, A
=
300 (medium consequences, see Figure 10) for the dashed-dotted scenarios, because N
rl
(
t
) exceeds N
v
h within [
1001,3000
] (s). Thus, the solid line scenario is more abrupt in its development towards failure and expected to have more catastrophic consequences, and, thus, more overall risk, than the dashed-dotted scenario, because the time between IE and the exceedance of N
v
h is shorter (i.e, less grace time).Figure 11
Probabilityp
(
t
), consequences c
(
t
), and risk r
(
t
) for two sequences of events.
(a)
(b)
(c)
(d)


#### 3.2.1. Features Selection

The identification of the Near Misses is treated as an unsupervised classification problem and addressed by clustering, where (i) the number of clusters is unknown and (ii) the features that enable the best clustering according to the risk-based characteristic profiles ofp
(
t
), c
t, and r
(
t
) of the accidental scenarios are unknown. Unsupervised clustering, thus, entails identifying the number K of clusters in which similar scenarios can be grouped according to similar values of some scenario features. To do this, from the profiles p
(
t
), c
t, and r
(
t
), we extract some statistical indicators as features [36]:(1)
mean value:(10)
μ
n
=
1
T
miss
∑
t
=
1
N
n
t
,
(2)
peak value:(11)
max
⁡
=
max
⁡
t
=
1,2
,
…
,
T
miss
n
t
,
(3)
standard deviation:(12)
σ
n
=
1
T
miss
-
1
∑
t
=
1
T
miss
n
t
-
μ
2
,
(4)
root mean square:(13)
RMS
=
1
T
miss
∑
t
=
1
T
miss
n
t
2
,
(5)
skewness:(14)
SK
=
∑
t
=
1
T
miss
n
t
-
μ
3
T
miss
-
1
,
(6)
kurtosis:(15)
KU
=
∑
t
=
1
T
miss
n
t
-
μ
4
T
miss
-
1
,where n
(
t
) is alternatively equal to p
(
t
), c
t, and r
(
t
) and, thus, the total number of features is equal to 6
×
3
=
18. Among these 18 available features, we search for those that are optimal for clustering the 64381 scenarios in Near Misses and safe scenarios.We resort to a wrapper framework [37, 38], whereby a modified binary differential evolution (MBDE) search engine [33, 39] searches candidate groups of features sets that are fed to a K-means clustering algorithm [40]; eventually, the wrapper evolves so that among these candidate groups, the group retained is that which makes the K-means clustering algorithm perform best (most compact and separate clusters). The idea behind the wrapper approach is shown in Figure 12. During the features search by MBDE, the K-means clustering is run on the N
=
0.80
×
64381
=
51505 (training) safe scenarios with sets of features (F) that are randomly selected by the MBDE algorithm. The optimal number (K) of clusters is also unknown and it is determined by looking at the clustering performance obtained by the K-means with reference to the Calinski-Harabasz (CH) index [41], which accounts for the ratio of the overall between-cluster variance (separation) and the overall within-cluster variance (compactness). The search proceeds iteratively until the CH index is maximised and the number of clusters K is fixed. Then, the results of the wrapper algorithm are evaluated on an independent test set (T
set), that is, the 0.2
×
643281
=
12876 safe scenarios that have been left out during the training phase.Figure 12
Wrapper approach for optimal feature subset selection based on a MBDE optimization algorithm and aK-means classifier.The CH index for a numberK of clusters, k
=
1,2
,
…
,
K is equal to [41](16)
CH
=
S
S
b
S
S
w
×
N
×
F
-
K
K
-
1
,where S
S
b is the overall between-cluster variance:(17)
S
S
b
=
∑
k
=
1
K
n
k
m
k
-
m
2
,and S
S
w is the overall within-cluster variance:(18)
S
S
w
=
∑
k
=
1
K
i
∑
x
∈
k
x
c
-
m
k
2
,where n
k is the number of scenarios x
c assigned to the kth cluster, m
k is the centroid of the kth cluster, that is, the mean of the selected features belonging to the kth cluster, m is the mean of the selected features, and m
k
-
m
2 and x
c
-
m
k
2 are the L
2 norms, that is, Euclidean distances, between the two vectors.The optimal features selection provides as best features the standard deviation ofc
t, the standard deviation of r
t, and the root mean square of r
t; the best performance is obtained with CH
=
9.35
e
+
04 and K
=
5.


#### 3.2.2. The Clustering Results

TheK
=
5 obtained clusters of the safe scenarios are shown in Figure 13 with reference to the features of mean risk (μ
risk) and time elapsed from the instant t
risk at which r
t starts to deviate from zero, that is, the time interval during which the system is exposed to risk. The rationale behind this choice is that the larger μ
risk and the longer t
risk, the more dangerous the scenarios. In Figure 13, clusters 3, 4, and 5 (triangles, crosses, and squares, resp.) are well separated, that is, the low level risk scenarios clusters are widened by the adoption of (8) for the quantification of the risk profile r
t. It is possible to distinguish the scenarios having the lowest risk level from the scenarios having low risk level, and, thus, the highest risk scenarios are well separated from the lower risk scenarios. The good performance obtained when (8) is adopted instead of other p
t profiles, for example, linear probability function (p
t
∝
N
rl
t) that would give the same importance to any level N
rl, for the quantification of the risk profile r
t, is due to the fact that (8) “filters out” (i.e., neglects) scenarios whose N
rl
t is close to N
ref and “mines” (i.e., weighs more) scenarios whose N
rl
t is close to N
high: the 332 circles in Figure 13 (listed in Appendix B) can, thus, be considered the Near Misses scenarios, that is, scenarios that incidentally keep the system into safe state, although in endangered and insecure operational conditions.Figure 13
Clustering results.


### 3.3. The “Backward” Approach

Once the Near Misses for the SG high level failure mode have been identified by clustering, we can search for similarities among them in terms of their multiple value sequences, that is, order and timing of event occurrences and deterministic process variables values. This “Backward” approach can lead us to finding the minimum conditions, that is, minimumμ
risk and minimum t
risk, that lead the system into a quasifault state. The problem can be framed as a multiobjective optimization problem (MOP) [42] that looks for the set of scenarios x
- to dominate any other scenarios with respect to the fitness function f:(19)
f
x
=
f
1
μ
risk
,
f
2
t
risk
,where f
1 and f
2 are the objectives functions of the defined MOP, that is, minimum μ
risk and t
risk, respectively. The solution of the MOP of (19) is the Pareto set shown in Figure 14, where 12 solutions are plotted (squares lined by continuous line) and listed in Table 3. These scenarios x
- of minimum f
x are expected to cover all failure of Near Misses scenarios cluster.Table 3
List of the Pareto-optimalx
- sequence vectors.
x
-
Steam valve failure time
Steam valve failure magnitude
Steam valve failure order
Safety valve failure time
Safety valve failure magnitude
Safety valve failure order
Sensor-PID failure time
Sensor-PID failure magnitude
Sensor-PID failure order
PID failure time
PID failure magnitude
PID failure order
1
0
0
0
2
3
1
4
1
2
4
3
3
2
3
1
4
1
2
3
1
1
2
1
4
1
3
3
1
4
2
2
3
2
1
1
2
4
2
4
3
1
3
2
2
2
4
1
4
1
4
1
5
3
1
4
3
3
3
2
1
2
1
4
1
6
3
2
4
1
2
2
1
1
3
1
4
1
7
3
2
4
2
3
1
2
1
2
3
3
3
8
4
1
3
3
2
2
4
1
4
2
4
1
9
4
2
3
1
3
2
0
0
0
1
3
1
10
4
2
4
2
1
3
2
1
1
2
4
2
11
4
2
3
4
2
2
0
0
0
4
4
1
12
4
2
4
4
2
3
2
1
1
4
4
2Figure 14
Pareto front for the cluster of Near Misses.The coverage can be verified by, first, identifying the most similar characteristics of the sequence vectors belonging to the Near Misses cluster with the Pareto set scenariosx
-, and, then, by solving a set covering problem (SCP) [43, 44].The most similar characteristics can be computed by coverage vectors (one for each scenario belonging tox
-): this entails calculating the Hamming distance [45] between each sequence vectors in x
- and each one of the other sequence vectors in the Near Misses cluster [46]. The entries of the coverage vector (in our case twelve entries, one for time, magnitude, and order of occurrence of each component failure; see Figure 5) are increased if the Hamming distance between one same entry of the considered scenario belonging to x
- and of the Near Misses vectors is equal to zero, as shown, without loss of generality, in Figure 15 for 1 sequence vector of x
- and only 2 Near Misses vectors.Figure 15
Coverage vector computation by Hamming distance.Table4 lists the 12 coverage vectors, where each entry is the percentage of Near Misses vectors having the same stochastic behavior of the optimal set x
- shown in Table 3. It can be seen that, for each scenario belonging to x
-, columns 8, 11, and 12 (e.g., sensor-PID communication failure magnitude, PID failure magnitude, and PID order of failure, resp.) have the largest values of the coverage vectors: this means that the majority of the sequence vectors of the Near Misses clusters can be well represented by (only) these failures. Furthermore, the analysis of the MVL values of the scenarios belonging to x
- (Table 3) where the largest coverage values of these columns are registered (i.e, 87%, 98.5%, and 85.2% for columns 8, 11, and 12, resp.) highlights that these failures are characterized by the same MVL values that can be summarized as follows:(i)
the failure of the communication between the sensor monitoring theN
rl and the PID controller;
(ii)
the failure of the PID controller with magnitude belonging to[
-
5
,
-
1
]
% of the Q
e value that should be provided at P
o, that is, magnitude equal to 4 in MVL framework, and it is the first accident occurring along the sequence of events in over 85% of the Near Misses scenarios.A SCP can, thus, be solved for verifying that these latest characteristics are the minimum set of stochastic event occurrences and deterministic process variables values of x
- that exhaustively describe the scenarios belonging to the Near Misses cluster: if a Near Miss sequence vector is characterized by (at least) one of the common characteristics, this is covered by the optimal set x
-. In the present application we have verified that all the scenarios belonging to the identified Near Misses cluster are covered by the minimal conditions that lead the system into a quasifault state, that is, the optimal set x
-. In conclusion, the occurrence of one of the common characteristics listed above is sufficient to lead the system into endangered and insecure operational conditions.Table 4
List of coverage vectors for each scenario belonging to the Pareto setx
-.
x
-
Steam valve failure time
Steam valve failure magnitude
Steam valve failure order
Safety valve failure time
Safety valve failure magnitude
Safety valve failure order
Sensor-PID failure time
Sensor-PID failure magnitude
Sensor-PID failure order
PID failure time
PID failure magnitude
PID failure order
1
11.4
11.4
11.4
8.1
30.4
0.6
25.3
87
46.4
1.2
1.5
0.6
2
30.1
44.3
27.1
1.2
19.9
38
19.6
87
46.4
68.1
98.5
85.2
3
30.1
44.3
27.1
8.1
19.9
38
22.3
87
14.2
27.4
98.5
14.2
4
30.1
44.3
42.2
8.1
19.9
19.9
25.3
87
12.7
68.1
98.5
85.2
5
30.1
44.3
27.1
28
30.4
38
22.3
87
46.4
68.1
98.5
85.2
6
30.1
36.1
27.1
1.2
19.9
19.9
19.6
87
13.9
68.1
98.5
85.2
7
30.1
36.1
27.1
8.1
30.4
0.6
22.3
87
46.4
3.3
1.5
0.60
8
52.4
44.3
42.2
28
19.9
19.9
25.3
87
12.7
27.4
98.5
85.2
9
52.4
36.1
42.2
1.2
30.4
19.9
13
13
13
68.1
1.5
85.2
10
52.4
36.1
27.1
8.1
9.04
38
22.3
87
14.2
27.4
98.5
14.2
11
52.4
36.1
42.2
51.2
19.9
19.9
13
13
13
1.2
98.5
85.2
12
52.4
36.1
27.1
51.2
19.9
38
22.3
87
14.2
1.2
98.5
14.2


## 3.1. Prime Implicants Identification

A PI is a set of variables that represents a minimal combination of accident component failures necessary for system failure and cannot be covered by a more reduced implicant [17, 18]. Note that in our case the “PIs” identification task may consider noncoherent structure functions, for which both failed and working states of the same components can lead the system to failure. In such circumstances, traditional methods, for example, based on minimal cut sets analysis, cannot be applied, whereas dynamic reliability methods need to be applied for the identification of the PIs [33, 34].The PIs identification among the whole set of 100509 possible scenarios is performed by means of the visual interactive method presented in [34]. The basic idea it relies on is that PIs are those scenarios with as few as possible events that are capable of leading the system into a failure state [35]; then, we first select as most important feature for the PIs identification the literal cost of the sequence vector (i.e., the number of components whose behavior is specified in the accident sequence) and then the accident sequences associated with the lowest literal cost are selected and stored as PIs. In fact, these are the most reduced sequences (i.e., with least number of events) that cannot be covered by any other implicant, and, thus, these are PIs by definition. The selected PIs and the implicants covered by them are deleted from the set of implicants and the procedure is repeated for the remaining implicants until all are covered. By so doing, 1255 PIs are identified for the high level failure mode, covering 36128 minterms. The total computational time approximately required for the identification of the PIs is 780 (s) on an Intel Core 2 Duo T9300 CPU @2.50 GHz.


## 3.2. The “Forward” Classification

Once the (1255) PIs for the SG high level failure mode have been identified, they are removed from the set of all possible scenarios, which is left with 64381 safe scenarios. For the identification of Near Misses among these, we resort to their definition as sequences of failure events that indeed keep the system in a safe condition but endangered (i.e., a quasifault system state). To this aim, we introduce a risk-based characterization of these remaining scenarios, calculating their associated risk, at each time instantt, as [23](7)
Risk
t
=
p
t
×
c
t
,where p
t is the probability that at time t the scenario can lead the system into an accidental scenario and c
(
t
) is the consequence that the developing scenario would cause to the system.In this view, we build a functional relationship such thatp increases as N
rl moves further away from the reference level N
ref, in a way that p
=
0 if N
rl is equal to N
ref and p
=
1 if N
rl reaches N
high. Such relationship is given in (8) below, assuming that scenarios whose N
rl
(
t
) approaches N
high are more prone to failure than those with N
rl
(
t
) close to N
ref; that is, (8) “filters out” (i.e., neglects) scenarios whose N
rl
t is close to N
ref and “mines” (i.e., weighs more) scenarios whose N
rl
t is close to N
high:(8)
p
t
=
φ
N
rl
t
-
μ
+
5
σ
σ
=
∫
N
ref
N
rl
t
1
2
π
σ
e
N
rl
t
-
μ
+
5
σ
2
/
2
σ
2
d
N
rl
,where φ is the cumulative probability function of the Gaussian distribution with mean μ
=
N
ref, and standard deviation σ
=
(
N
high
-
N
ref
)
/
5. Figure 9 shows the trend of p
(
t
).Figure 9
Probability functionp
(
t
) for the definition of risk.The consequencec
(
t
) of a scenario increases as N
rl approaches the failure threshold N
high, and c
(
t
) can be calculated at time t as [23] (9)
c
t
=
A
NRL
t
-
μ
+
3
σ
/
NRL
t
-
μ
,where A is the intensity coefficient that accounts for the closeness of N
rl to the thresholds N
h
l, N
v
h, and N
high, and for the exceedance time between the first event of the failure sequence (hereafter called initiating event (IE)) and the time of exceeding the threshold: the shorter this time, the more critical the scenario. Thus, the larger A is, the faster and closer N
rl approaches a threshold; we assume A
=
100 (no consequences) if no threshold is exceeded; A
=
200 (low consequences) if N
rl exceeds N
h
l after at least 2001 (s) from IE or if N
rl exceeds N
v
h after at least 3001 (s) from IE; A
=
300 (medium consequences) if N
rl exceeds N
h
l within 2000 (s) from IE, if N
rl exceeds N
v
h and the elapsed time is in [
1001,3000
] (s), and if N
rl exceeds N
high after at least 2001 (s) from IE; A
=
400 (catastrophic consequences) if N
rl exceeds N
v
h within 1000 (s) from IE or if N
rl exceeds N
high and the elapsed time from IE is in [
1,2000
] (s). A matrix representation of the intensity coefficient is shown in Figure 10.Figure 10
Matrix representation of the intensity coefficientA.By so doing, the available 64381 remaining safe scenarios are fully described at each time instantt
=
1,2
,
…
,
4000 [s] by their values of probability p
(
t
), consequence c
t, and overall risk r
(
t
). An example of the p
(
t
), c
(
t
), and r
(
t
) evolutions for two generic trends of N
rl
(
t
) is shown in Figure 11. More specifically, the N
rl
(
t
) behaviors represented in Figure 11(a) are due to(i)
solid line: the PID controller fails at 100 (s) with magnitude 4 and the safety relief valve fails at 190 (s) with magnitude 2;
(ii)
dashed-dotted line: the safety relief valve fails at 100 (s) with magnitude 1, the communication between the sensor measuringN
rl and the PID controller is interrupted at 136 (s), and the PID controller fails at 3917 (s) with magnitude 5.It is worth analysing the behavior of p
(
t
), c
(
t
), and, thus, r
(
t
) (Figures 11(b), 11(c), and 11(d), resp.): all three abovementioned functions increase as N
rl
(
t
) moves further away from N
ref and decrease as N
rl
(
t
) approaches N
ref. The steps shown in the consequences and risk plots (around 800 [s] for the solid line scenario and around 3500 [s] for the dashed-dotted line scenario) are due to the change of the discrete consequence intensity coefficient A along the scenarios. The solid line scenario is faster than the dashed-dotted line scenario (upper plot) and, thus, the value of the parameters A for the former scenario is 400 (catastrophic consequences, see Figure 10), due to the fact that N
rl
(
t
) exceeds N
v
h within 1000 (s), whereas, A
=
300 (medium consequences, see Figure 10) for the dashed-dotted scenarios, because N
rl
(
t
) exceeds N
v
h within [
1001,3000
] (s). Thus, the solid line scenario is more abrupt in its development towards failure and expected to have more catastrophic consequences, and, thus, more overall risk, than the dashed-dotted scenario, because the time between IE and the exceedance of N
v
h is shorter (i.e, less grace time).Figure 11
Probabilityp
(
t
), consequences c
(
t
), and risk r
(
t
) for two sequences of events.
(a)
(b)
(c)
(d)


### 3.2.1. Features Selection

The identification of the Near Misses is treated as an unsupervised classification problem and addressed by clustering, where (i) the number of clusters is unknown and (ii) the features that enable the best clustering according to the risk-based characteristic profiles ofp
(
t
), c
t, and r
(
t
) of the accidental scenarios are unknown. Unsupervised clustering, thus, entails identifying the number K of clusters in which similar scenarios can be grouped according to similar values of some scenario features. To do this, from the profiles p
(
t
), c
t, and r
(
t
), we extract some statistical indicators as features [36]:(1)
mean value:(10)
μ
n
=
1
T
miss
∑
t
=
1
N
n
t
,
(2)
peak value:(11)
max
⁡
=
max
⁡
t
=
1,2
,
…
,
T
miss
n
t
,
(3)
standard deviation:(12)
σ
n
=
1
T
miss
-
1
∑
t
=
1
T
miss
n
t
-
μ
2
,
(4)
root mean square:(13)
RMS
=
1
T
miss
∑
t
=
1
T
miss
n
t
2
,
(5)
skewness:(14)
SK
=
∑
t
=
1
T
miss
n
t
-
μ
3
T
miss
-
1
,
(6)
kurtosis:(15)
KU
=
∑
t
=
1
T
miss
n
t
-
μ
4
T
miss
-
1
,where n
(
t
) is alternatively equal to p
(
t
), c
t, and r
(
t
) and, thus, the total number of features is equal to 6
×
3
=
18. Among these 18 available features, we search for those that are optimal for clustering the 64381 scenarios in Near Misses and safe scenarios.We resort to a wrapper framework [37, 38], whereby a modified binary differential evolution (MBDE) search engine [33, 39] searches candidate groups of features sets that are fed to a K-means clustering algorithm [40]; eventually, the wrapper evolves so that among these candidate groups, the group retained is that which makes the K-means clustering algorithm perform best (most compact and separate clusters). The idea behind the wrapper approach is shown in Figure 12. During the features search by MBDE, the K-means clustering is run on the N
=
0.80
×
64381
=
51505 (training) safe scenarios with sets of features (F) that are randomly selected by the MBDE algorithm. The optimal number (K) of clusters is also unknown and it is determined by looking at the clustering performance obtained by the K-means with reference to the Calinski-Harabasz (CH) index [41], which accounts for the ratio of the overall between-cluster variance (separation) and the overall within-cluster variance (compactness). The search proceeds iteratively until the CH index is maximised and the number of clusters K is fixed. Then, the results of the wrapper algorithm are evaluated on an independent test set (T
set), that is, the 0.2
×
643281
=
12876 safe scenarios that have been left out during the training phase.Figure 12
Wrapper approach for optimal feature subset selection based on a MBDE optimization algorithm and aK-means classifier.The CH index for a numberK of clusters, k
=
1,2
,
…
,
K is equal to [41](16)
CH
=
S
S
b
S
S
w
×
N
×
F
-
K
K
-
1
,where S
S
b is the overall between-cluster variance:(17)
S
S
b
=
∑
k
=
1
K
n
k
m
k
-
m
2
,and S
S
w is the overall within-cluster variance:(18)
S
S
w
=
∑
k
=
1
K
i
∑
x
∈
k
x
c
-
m
k
2
,where n
k is the number of scenarios x
c assigned to the kth cluster, m
k is the centroid of the kth cluster, that is, the mean of the selected features belonging to the kth cluster, m is the mean of the selected features, and m
k
-
m
2 and x
c
-
m
k
2 are the L
2 norms, that is, Euclidean distances, between the two vectors.The optimal features selection provides as best features the standard deviation ofc
t, the standard deviation of r
t, and the root mean square of r
t; the best performance is obtained with CH
=
9.35
e
+
04 and K
=
5.


### 3.2.2. The Clustering Results

TheK
=
5 obtained clusters of the safe scenarios are shown in Figure 13 with reference to the features of mean risk (μ
risk) and time elapsed from the instant t
risk at which r
t starts to deviate from zero, that is, the time interval during which the system is exposed to risk. The rationale behind this choice is that the larger μ
risk and the longer t
risk, the more dangerous the scenarios. In Figure 13, clusters 3, 4, and 5 (triangles, crosses, and squares, resp.) are well separated, that is, the low level risk scenarios clusters are widened by the adoption of (8) for the quantification of the risk profile r
t. It is possible to distinguish the scenarios having the lowest risk level from the scenarios having low risk level, and, thus, the highest risk scenarios are well separated from the lower risk scenarios. The good performance obtained when (8) is adopted instead of other p
t profiles, for example, linear probability function (p
t
∝
N
rl
t) that would give the same importance to any level N
rl, for the quantification of the risk profile r
t, is due to the fact that (8) “filters out” (i.e., neglects) scenarios whose N
rl
t is close to N
ref and “mines” (i.e., weighs more) scenarios whose N
rl
t is close to N
high: the 332 circles in Figure 13 (listed in Appendix B) can, thus, be considered the Near Misses scenarios, that is, scenarios that incidentally keep the system into safe state, although in endangered and insecure operational conditions.Figure 13
Clustering results.


## 3.2.1. Features Selection

The identification of the Near Misses is treated as an unsupervised classification problem and addressed by clustering, where (i) the number of clusters is unknown and (ii) the features that enable the best clustering according to the risk-based characteristic profiles ofp
(
t
), c
t, and r
(
t
) of the accidental scenarios are unknown. Unsupervised clustering, thus, entails identifying the number K of clusters in which similar scenarios can be grouped according to similar values of some scenario features. To do this, from the profiles p
(
t
), c
t, and r
(
t
), we extract some statistical indicators as features [36]:(1)
mean value:(10)
μ
n
=
1
T
miss
∑
t
=
1
N
n
t
,
(2)
peak value:(11)
max
⁡
=
max
⁡
t
=
1,2
,
…
,
T
miss
n
t
,
(3)
standard deviation:(12)
σ
n
=
1
T
miss
-
1
∑
t
=
1
T
miss
n
t
-
μ
2
,
(4)
root mean square:(13)
RMS
=
1
T
miss
∑
t
=
1
T
miss
n
t
2
,
(5)
skewness:(14)
SK
=
∑
t
=
1
T
miss
n
t
-
μ
3
T
miss
-
1
,
(6)
kurtosis:(15)
KU
=
∑
t
=
1
T
miss
n
t
-
μ
4
T
miss
-
1
,where n
(
t
) is alternatively equal to p
(
t
), c
t, and r
(
t
) and, thus, the total number of features is equal to 6
×
3
=
18. Among these 18 available features, we search for those that are optimal for clustering the 64381 scenarios in Near Misses and safe scenarios.We resort to a wrapper framework [37, 38], whereby a modified binary differential evolution (MBDE) search engine [33, 39] searches candidate groups of features sets that are fed to a K-means clustering algorithm [40]; eventually, the wrapper evolves so that among these candidate groups, the group retained is that which makes the K-means clustering algorithm perform best (most compact and separate clusters). The idea behind the wrapper approach is shown in Figure 12. During the features search by MBDE, the K-means clustering is run on the N
=
0.80
×
64381
=
51505 (training) safe scenarios with sets of features (F) that are randomly selected by the MBDE algorithm. The optimal number (K) of clusters is also unknown and it is determined by looking at the clustering performance obtained by the K-means with reference to the Calinski-Harabasz (CH) index [41], which accounts for the ratio of the overall between-cluster variance (separation) and the overall within-cluster variance (compactness). The search proceeds iteratively until the CH index is maximised and the number of clusters K is fixed. Then, the results of the wrapper algorithm are evaluated on an independent test set (T
set), that is, the 0.2
×
643281
=
12876 safe scenarios that have been left out during the training phase.Figure 12
Wrapper approach for optimal feature subset selection based on a MBDE optimization algorithm and aK-means classifier.The CH index for a numberK of clusters, k
=
1,2
,
…
,
K is equal to [41](16)
CH
=
S
S
b
S
S
w
×
N
×
F
-
K
K
-
1
,where S
S
b is the overall between-cluster variance:(17)
S
S
b
=
∑
k
=
1
K
n
k
m
k
-
m
2
,and S
S
w is the overall within-cluster variance:(18)
S
S
w
=
∑
k
=
1
K
i
∑
x
∈
k
x
c
-
m
k
2
,where n
k is the number of scenarios x
c assigned to the kth cluster, m
k is the centroid of the kth cluster, that is, the mean of the selected features belonging to the kth cluster, m is the mean of the selected features, and m
k
-
m
2 and x
c
-
m
k
2 are the L
2 norms, that is, Euclidean distances, between the two vectors.The optimal features selection provides as best features the standard deviation ofc
t, the standard deviation of r
t, and the root mean square of r
t; the best performance is obtained with CH
=
9.35
e
+
04 and K
=
5.


## 3.2.2. The Clustering Results

TheK
=
5 obtained clusters of the safe scenarios are shown in Figure 13 with reference to the features of mean risk (μ
risk) and time elapsed from the instant t
risk at which r
t starts to deviate from zero, that is, the time interval during which the system is exposed to risk. The rationale behind this choice is that the larger μ
risk and the longer t
risk, the more dangerous the scenarios. In Figure 13, clusters 3, 4, and 5 (triangles, crosses, and squares, resp.) are well separated, that is, the low level risk scenarios clusters are widened by the adoption of (8) for the quantification of the risk profile r
t. It is possible to distinguish the scenarios having the lowest risk level from the scenarios having low risk level, and, thus, the highest risk scenarios are well separated from the lower risk scenarios. The good performance obtained when (8) is adopted instead of other p
t profiles, for example, linear probability function (p
t
∝
N
rl
t) that would give the same importance to any level N
rl, for the quantification of the risk profile r
t, is due to the fact that (8) “filters out” (i.e., neglects) scenarios whose N
rl
t is close to N
ref and “mines” (i.e., weighs more) scenarios whose N
rl
t is close to N
high: the 332 circles in Figure 13 (listed in Appendix B) can, thus, be considered the Near Misses scenarios, that is, scenarios that incidentally keep the system into safe state, although in endangered and insecure operational conditions.Figure 13
Clustering results.


## 3.3. The “Backward” Approach

Once the Near Misses for the SG high level failure mode have been identified by clustering, we can search for similarities among them in terms of their multiple value sequences, that is, order and timing of event occurrences and deterministic process variables values. This “Backward” approach can lead us to finding the minimum conditions, that is, minimumμ
risk and minimum t
risk, that lead the system into a quasifault state. The problem can be framed as a multiobjective optimization problem (MOP) [42] that looks for the set of scenarios x
- to dominate any other scenarios with respect to the fitness function f:(19)
f
x
=
f
1
μ
risk
,
f
2
t
risk
,where f
1 and f
2 are the objectives functions of the defined MOP, that is, minimum μ
risk and t
risk, respectively. The solution of the MOP of (19) is the Pareto set shown in Figure 14, where 12 solutions are plotted (squares lined by continuous line) and listed in Table 3. These scenarios x
- of minimum f
x are expected to cover all failure of Near Misses scenarios cluster.Table 3
List of the Pareto-optimalx
- sequence vectors.
x
-
Steam valve failure time
Steam valve failure magnitude
Steam valve failure order
Safety valve failure time
Safety valve failure magnitude
Safety valve failure order
Sensor-PID failure time
Sensor-PID failure magnitude
Sensor-PID failure order
PID failure time
PID failure magnitude
PID failure order
1
0
0
0
2
3
1
4
1
2
4
3
3
2
3
1
4
1
2
3
1
1
2
1
4
1
3
3
1
4
2
2
3
2
1
1
2
4
2
4
3
1
3
2
2
2
4
1
4
1
4
1
5
3
1
4
3
3
3
2
1
2
1
4
1
6
3
2
4
1
2
2
1
1
3
1
4
1
7
3
2
4
2
3
1
2
1
2
3
3
3
8
4
1
3
3
2
2
4
1
4
2
4
1
9
4
2
3
1
3
2
0
0
0
1
3
1
10
4
2
4
2
1
3
2
1
1
2
4
2
11
4
2
3
4
2
2
0
0
0
4
4
1
12
4
2
4
4
2
3
2
1
1
4
4
2Figure 14
Pareto front for the cluster of Near Misses.The coverage can be verified by, first, identifying the most similar characteristics of the sequence vectors belonging to the Near Misses cluster with the Pareto set scenariosx
-, and, then, by solving a set covering problem (SCP) [43, 44].The most similar characteristics can be computed by coverage vectors (one for each scenario belonging tox
-): this entails calculating the Hamming distance [45] between each sequence vectors in x
- and each one of the other sequence vectors in the Near Misses cluster [46]. The entries of the coverage vector (in our case twelve entries, one for time, magnitude, and order of occurrence of each component failure; see Figure 5) are increased if the Hamming distance between one same entry of the considered scenario belonging to x
- and of the Near Misses vectors is equal to zero, as shown, without loss of generality, in Figure 15 for 1 sequence vector of x
- and only 2 Near Misses vectors.Figure 15
Coverage vector computation by Hamming distance.Table4 lists the 12 coverage vectors, where each entry is the percentage of Near Misses vectors having the same stochastic behavior of the optimal set x
- shown in Table 3. It can be seen that, for each scenario belonging to x
-, columns 8, 11, and 12 (e.g., sensor-PID communication failure magnitude, PID failure magnitude, and PID order of failure, resp.) have the largest values of the coverage vectors: this means that the majority of the sequence vectors of the Near Misses clusters can be well represented by (only) these failures. Furthermore, the analysis of the MVL values of the scenarios belonging to x
- (Table 3) where the largest coverage values of these columns are registered (i.e, 87%, 98.5%, and 85.2% for columns 8, 11, and 12, resp.) highlights that these failures are characterized by the same MVL values that can be summarized as follows:(i)
the failure of the communication between the sensor monitoring theN
rl and the PID controller;
(ii)
the failure of the PID controller with magnitude belonging to[
-
5
,
-
1
]
% of the Q
e value that should be provided at P
o, that is, magnitude equal to 4 in MVL framework, and it is the first accident occurring along the sequence of events in over 85% of the Near Misses scenarios.A SCP can, thus, be solved for verifying that these latest characteristics are the minimum set of stochastic event occurrences and deterministic process variables values of x
- that exhaustively describe the scenarios belonging to the Near Misses cluster: if a Near Miss sequence vector is characterized by (at least) one of the common characteristics, this is covered by the optimal set x
-. In the present application we have verified that all the scenarios belonging to the identified Near Misses cluster are covered by the minimal conditions that lead the system into a quasifault state, that is, the optimal set x
-. In conclusion, the occurrence of one of the common characteristics listed above is sufficient to lead the system into endangered and insecure operational conditions.Table 4
List of coverage vectors for each scenario belonging to the Pareto setx
-.
x
-
Steam valve failure time
Steam valve failure magnitude
Steam valve failure order
Safety valve failure time
Safety valve failure magnitude
Safety valve failure order
Sensor-PID failure time
Sensor-PID failure magnitude
Sensor-PID failure order
PID failure time
PID failure magnitude
PID failure order
1
11.4
11.4
11.4
8.1
30.4
0.6
25.3
87
46.4
1.2
1.5
0.6
2
30.1
44.3
27.1
1.2
19.9
38
19.6
87
46.4
68.1
98.5
85.2
3
30.1
44.3
27.1
8.1
19.9
38
22.3
87
14.2
27.4
98.5
14.2
4
30.1
44.3
42.2
8.1
19.9
19.9
25.3
87
12.7
68.1
98.5
85.2
5
30.1
44.3
27.1
28
30.4
38
22.3
87
46.4
68.1
98.5
85.2
6
30.1
36.1
27.1
1.2
19.9
19.9
19.6
87
13.9
68.1
98.5
85.2
7
30.1
36.1
27.1
8.1
30.4
0.6
22.3
87
46.4
3.3
1.5
0.60
8
52.4
44.3
42.2
28
19.9
19.9
25.3
87
12.7
27.4
98.5
85.2
9
52.4
36.1
42.2
1.2
30.4
19.9
13
13
13
68.1
1.5
85.2
10
52.4
36.1
27.1
8.1
9.04
38
22.3
87
14.2
27.4
98.5
14.2
11
52.4
36.1
42.2
51.2
19.9
19.9
13
13
13
1.2
98.5
85.2
12
52.4
36.1
27.1
51.2
19.9
38
22.3
87
14.2
1.2
98.5
14.2


## 4. Conclusions

In this paper, a risk-based clustering approach for Near Misses identification has been proposed. The approach includes a risk-based feature selection task, where by each safe scenario it is described in terms of probability, consequence, and overall risk. The optimal features set is identified by a wrapper approach based on the combination of a MBDE algorithm withK-means clustering. The characteristics of the Near Misses scenarios are, then, identified solving a multiobjective optimization problem and Hamming distance as a measure of similarity.The application of the approach to a case study of IDPSA of a UTSG has shown the possibility of retrieving relevant information for risk monitoring.



---

*Source: 693891-2015-08-04.xml*