A Distributed Cyber-Physical Framework for Sensor Fault Diagnosis of Marine Internal Combustion Engines

This article proposes a distributed model-based methodology for the diagnosis of faults affecting multiple sensors used for condition monitoring and control of marine internal combustion engines (ICEs). To handle the complexity of the ICE, we consider it as a set of interconnected physical subsystems that constitute the physical layer. For every subsystem, the detection of sensor faults relies on the design of cyber agents, where every agent monitors one subsystem. To handle the heterogeneous dynamics of each subsystem in the fault detection decision-making process, each agent uses differential and algebraic residuals alongside adaptive bounds. For isolation purposes, a combinatorial decision logic is employed, realized in two cyber levels: the local and the global decision logic. The first aims at the recognition of all sensor fault patterns that might have affected the engine based on the local agent fault signatures and certain binary decision matrices. The latter is used to capture the propagation of sensor faults between the different monitoring agents. Simulation results are used to showcase the proposed methodology’s efficiency in tackling the problem and its applicability.

Abstract-This article proposes a distributed model-based methodology for the diagnosis of faults affecting multiple sensors used for condition monitoring and control of marine internal combustion engines (ICEs).To handle the complexity of the ICE, we consider it as a set of interconnected physical subsystems that constitute the physical layer.For every subsystem, the detection of sensor faults relies on the design of cyber agents, where every agent monitors one subsystem.To handle the heterogeneous dynamics of each subsystem in the fault detection decisionmaking process, each agent uses differential and algebraic residuals alongside adaptive bounds.For isolation purposes, a combinatorial decision logic is employed, realized in two cyber levels: the local and the global decision logic.The first aims at the recognition of all sensor fault patterns that might have affected the engine based on the local agent fault signatures and certain binary decision matrices.The latter is used to capture the propagation of sensor faults between the different monitoring agents.Simulation results are used to showcase the proposed methodology's efficiency in tackling the problem and its applicability.
Index Terms-Differential algebraic equations, fault diagnosis, fault resilience metrics, interconnected systems, marine safety, nonlinear systems.

I. INTRODUCTION
M Arine vessels are currently used for the transport of millions of passengers and tons of cargo worldwide.Securing their safety through improved monitoring of the vessel's vital systems and sensors is one of the main goals in the maritime field.Despite this, recent reports from the European Maritime Safety Agency [1] suggest that 37 % of vessel accidents for the years 2014-2020 are attributed to system or equipment failure.In addition, 28% of all vessel casualties for the year 2020 are attributed to loss of propulsion, an event that leaves the vessel ungoverned in the highly uncertain sea environment.As a result, the propulsion system of vessels is considered critical for overall safety [2], [3].In the majority of marine vessels, the propulsion system includes at least one internal combustion engine (ICE).
Another relevant concern is that of resiliency against multiple vulnerabilities (e.g., faults) affecting marine systems [4] such as the marine ICEs.Toward realizing the vision for ship autonomy, the integration of cutting-edge technologies such as novel sensors and more sophisticated monitoring architectures is expected to significantly increase, leading to cyber-enabled ships [5].Although greater operational dependence of vessels on interconnected cyber-physical infrastructures can contribute to assisting decision-making and reducing the effects of human errors [6], it can also lead to easier propagation of faults and, as a result, flawed decisions.The fault diagnosis task, thus, becomes more and more important for vessel-wide resiliency.
Current literature on marine ICEs has addressed the diagnosis of process faults while sensor faults have mostly been overlooked.However, one or more sensors reporting erroneous information could result in masking of process faults, unneeded maintenance, and flawed decision-making by operators [7].A misstep in the onboard decision-making during the voyage will in turn affect the lives onboard the vessel, the transported cargo, and possibly the environment.In addition, due to the larger size of engines, limited operator access, as well as the highly uncertain sea environment, the consequences of a wrong decision regarding faults are in most cases greater in a marine system than a land-based system [8].
The diagnosis of marine ICEs typically takes into account only single sensor values while a marine ICE can incorporate 15-17 sensors [9].Making use of the complete sensor network has been proven useful in effectively supporting decisionmaking onboard marine vessels [10].Moreover, the occurrence of multiple sensor faults has become a significant problem to tackle, due to the large number of sensors distributed in marine systems [11].Although more realistic, this problem has not yet received much attention in the relevant literature.
In the area of sensor fault diagnosis for marine ICEs there has already been some research activity regarding data-driven techniques.In [12], an artificial neural network technique is applied to validate the coherence of sensor measurements, isolate faulty sensors, and recover the lost information from healthy measurements for a six-cylinder engine.However, only single faults are considered and a massive amount of data is needed for training purposes in order to include both healthy and faulty conditions.Perera [13] in their work propose fault detection in two levels regarding sensors used for ship performance analysis (e.g., shaft speed and fuel consumption sensors).The first level's purpose is to detect faults like repeated data points and data points outside selected fixed thresholds.Using fixed thresholds may increase the conservativeness of the approach and careful selection is essential to minimize false alarms.Then, the second level aims to detect in-range sensor faults by using localized models constructed by the data and principal component analysis.The models though only describe three regions of operation for the engine.Tsaganos et al. [14] focus on the diagnosis of faults related to multiple sensors (e.g., shaft speed, pressure, and temperature sensors) of marine diesel engines.To this end, data from an engine simulator are used based on different fault scenarios, and multiple model-free algorithms are tested and compared to these datasets.The results showcase the performance-wise superiority of the AdaBoost Algorithm with a Simple Cart base classifier.The limitation in this case is that the algorithm requires more time for classification.
Regarding model-based methods, in [15], a sensor fault detection methodology is suggested for temperature sensors of selective catalyst reduction systems built upon diesel engines using a suitable temperature model and an extended Kalman filter.In the decision-making process though, the authors choose an arbitrary threshold to perform fault detection, which may lead to missed detection of sensor faults or false alarms if it is not well selected.In [16], a model-based sensor fault diagnosis method is proposed for engine test beds using a multistage geometric analysis of the extracted residuals.However, this method considers static models and single-sensor fault occurrence.Finally, Stoumpos and Theotokatos [17] discuss a unified digital system for diagnosis and health management for dual-fuel marine engines by combining a thermodynamic model with data-driven methods to diagnose sensor faults and compensate for their effects.More precisely, a feed-forward neural network approach is used, enhanced by expert knowledge to improve isolation while considering two case studies where faults occur at two engine sensors simultaneously.This work lacks the performance analysis or experimental validation of the proposed approach.
Quantifying resilience for cyber-physical systems through the use of well-defined metrics has been an active topic in recent years and is necessary for the certification of these systems.Resilient control of complex cyber-physical systems such as marine ICEs requires the timely and effective capture of sensor fault events so that their effects can be properly identified and mitigated.In the fault diagnosis community, metrics to characterize the performance of the method are mostly related to the detectability such as the minimum detectable fault magnitude, the missed detection rate, and the detection delay [18], [19].The performance of the ability to isolate faults has attracted less attention, especially considering multiple sensor faults and interconnected differential-algebraic systems.
The goal of this article is the design and performance analysis of a distributed model-based cyber-physical framework for the diagnosis of faults affecting multiple sensors used for condition monitoring and control purposes of marine ICEs.A mean value first principle (MVFP) engine model [20], [21] described by differential-algebraic equations (DAEs) is considered to represent the physical layer.Sensor fault detection is conducted by designed cyber agents using both state and algebraic estimators alongside sensor measurements to form the various residuals, which are then compared to designed adaptive thresholds.As for the isolation of sensor faults, a combinatorial decision logic [22] is employed consisting of two cyber levels.At the first local level, sensor fault detection and isolation (SFDI) agents are designed to monitor specific sensor sets and diagnose faults in those sets.Then, at a higher level, a global agent is designed to isolate multiple sensor faults that may be propagated between the local monitoring agents.The design of the various agents is assessed using performance analysis tools [23], extended, however, to the needs of the differential-algebraic nature of the system.This assessment of performance can pave the way for certification of the monitoring system design and, in turn, provide resilience assertions.
Compared to the authors' previous work [21], the present article also includes the performance analysis of the designed monitoring scheme, proposing new metrics regarding the detectability and ability to isolate sensor faults which are supported by extensive simulation results.Moreover, we elaborated on the design of the algebraic estimators and the adaptive thresholds where set inversion via interval analysis (SIVIA) [24] has been used.This method allows us to decrease the conservativeness in the algebraic adaptive threshold design and further improve the diagnosis process.
The main contributions of this article are listed.
1) The development of the distributed SFDI architecture for marine ICEs (see Section III).2) In addition, the use of adaptive thresholds in Section IV for detection purposes allows for the reduction of conservativeness in decision-making, excluding false alarms.
3) The fact that the MVFP model, provided in the Appendix, can be easily reconfigured in both parameters and dynamics to host more subsystems associated with the ICE allows the generalization of the results and enhances the applicability of the method.4) This article considers a system described by nonlinear DAEs, with the relevant formulation given in Section II, while most papers in SFDI literature deal with systems described by ordinary differential equations (ODEs) [25], [26], [27], [28].5) Furthermore, the present article also includes performance metrics regarding the propagation of sensor faults between different monitoring modules, in Section VI-A, in addition to the local ones found extensively in literature [18], [29].6) Another contribution lies with the performance assessment of the ability to isolate sensor faults in Section VI-B, a property which has been scarcely discussed in relevant works in literature.Simulation results using MATLAB are shown in Section VII, followed by some concluding remarks in Section VIII.

II. PROBLEM FORMULATION
Marine ICEs are complex systems incorporating components characterized by heterogeneous dynamics and inherent interconnections.In Fig. 1, a simplified representation of a marine ICE is shown, where the different parts are grouped in four distinct subsystems and a total of ten sensors are deployed for condition monitoring.In most cases, the marine engine is controlled by a regular PI controller [20].Subsystem 1 incorporates the fuel pump, used to supply fuel to the engine's cylinders and its output fuel injection sensor.Subsystem 2 is the thermomechanical process which refers to both the thermal and mechanical processes occurring inside the engine.As fuel oil, air, and exhaust gases are all present in this block, it is characterized by heterogeneous dynamics.One pressure, one temperature, and one torque sensor are considered for this subsystem.Subsystem 3 consists of the exhaust manifold and the turbine, both handling the exhaust gas, having therefore similar dynamics.As for the sensors, a pressure sensor after the exhaust manifold and two temperature sensors before and after the turbine are considered.Finally, subsystem 4 includes the compressor and the intercooler of the engine, both dealing with the air supply.A pressure sensor is considered after the compressor as well as two temperature sensors before and after the intercooler.This work aims to diagnose faults affecting multiple sensors in the marine fuel engine system.
Given the heterogeneous dynamics and interconnections of the subsystems in marine ICEs, every physical system (I ) , I = {1, . . ., N } (N = 4) is described by a set of DAEs [29] where x (I ) ∈ R n I −r I is the state variable vector, z (I ) ∈ R r I is the algebraic variable vector, χ (I )  ∈ R k I are the interconnection variables from the neighboring subsystems, u (I )  ∈ R l I is the control input vector, γ (I )  : R n I −r I × R l I → R n I −r I represents the known nonlinear system dynamics, and h (I )  : represents the known interconnection dynamics with the neighboring subsystems, ξ (I ) : R n I × R k I × R l I → R n I −r I is a smooth vector field.The term A (I ) x (I ) represents the linear part of the system's (I ) dynamics, where A (I ) ∈ R (n I −r I )×(n I −r I ) is assumed known.Each system incorporates a set of sensors S (I ) = n I j=1 S (I ) { j} described as S (I ) : where y (I ) , where f (I ) j (t), j ∈ J (I ) = {1, . . ., n I } denotes the change in the output due to a fault in the jth sensor.Permanent abrupt faults can be modeled as follows [19]: where T (I ) f j is the time instant of occurrence of the jth fault and φ(I) j is its fault signature.The objective of this article is to design a methodology for the detection and isolation of multiple, permanent sensor faults for nonlinear DAE interconnected subsystems described by ( 1) and ( 2), subject to the following assumptions.
Assumption 1: The measurement noise of each sensor is unknown but uniformly bounded, meaning Assumption 2: The nonlinear vector fields γ (I ) , h (I ) are locally Lipschitz in x ∈ X , z ∈ Z for all u ∈ U and t ⩾ 0 with Lipschitz constants λ γ I , λ h I , respectively.
The differential-algebraic model for marine ICEs, such as the one shown in Fig. 1, can be found in the Appendix.Note that hereby the time dependency may be dropped for simplicity of mathematical expressions.

III. DISTRIBUTED SENSOR FAULT DIAGNOSIS ARCHITECTURE
In this section, the architecture of the proposed sensor fault diagnosis methodology is described.As can be seen from Fig. 2, for each subsystem (I ) , I = 1, . . ., N , a monitoring agent M (I ) is designed to monitor the health of the sensor set S (I ) belonging to the specific subsystem.Internally, each agent is composed by q I monitoring modules denoted as M (I,q) , q = 1, . . ., q I (q 1 = 1, q 2 = q 3 = q 4 = 3).The qth local sensor subset, denoted as S (I,q) ⊆ S (I ) , is assigned to be monitored by the module M (I,q) and is expressed as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.where y (I,q) , d (I,q) , f (I,q)  ∈ R m (I,q) , and m (I,q) ≤ n I denotes the cardinality of the sensor subset S (I,q) .Each monitoring module M (I,q) , q = 1, . . ., q I internally compares certain residuals to adaptive thresholds, with the design of both being discussed in Section IV.Communication of sensor values, denoted as y (I,q) χ , between the agents is essential in the calculation of the various adaptive thresholds.Thus, a distributed monitoring architecture is implemented.Given the differential-algebraic formulation (1) of systems considered in this article, a mixed differential-algebraic diagnosis scheme is proposed.
The isolation of sensor faults is then realized in two steps.First, a local decision logic is used to indicate the presence of one or more faults, locally affecting the sensors S (I ) monitored by the agents M (I ) .A global agent is also implemented to enhance the decision process by providing information on whether sensor faults have been propagated from or to the agent M (I ) from the neighboring agents.The isolation process applies a combinatorial decision logic and diagnostic reasoning between the local and global levels and is described in Section V.

A. Residual Generation
The residual vector of the module M (I,q) is defined by where x(I,q) denotes the estimation of the differential state x (I ) involved with the module M (I,q) and is calculated using a standard nonlinear Luenberger estimator as follows: x(I,q) (t) + γ (I ) ( x(I,q) (t), y (I,q) z (t), u (I ) (t)) + h (I ) ( x(I,q) (t), y (I,q) z (t), y (I,q) χ (t), u (I ) (t)) + L (I,q) (y (I,q) x (t) − x(I,q) (t)) where L (I,q) ∈ R (n I −r I )×(n I −r I ) is chosen such that the matrix − L (I,q) is Hurwitz.Subtracting (1) from ( 6) and substituting from (2) yields where can then also be expressed as ϵ

B. Computation of Adaptive Thresholds
The design of adaptive thresholds takes into account the need for them to bind the respective residuals under healthy sensor conditions.Mathematically, the aforementioned design principle can be expressed as Under the Assumptions 1 and 2 and after some mathematical manipulations of ( 7), the adaptive thresholds for Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The thresholds corresponding to the algebraic residuals ϵ (I,q) y z j , ε(I,q) y z j shown in (9) are computed using SIVIA [24].This numerical method benefits from the property of monotonic convergence.Hereafter, the notation [•] will be used to denote the interval of potential values of the variable (•).Under Assumption 1 and using (1b) with x ∈ X (I,q) , z ∈ Z (I,q) , χ (I,q) ∈ X (I,q) , and u (I ) ∈ U (I,q) , the following algebraic estimator can be constructed: = ξ (I,q) (X (I,q) , Z (I,q) , X (I,q) , U (I ) ) (11) where ) , ū(I) ] are known intervals and Z (I,q) is unknown.According to (1b), the inversion set [24] is defined as Y (I,q) = {0}.The target of the numerical method is to estimate the unknown interval box [z (I,q)  ] = Z (I,q) so that (I,q)  ⊆ Y (I,q) .SIVIA is used to dynamically update the unknown interval estimation, starting from an initial estimation Z (I,q) 0 , as can be seen in Fig. 3.The high nonlinearity of the ICE's systems is amended by using forward-backward propagation contractors integrated with SIVIA [24].Then, based on (2) and ( 5) and some mathematical manipulations, the residual interval under healthy sensor conditions can be calculated as follows: Thus, the adaptive thresholds ϵ (I,q) y z j , ε(I,q) y z j are computed as ϵ (I,q) y z j , ε(I,q) y z j = min[ϵ (I,q) y z j ], max[ϵ (I,q) y z j ] .

C. Detection Logic
Sensor fault detection in S (I,q) by the monitoring modules M (I,q) occurs by comparing the previously defined residuals to the designed adaptive thresholds, based on a set of analytical redundancy relations (ARRs).The jth ARR can be defined as for the monitoring modules using the residual expression ϵ (I,q) y x defined in ( 5) and the threshold expression of (10a).Otherwise, the jth ARR is defined as follows: The set of ARRs based on which the module decides on the presence of local sensor faults is defined as , where J (I,q) is an index set, defined as J (I,q) = { j : S (I ) { j} ∈ S (I,q) }.The first time instant that ( 14) or ( 15) is invalid for at least one j ∈ J (I,q) signifies the time instant of fault detection T (I,q) D j by the local SFDI module M (I,q) , defined as y z j (t), ε(I,q) y z j (t)]} accordingly.Until this instant, the local sensing subsystem S (I,q) is considered nonfaulty meaning that either no fault exists or that faults exist but remain undetected.
The output of M (I,q) is denoted by D (I,q) and in the case of permanent sensor faults, which can be defined as with : j ∈ J (I,q) }.

A. Local Decision Logic
From the detection logic step, a binary decision vector = [D (I,1) , . . ., D (I,q I ) ] can be obtained for the monitoring agent M (I ) and compared to the columns of a binary fault signature matrix (FSM) F (I ) , consisting of N I rows and N C I + 2 columns, where N C I = 2 n I − 1.The design of this matrix will be described in the simulation results section for the easiness of the analysis.While D (I ) (t) = 0 N I , the diagnosis set D (I )  s is empty.In addition, if D (I,q) = F (I ) qi ∀q ∈ 1, . . ., N I , then the observed pattern D (I ) (t) is said to be consistent with the theoretical pattern F (I ) i and the diagnosis set is defined as D (I )  s (t) = {F (I ) ci : i ∈ I (I ) D (t)} where I (I ) D (t) is the consistency index set defined as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

B. Global Decision Logic
Together with the local diagnosis set D (I )  s , the agent M (I ) , I ∈ {1, . . ., N } also provides a decision on the propagation of sensor faults from the interconnected subsystems, denoted as D (I )  χ (t) with where f (I ) p ∈ R n * I , n * I ⩽ n I , collectively amounts for the sensor faults that are propagated from the agent M (I ) to its neighboring agents due to the exchange of sensor information and f (I ) χ corresponds to the sensor faults propagated to the agent from the neighboring agents.The global decision logic serves to isolate sensor faults propagated through the interconnections between the monitoring agents.As shown in Fig. 2, a global agent G collects the decisions on the propagation of sensor faults from the N local agents D χ (t) = [D (1)  χ (t), . . ., D (N ) χ (t)] and compares them with the columns of a global binary sensor FSM F χ consisting of N rows and N C = 2 p − 1 columns ( p ⩽ N I =1 { p I }, and p I is the length of f (I ) χ ).A (*) is used in F χ instead of 1 in case the sensor fault is propagated to the agent M (I ) from the other agents M (J ) , J ∈ {1, . . ., N }, J ̸ = I to indicate that a set of ARRs is less sensitive to this fault.The attenuated sensitivity may be due to the propagation effects of sensor faults or the algebraic nature of ARRs.The isolation of sensor faults in S (I,q) requires the combination of the local and global decision logic levels.More precisely, the output D χ s (t) of the global agent is used to update the diagnosis sets D (I ) s (t) of the local monitoring agents M (I ) .The occurrence of f (I )  χ and its combinations can be excluded from the local diagnosis sets, if Moreover, the intersection of the updated sets is considered the diagnosis results.The formal mathematical expression of the resulting global diagnosis set is the following: The isolated faulty sensors set is then denoted as and is a superset of the actual faulty sensors set S F 0 .

VI. FAULT RESILIENCE QUANTIFICATION METRICS
The performance analysis of the proposed distributed monitoring scheme in this article concerns the detectability and ability to isolate sensor faults.This section provides the definitions of the key performance indicators (KPIs) based on which the fault resilience of the proposed monitoring scheme is assessed.

A. Distributed Sensor Fault Detectability
Having designed the adaptive thresholds according to (10a) and ( 13), the occurrence of false alarms is excluded.Thus, the KPIs of interest for sensor fault detectability include the following.
1) Minimum Detectable Sensor Fault Magnitude: The minimum detectable sensor fault magnitude for each of the local monitoring modules M (I,q) can be expressed as 2) Missed Detection Rate: In comparison to existing literature [18], [30], the missed detection rate metric is assessed both locally (local monitoring modules) as well as globally (interconnected monitoring modules).According to [18], the missed detection rate (MDR) for each of the local monitoring agents M (I,q) , I = 1, . . ., N I , q = 1, . . ., q I is expressed as Following a similar rationale, the MDR concerning the propagation of sensor faults will be defined as follows: 3) Detection Delay: The detection delay for each of monitoring modules M (I,q) is expressed as Thus, taking into consideration the aforementioned KPIs, good monitoring performance regarding distributed architectures would be associated with low values for all MDF (I,q) , MDR (I,q) , MDR, and DL (I,q) .

B. Ability to Isolate Sensor Faults
Considering the existing literature on fault diagnosis, no clear metrics on the ability to isolate sensor faults have been proposed so far.In this article, we define the uncertainty U G and the exoneration efficiency E G of the global diagnosis set as follows: where |•| denotes the cardinality of the set (•).The isolation (P U G ) and exoneration (P E G ) performance of the sensor fault isolation process can be then quantified as follows: where T is the running time of the system and T D0 = min{T (I,q) D : I = 1, . . ., N , q = 1, . . ., q I } is the time the first sensor fault gets detected.

VII. SIMULATION RESULTS
In this section, the distributed sensor fault diagnosis architecture described in Sections III-V is applied on the marine ICE model, described in the Appendix.The data needed for the model and the controller are extracted from [20] while the noise bound of each sensor is assumed to be equal to 5 % of the amplitude of the sensor value.
Tables I-IV represent the local and global sensor fault signature matrices.In Table I, the three first columns provide the theoretically expected local sensor fault patterns while the rest of the columns show the propagated sensor fault patterns.For brevity purposes, only the single sensor fault columns are shown while this article deals with multiple sensor faults.To obtain the multiple fault columns, a logical conjunction between the respective single fault columns needs to be applied.As for the rows of this matrix, each of them corresponds to a local module M (2,q) making use of the set of ARRs E (2,q) , q = {1, 2, 3}.A similar approach is also followed for the creation of Tables II and III.The reason for having two modules M (3,1) , M (3,2) monitoring the same set of sensors (S (3,1)  = S (3,2) ) is to improve the ability to isolate sensor faults.For instance, a local decision vector    occurrence of sensor faults f (3)  1 , f (3) 2 .In the case that there were two modules, for example, M (3,1) and M (3,3) , then an observed pattern D (3)  = [D (3,1) From the simulation results shown in Fig. 4, the diagnosis process proceeds as follows.For t < 20 s, the diagnosis set is empty (D = {}), as all detection decisions are zero.For 20 ⩽ t < 45.06 s, the second monitoring agent M (2)  outputs the decision vector D (2) (t) = [0 0 1] ⊤ as can be seen in Fig. 4(b)-(d).Since D (1) 1 , f (3) 1 , and f (4)   1 are excluded from the diagnosis sets.Comparing D (2) with Table I, the resulting local diagnosis set is 4 } .The global decision vector for this time interval is D χ (t) = [0 1 0 0] ⊤ and when compared to Table IV Then, for 45.06 ⩽ t < 69.9 s, the decision vector for the monitoring agent M (2) remains the same, thus resulting in the same local diagnosis set D (2)  s (t).The rest of the local agent decision vectors are D (1)  = 0, D (4)  = [0 0 0] ⊤ , therefore excluding the sensor faults f (1)  1 , f 1 from the diagnosis sets.However, the decision vector for agent M (3) is now equal to D (3)  = [1 0 0] ⊤ , also seen in Fig. 4(e)-(g) and comparing Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.  1) , (c) M (2,2) , (d) M (2,3) , (e) M (3,1) , (f) M (3,2) , (g) M (3,3) , (h) M (4,1) , (i) M (4,2) , and (j) M (4,3) .
that to Table II, the resulting local diagnosis set is 3 } , • • • .The decision vector on the propagation of sensor faults is now D χ (t) = [0 1 1 0] ⊤ and comparing with Table IV, we obtain the diagnosis set on fault propagation 3 } .For 69.9 ⩽ t < 70.01 s, only the decision vector of agent M (3) changes to D (3)  = [1 0 1] ⊤ as can be seen in Fig. 4(e)-(g).Following the same approach, the resulting global diagnosis set changes to 3 } .Finally, for t ≥ 70.01 sec, the local decision vector of agent M (4) is D (4)  = [0 1 0] ⊤ and when compared to Table III, results in the local diagnosis set D (4) 3 } .Based on the above analysis, the proposed methodology managed to isolate sensor faults in S (4)  {2} and combinations of the sensors S (2)  {2}, S (2)  {3}, S {1}, S (3)  {2}, and S (3)  {3}.It is also useful to note that the resulting diagnosis set managed to include the sensors with simulated faults in our scenario.
To investigate the effectiveness of the proposed distributed SFDI scheme, the KPIs mentioned in Section VI were calculated.Uniform noise bands for all sensors were chosen to provide a common point of reference for the analysis.Regarding the minimum detectable sensor fault magnitudes, defined in (19), the respective results for a noise bound of 5% are given in Table V.The algebraic modules M (2,1) , M (2,2) , M (2,3) , M (3,2) , M (3,3) , M (4,2) , and M (4,3) exhibit an MDF (I ) j value according to the complexity of the model equations and the number of interconnection variables χ (I )  ∈ R k I .For instance, the system dynamics of ( 30), (33), or (49b), used by the modules M (2,3) and M (4,3) , respectively, are more complex and include more interconnections than the dynamics of (49a) used by the module M (4,2) .As a result, the MDF (I ) j indicator has a substantially higher value for those modules.Regarding the rest of the modules that use is comparable to the noise level as can be seen in Table V.The results regarding the missed detection rates and detection delays are given for two fault cases affecting the following sensors, namely 1) the fuel mass sensor f (1)   1 ̸ = 0 and 2) the engine torque sensor f (2)   3 ̸ = 0.This particular choice of sensors was made taking into consideration the criticality of these two sensors for the control of marine ICEs.Extensive simulation was used to obtain the data points for ten levels of sensor noise affecting all system sensors (1%-10%) [30].At each noise level, the simulation is run 100 times, with varying seeds.The simulations were carried out using TU Delft's Blue Supercomputer [31].The simulated faults in each of the two different scenarios occur at the time instant T (1)  f 1 = T (2) f 3 = 20s and are again considered permanent, abrupt and offset.As for the magnitudes of the sensor faults, in scenario 1) we consider a fault magnitude of φ(1) 1 = 7%x (1)  nom while in the scenario and 2) a sensor fault magnitude φ(2) 3 = 65%z (2)  3,nom .Moreover, as both MDR (I,q) and MDR are defined as probabilities in expressions (20) and (21), aside from the mean values found by simulation for these indicators, the respective confidence intervals for the mean value are also depicted in Fig. 5 using a shaded area representation (95% confidence level).As illustrated in Fig. 5(a) and (b), both MDR and MDR (I,q)  are increasing or remain constant as the variance of the sensor noise is increasing as well.Moreover, MDR (I,q)  ≤ MDR, which indicates the persistency of the local module design to detect the sensor fault.In addition, Fig. 5(b) suggests that between a noise level of 6%-7%, although the local module misses the fault (MDR (2,3)  = 100%), the interconnected agents manage to still detect the fault effects through propagation.This strength of the distributed architecture is particularly useful for the diagnosis process.As for the detection delay, the values of the indicator DL (I,q) , defined in (22), are given for the same fault scenarios in Table VI both for the local monitoring modules and the interconnected modules.The local module detection delays are relatively small, indicating a solid design of the residual and adaptive thresholds.Again, as the noise level increases so does the detection delay as can be seen in both scenarios.
Regarding sensor fault isolation, the initial scenario with three sensor faults described at the beginning of this section was used.Using TU Delft's Blue supercomputer, the simulation was run 100 times with a noise level of 5% affecting all system sensors and a varying seed.Based on the acquired results, implementing the metrics defined in ( 25) and ( 26) yields P U G = (51.16± 0.81)% and P E G = (57.83± 2.93)% (95% confidence level) suggesting a fair isolation performance for the three consecutive faults scenario.

VIII. CONCLUSION
In this article, we illustrated a distributed cyber-physical framework for isolating sensor faults in marine ICEs.The goal of the proposed methodology was the isolation of sensor faults affecting multiple sensors of the engine using the information exchanged between its different subsystems.The core of the diagnosis approach consists of two cyber layers: one based on a bank of local monitoring agents monitoring specific sensor sets and a global decision logic layer that decides on the propagation of faults between the different subsystems.Each monitoring agent was composed of one or more monitoring modules.Fault detection was carried out by the different modules using differential or algebraic residuals and comparing them to designed adaptive bounds.Sensor faults were isolated by comparing the decisions of the various monitoring agents to combinations of sensor faults and by applying diagnostic reasoning.Future work will focus on expanding the current cyber layer to account for other types of vulnerabilities such as process faults and cyberattacks.

APPENDIX MARINE ICE MODELING
To describe the operation of the ICE, the physical MVFP model described in [20] and [21] is used.As shown in Fig. 1, the MVFP system model is decomposed in four nonlinear DAE interconnected subsystems, formulated using the general formulation shown in ( 1) and (2).
A. Fuel Pump ( (1) ) Subsystem 1 is expressed as where x (1) (t) ∈ R is the amount of fuel injected per cylinder per engine cycle in [kg], x (1)  nom ∈ R signifies the same quantity under nominal engine conditions, u(t) ∈ R is the fuel injection setting in [%], and τ X = (1/(4n nom f e )) is the fuel injection time delay in [sec].The nominal fuel injection amount x (1)   nom is expressed as x (1)  nom = SFC nom P nom f e k e i e n nom f e (28) where n nom f e is the nominal rotational engine speed in [RPS], SFC nom is the nominal fuel consumption of the engine in [kg/Wh], P nom f e denotes the nominal power output of the engine in W, i e is the number of engine cylinders, and k e denotes the number of crank revolutions per engine cycle (k e = 1 for a two-stroke engine and k e = 2 for a four-stroke engine).The output of the fuel injection sensor y (1)  ∈ R is described by B. Thermomechanical Process ( (2) ) This subsystem has three algebraic variables, namely the pressure (z (2)  1 (t)) in [Pa] and the temperature (z (2) 2 (t)) in [K] inside the engine's cylinders and the engine's shaft torque (z (2)  3 (t)) in [Nm].The mathematical representation of the system is 1 − ξ (2)  z1 (x (1) , x (4) , z (4) 1 ) z (2)  2 − ξ (2)  z2 (x (1) , x (4) , z (4) 1 ) z (2)  3 − ξ (2)  z3 (x (1) , z (2) 3 , x (4) , z   = ξ (2) (x (1) , z (2) , x (4) , z (4) 1 ) where the functions ξ (2)  z1 , ξ (2)  z2 , ξ (2)  z3 ∈ R can be modeled using the Seilinger thermodynamic cycle as follows [20]: ξ (2)  z3 = v 1 i e x (4)   2π k where n f e = ((2π/c)z ( 2) is the nominal constant volume portion, X grad cv is the gradient of the constant volume portion, X nom ct is the nominal constant temperature portion, η is the thermal efficiency incorporating both the combustion and heat release processes, h L is the lower heating value of fuel at ISO conditions in [J/kg], R a is the gas constant of air in [J/kgK], v 1 is the cylinder volume at start of compression in [m 3  ], r c is the effective compression ratio determined by the inlet valve timing, κ a is the specific heat ratio of the air, r eo is the ratio of the volume at Seiliger point 6, n exp is the polytropic exponent for expansion, c p,a is the specific heat at constant pressure for the scavenge air in [J/kgK], c v,a is the specific heat at constant volume for the scavenge air in [J/kgK], Q nom loss denotes the nominal mechanical losses of the engine in [Nm], Q grad loss denotes the gradient of mechanical losses of the engine in [Nm], and c is a constant.

x∈
R n I −r I denotes the sensor values of state variables, y (I ) z ∈ R r I denotes the sensor values of the algebraic variables, d (I ) x ∈ R n I −r I , d (I ) z ∈ R r I are the measurement noise vectors, and f (I ) x ∈ R n I −r I , f (I ) z ∈ R r I are the sensor fault vectors.Each fault vector is given by

Fig. 2 .
Fig. 2. Distributed SFDI architecture for marine ICEs.The distributed aspect is due to the communication between monitoring agents, imitating the interconnections in the physical layer (see the Appendix).The interconnection variables between the agents, denoted as y (I,q) χ , are shown with red, dotted lines.

Fig. 5 .
Fig. 5. Missed detection rates with respect to increasing sensor noise levels d(I)j for all sensors.Results for two single sensor fault scenarios are shown (blue/continuous line: MDR (I,q) , magenta/ dash-dotted line: MDR).Each point on the curves corresponds to the times that the corresponding diagnosis agent failed to detect the presence of the sensor fault out of the 100 simulations obtained for each sensor noise variance d (I ) j (t).The confidence intervals are also shown as a shade (95% confidence level).(a) f(1)  1 ̸ = 0. (b) f(2)  3 ̸ = 0.

ψ 4 = 1 +η
η tur (x(4) )( tur − 1tur (x(4) ) = a tur + b tur x(4)  + c tur (x(4) )2  (43) τ p d is the time delay for filling the exhaust receiver in [sec], p ex is the pressure after the turbocharger in [Pa] assumed equal to the atmospheric pressure, a Z is the Zinner turbine area decrease factor assumed 1 for a constant pressure turbocharger, A eff is the turbine's effective area in [m 2 (t) z (3) (t) ⊤ + d (3) (t) + f (3) (t).(44) D. Air Path ( (4) ) z1 = T c − ϵ inl (T inl − T c ) (49a) ξ (4) z2 = T amb + χ g η tur (δ f (t) + η com )x (3) z (3) 2 − z (3) 1 (49b) τ T C is the compressor time delay in [sec], p amb is the ambient pressure in [Pa], T amb is the ambient temperature in [K],a η , b η , c η are the polynomial coefficients of the turbocharger for estimating its efficiency, η com is the mechanical efficiency of the compressor which can be considered constant, T c is the charge air temperature after the intercooler in [K], ϵ inl is the parasitic effectiveness of the heat exchange between the inlet duct and the air, and T inl is the temperature of the inlet duct that heats the inducted air in [K].

TABLE I PART
OF THE SENSOR FSM OF M(2)(SINGLE FAULTS)

TABLE II PART
OF THE SENSOR FSM OF M(3)(SINGLE FAULTS)

TABLE III PART
OF THE SENSOR FSM OF M (4) (SINGLE FAULTS)

TABLE IV PART
OF THE GLOBAL SENSOR FSM (SINGLE FAULTS) , results in the diagnosis set on fault propagation D χ s

TABLE V MINIMUM
DETECTABLE SENSOR FAULT MAGNITUDES BY RESPECTIVE LOCAL MONITORING MODULES (NOISE BOUND 5%)