Sensor data play an instrumental role in system health prognostics, degradation, fault detection, maintenance and control [1][2][3][4][5][6][7][8]. Sensor monitoring signals, arising from a complex engineering system, are not only statistically-correlated but also physically-dependent through unequivocal system working principles, governing physics, system configuration, etc. Very often, true system health states are not directly observable and need to be estimated from sensor monitoring signals. Known system physics imposes fundamental constraints and regularizations on how sensor data can be used to estimate hidden system health states. When a solid connection is established between sensor monitoring data and system hidden health states through system working principles, the interpretability of data-driven system prognostics can be significantly strengthened. The objective of this paper is to propose a physics-regularized data-driven approach for the health prognostics of complex engineered systems with dependent health states, using sensor monitoring data and system working principles. to cool the air in server rooms of a DC. Used water is re-circulated back to the chiller to be cooled again. Fig. 1 shows the schematic of the reciprocating chiller-one of the most common commercial chillers. The essential components of a chiller system include the compressor, expansion device, condenser and evaporator. In a refrigeration cycle, low pressure and low temperature vapor is fed to the compressor. The compressor increases both the pressure and temperature of the vapor. High pressure and high temperature vapor is then passed to the condenser and is cooled by giving up its latent heat. As a result, the vapor condenses back to its liquid form. The high pressure liquid from the condenser is then expanded through an expansion valve. At this point the refrigerant is at a low pressure and is mostly liquid with low boiling point. When the low-pressure liquid refrigerant enters the evaporator coils, it boils and absorbs the latent heat of evaporation from the surrounding air. The vapor at low pressure and low temperature then passes to the compressor and the whole refrigeration cycle repeats itself.
Due to the importance of DC cooling systems, critical operating parameters of chillers in DC are closely monitored by sensors. For example, Fig. 2 shows the observed daily Coefficient of Performance (COP), condenser coolant inlet temperature 𝑇 𝑐 , evaporator coolant outlet temperature 𝑇 𝑒 , and cooling capacity 𝑄 𝑒 over a 57-day study period (the data are provided by a major DC operator). Here, COP is an overall indicator of a chiller's energy efficiency, defined as the ratio, COP = 𝑄 𝑒 ⋅ 𝑃 -1 , between the cooling capacity 𝑄 𝑒 (i.e., the rate of heat withdrawn from the data center server room) and the power input 𝑃 (i.e., energy consumption rate of a cooling system) [10]. It is seen from Fig. 2 that the daily COP gradually degrades over the 57-day period (a higher COP equates to higher energy efficiency and lower operating cost). The daily condenser coolant inlet temperature 𝑇 𝑐 varies between 297.5𝐾 to 300𝐾. This parameter is often affected by not only the temperature of the chilled water produced by the chiller, but also other external factors such as room temperature and computing load of the servers in the computer room. The variation of the daily evaporator coolant outlet temperature 𝑇 𝑒 is extremely small (less than 0.5𝐾) because 𝑇 𝑒 is the temperature of the cooled water that the cooling system is supposed to supply. The cooling capacity 𝑄 𝑒 varies between 200𝐾𝑊 and 250𝐾𝑊 . Hence, there exist both practical need and theoretical interest to answer a fundamental question: how can the temporal dynamics (e.g., degradation) of hidden system health states be estimated from multiple sensor monitoring data?
Addressing the question above is confronted with multiple challenges (which apply to not only the motivating application above, but also many other health prognostics problems for engineering systems):
• System physics imposes fundamental modeling constraints and regularization, which need to be integrated into data-driven system health prognostics. In the motivating example, the governing physics between COP and other critical operating conditions can be described by the first law of thermodynamics [11]:
where the condenser coolant inlet temperature 𝑇 𝑐 gives the temperature of the water cycled back to the cooling system from the computer room, the evaporator coolant outlet temperature 𝑇 𝑒 gives the temperature of the cooled water produced by the chiller, and the cooling capacity 𝑄 𝑒 measures the rate of heat withdrawn from the computer room which can be calculated from other parameters such as the measured water flow rate, pipe diameter, etc. The three parameters of the governing physics, (𝛾 1 , 𝛾 2 , 𝛾 3 ), characterize the internal irreversibilities states of a particular chiller. Hence, the thermodynamics model (1), as a fundamental system working principle, determines the critical relationship between COP (i.e., energy efficiency) and multiple observed operating parameters. Such a relationship can hardly be faithfully recovered by black-box approaches driven by the statistical correlation among sensor signals, calling for physics-informed statistical health prognostics approaches.
• The true health states of a complex engineering system are usually hidden and not directly observed by sensors. Having sensor measurements does not automatically imply that one has measured the right variables. It is often necessary to properly define and estimate hidden system health states from observable sensor signals which are dependent on each other due to some fundamental system physics. In the motivating example, it is meaningful to treat the internal irreversibility states (𝛾 1 , 𝛾 2 , 𝛾 3 ) in (1), or some functions of these parameters as system health state variables, and estimate the defined health state variables from multiple sensor signals by invoking the system physics (1). In this case, each system state variable possesses an interpretable physical meaning.
In fact, individual sensor signals are rarely ideal measures of the health state of an engineering system. Each sensor monitors a single parameter (dimension) which only reflects the ''local behavior'' of the chiller. The energy efficiency (COP) of a chiller depends on various external environmental factors and internal system health state. The change of working load and outdoor temperature may cause the drop of the observed COP, which does not necessarily imply that the system health state has degraded. Hence, instead of focusing on a single monitored parameter, a much more meaningful approach is to incorporate multiple sensor data streams and investigate if the deteriorated COP is due to the degradation of system internal health state, rather than the variation of some uncontrollable external factors. Incorporating system physics into data-driven models is essential to address this challenge.
• System health state variables are correlated and subject to complex temporal dynamics. During the operation of an engineering system, the internal health states of the system gradually deviate away from their nominal values. In the thermodynamic model (1), for example, although the nominal values of the system health states (𝛾 1 , 𝛾 2 , 𝛾 3 ) are designed into the system, the actual values of (𝛾 1 , 𝛾 2 , 𝛾 3 ) inevitably degrade over time, leading to deteriorated system performance. Hence, advanced stochastic models are needed to capture the correlated temporal dynamics among multiple system states in the absence of sufficient physical knowledge. The correlation among system states is often relevant when the temporal dynamics of these states are driven by some common, but unknown, underlying operating conditions, environmental process, or external shocks [12,13].
The aforementioned challenges are commonly faced by the health prognostics of a wide range of engineering systems, where the integration of fundamental system physics with data-driven models is essential for generating transparent, interpretable and actionable engineering insights.
The modeling of degradation data has been extensively investigated when degradation signals are directly observed. Meeker et al. [14] described the random-coefficients General Path Model (GPM) for degradation data. Based on GPM, Hong et al. [15] proposed a statistical method for degradation data modeling with dynamic covariates and presented an application to outdoor weathering data. Recently, Kim and Liu [16] proposed a deep learning framework that incorporates the general characteristics of degradation processes and provides the interval estimation of remaining useful life. Following the early work of Birnbaum and Saunders [17], Bhattacharyya and Fries [18], Doksum and Hoyland [19], stochastic processes have also been utilized to approximate real-world degradation processes; see e.g., [20][21][22][23][24]. The modeling of degradation data under dynamic environments has also received much attention [15,[25][26][27][28][29][30]. Comprehensive reviews of existing models are available from [30,31]. In our case, however, system health states are not directly observed and need to be estimated from sensor signals while invoking system working principles (i.e., the first and second challenges above). Hence, the above-mentioned degradation models do not automatically apply. If the degradation of hidden system state can be firstly estimated, then, one may apply the existing approaches for follow-up actions. For example, [32] investigated the optimization of on-condition failure thresholds and inspection intervals for multi-component systems with each component experiencing multiple failure processes due to simultaneous exposure to degradation and shock loads. We also note that, there exist approaches to fuse multiple signals to construct a composite Health Index that can then be modeled by degradation models [33][34][35]. In our case, however, different signals monitor different system operating parameters with different physical meanings. Hence, it is no longer appropriate to directly fuse these sensor signals into a univariate health index, and the physical connections among these sensor signals are lost during this process.
To tackle the first two challenges above, this paper proposes a physics-regularized framework for health prognostics of complex engineered systems with multiple hidden health states. The approach consists of two layers: a data layer and a physics layer. The data layer captures the temporal dynamics (e.g., possible degradation or drift) of multiple system health states by a 𝑛th order Linear Time-Invariant (LTI) Stochastic Differential Equation. The physics layer, on the other hand, imposes regularization over system health states, by invoking the governing relationship among the distributions of observables (i.e., sensor monitoring data) and system health states. In Fig. 3. A physics-regularized data-driven approach for the health prognostics of complex engineered systems with dependent health states. particular, this layer establishes the conditional distribution of (observed) system performance, given (hidden) system health states as well as (observed) operation parameters. The idea is sketched in Fig. 3. Under this framework, the integration of system physics and sensor data is achieved in a non-intrusive manner in the sense that system physics serves as a soft constraint or regularization.
The framework above leads to a dynamic model (or, state-space) model to be described in the next Section. In the literature, [36] proposed a second-order polynomial dynamic linear model to characterize the growth of the depth of corrosion defects on energy pipelines. The model does not consider multiple sensor signals and is purely datadriven. Wang et al. [37] investigated the modeling and forecasting of temperature-induced strain of a long-span bridge using an improved Bayesian dynamical (state-space) linear model that involves autoregressive, trend, seasonal and regression components. This approach is not used for estimating hidden system health degradation by utilizing multiple sensor signals and does not require system governing physics to be integrated. Li et al. [38] proposed a two-factor state-space model for remaining useful life prediction under time-varying operating conditions. A single state variable is considered and the governing physics is not explicitly used to construct the state-space model. Veloso and Loschi [39] utilized a dynamic linear degradation model to deal with the heterogeneity in degradation paths. The model can be applied to degradation modeling where a univariate degradation signal is directly observed (which is not our case), and does not consider multiple sensor signals and the physics that links the monitored parameters. Skordilis and Moghaddass [40] proposed a novel generative framework for failure prognosis utilizing a hybrid state-space model that represents the evolution of system operating condition and its degradation over time. A single-layer feed-forward neural network is employed to model the nonparametric relationship between the multi-dimensional observation process and system dynamics.
Unlike the approaches reviewed above, the proposed physics-regularized framework leverages the governing system physics to directly construct the measurement equation that links multiple sensor signals. The degradation of multiple hidden system states are captured by Stochastic Differential Equations which give rise to the state equation. In addition, we particularly consider the statistical dependency among multiple system states, and model the dependency using a non-parametric approach based on the Archimedean family of copulas [41]. Unlike the existing work where a specific parametric copula function is often used [42,43], the Archimedean family of copulas includes the most commonly used copula functions (e.g., Clayton, Gumbel, Frank, Joe, etc.), and thus provide more flexible models considering potentially complex dependence structures among hidden system states (which may not be adequately captured by a specific parametric copula function). On the other hand, the use of non-parametric copula functions increases the computational complexity as more parameters need to be estimated. A hybrid Gibbs sampler based on the Forward Filtering Backward Sampling (FFBS) is developed to perform the statistical inference.
Finally, the connection between the proposed framework and the Gaussian Process regression is presented, connecting the proposed approach to a large body of literature in machine learning. The proposed approach is applied to solve a real problem with real datasets, demonstrating the significant potential of physics-informed machine learning for reliability and safety-the main theme of this special issue.
The remainder of the paper is organized as follows. Section 2 presents the proposed framework. A case study based on a real dataset is presented in Section 3 to illustrate the application of the proposed approach. Section 4 concludes the paper.
This section presents the physics-regularized statistical modeling framework for health prognostics of complex engineered systems with multiple hidden health states, utilizing both system physics and sensor monitoring data. In particular, we let 𝛼 𝑖 (𝑡), 𝑖 = 1, 2, … , 𝑚, denote the 𝑖th hidden system health state, and let 𝜶(𝑡) = (𝛼 1 (𝑡), 𝛼 2 (𝑡), … , 𝛼 𝑚 (𝑡)) 𝑇 be a 𝑚-dimensional continuous-time time-series that contains all state variables. For any state variable 𝑖, we further define a vector 𝒂 𝑖 (𝑡) = (𝛼 𝑖 (𝑡), 𝑑 𝑑𝑡 𝛼 𝑖 (𝑡), ⋅, 𝑑 𝑛-1 𝑑𝑡 𝑛-1 𝛼 𝑖 (𝑡)) 𝑇 that consists of the 𝑖th state variable 𝛼 𝑖 (𝑡) and its derivatives. A collection of 𝒂 𝑖 (𝑡) for all 𝑖 = 1, 2, … , 𝑚 is denoted by
The proposed framework consists of two layers: a data layer and a physics layer.
The data layer captures how the (unobserved) system health state 𝜶(𝑡) evolves over time during the operation of an engineering system. Degradation, for example, is one of the main reasons that causes the system health state variables to drift away from their nominal values, leading to deteriorated system performance. If the system is properly working, we consider a generic scenario where the temporal dynamics of the 𝑖th hidden state, i.e., 𝛼 𝑖 (𝑡) for 𝑖 = 1, 2, … , 𝑚, is governed by a 𝑛th order Linear Time-Invariant (LTI) Stochastic Differential Equation [44]:
for 𝑖 = 1, 2, … , 𝑚. Here, 𝑤 𝑖 (𝑡) ∼ 𝑁(0, 𝜎 2 𝑤 ) is the Gaussian white noise, and {𝜇 𝑗 } 𝑛-1 𝑗=0 are the coefficients. Note that, in Section 2.4, we will introduce the statistical correlation among 𝑤 𝑖 (𝑡) for 𝑖 = 1, 2, … , 𝑚, such that the degradation paths of multiple health variables are correlated. The differential Eq. ( 2) captures the dynamics of the hidden system state and has been widely applied to a spectrum of engineering applications such as image processing, vibration, circuits, signal processing and control [44,45]. In a special case when 𝑛 = 1 and 𝜇 0 = 0, Eq. ( 2) reduces to a simpler form of a stochastic model, α𝑖 (𝑡) = 𝑤 𝑖 (𝑡), where 𝛼 𝑖 (𝑡) becomes a Wiener process-a widely adopted model for univariate degradation processes.
When the state dynamics drifts away from its nominal condition during the operation of the system (such as aging, malfunction of certain components, etc.), the system performance is expected to deteriorate. To capture such a drift, we introduce a term 𝛽 𝑖 (𝑡) on the right side of (2), and obtain [Abnormal state dynamics with shift]:
It is important to note that, let
for each system state 𝑖, 𝑖 = 1, 2, … , 𝑚, the 𝑛th order LTI differential Eq. ( 2) has a state-space representation as follows:
where 𝑮 𝑖 is the feedback matrix,
𝑳 = diag(0, 0, … , 1) is the noise matrix, 𝒃 𝑖 (𝑡) is a 𝑛-dimensional vector that captures the potential deviation of state dynamics from the nominal condition, 𝒘 𝑖 (𝑡) ∼ 𝑁(𝟎, 𝜮 (𝑖) 𝒘 ) is the 𝑛-dimensional Gaussian white noise, and the system health state, 𝛼 𝑖 (𝑡), is recovered from 𝒂 𝑖 (𝑡) through a 1 × 𝑛 matrix 𝑯 = [1, 0, … , 0]. Note that, the statistical dependency between multiple system health state will be formally introduced in Section 2.4.
Because 𝑮 𝑖 linearly operates on 𝒂 𝑖 (𝑡) in ( 5), the differential equation in the first line of (5) can be solved at discrete times, and we obtain:
where exp(⋅) is the matrix exponential, 𝒃 𝑖 (𝑡)𝛥 is the first-order approximation of the total amount of shift ∫ 𝑡+𝛥 𝑡 𝒃 𝑖 (𝜏)𝑑𝜏 over a time interval with length 𝛥, 𝒒 𝑖 (𝑡) ∼ 𝑁(𝟎, 𝜮 (𝑖) 𝒒 ) and
The physics layer imposes regularization on how individual system health variables, 𝛼 𝑖 (𝑡) for 𝑖 = 1, 2, … , 𝑚, are physically connected following some fundamental system working principles. In particular, the main goal of the physics layer is to link the distribution of the (monitored) system responses to the (monitored) operating parameters, given (unobserved) system health state and (known) system physics:
where [⋅|⋅] represents the conditional density, 𝒚(𝑡) is a vector of system responses (observed), 𝒙(𝑡) = (𝑥 1 (𝑡), 𝑥 2 (𝑡), ⋅, 𝑥 𝑑 (𝑡)) 𝑇 is the 𝑑-dimensional streaming observations of critical system operating parameters (observed), and 𝜶(𝑡) = (𝛼 1 (𝑡), 𝛼 2 (𝑡), ⋅, 𝛼 𝑚 (𝑡)) 𝑇 is a 𝑚-dimensional system health state at time 𝑡 (not observed). Eq. ( 9) outlines a generic case which is universally relevant to almost all designed engineering systems.
The specification of the physics layer requires the construction of a mapping, 𝑓 , such that:
where 𝒗(𝑡) captures the measurement error. The mapping 𝑓 is constructed from known system physics and see Section 1.1 for a real example.
Let ã(𝑡) = (𝒂 1 (𝑡), … , 𝒂 𝑚 (𝑡)) 𝑇 , b(𝑡) = (𝒃 1 (𝑡)𝛥, … , 𝒃 𝑚 (𝑡)𝛥) 𝑇 and q(𝑡) = (𝒒 1 (𝑡), … , 𝒒 𝑚 (𝑡)) 𝑇 , we obtain a dynamic model by integrating the physics layer and data layer:
where G = diag(exp(𝑮 1 𝛥), exp(𝑮 2 𝛥), … , exp(𝑮 𝑚 𝛥)), 𝜮 q(𝑡) (𝑡) = diag(𝜮 (1) 𝒒 , 𝜮 (2) 𝒒 , … , 𝜮 (𝑚) 𝒒 ), H = diag{diag(𝑯)} and 𝜮 𝝐 = diag(𝜎 2 𝜖 ). If the mapping, 𝑓 , is linear or can be approximately by a linear operation such that 𝒚(𝑡) = 𝒙 𝑇 (𝑡)𝜶(𝑡) + 𝒗(𝑡) ≡ 𝑭 (𝑡)𝜶(𝑡) + 𝒗(𝑡), we re-write (11) as
where F (𝑡) = 𝑭 (𝑡) H. In fact, by defining a mapping 𝑔(𝑡) = F ã(𝑡), the dynamic model ( 12) is the state-space representation of a Gaussian Process () regression problem with the following form [45]:
where the function, 𝑔(𝑡), is a realization of a random prior with a specified covariance function 𝑘(⋅, ⋅). From the function-space perspective, a is a collection of random variables and any finite number of which have a joint Gaussian distribution [46]. The covariance function, 𝑘(⋅, ⋅), at the stationary state, can be computed by:
where 𝑬(𝑡) = exp( G(𝑡)) and 𝑷 ∞ solves the Riccati equation
Hence, the dynamic model proposed in this paper can be interpreted as a regression problem (13) which provides a powerful modeling approach in both statistics and machine learning. On the other hand, the dynamic model (12) provides major computational advantages rooted in its conditional structure.
Finally, we establish the statistical dependency among hidden system state variables, i.e., 𝛼 1 (𝑡), 𝛼 2 (𝑡), … , 𝛼 𝑚 (𝑡). In the model above, the dynamics of each system state is governed by the differential Eqs. (3). Hence, by introducing statistical dependency between 𝑤 1 (𝑡), 𝑤 2 (𝑡), … , 𝑤 𝑚 (𝑡), the statistical correlation among the system states can be naturally established.
Let 𝑋 1 , 𝑋 2 , . . . , 𝑋 𝑚 be random variables with joint distribution function 𝐹 (joint) and marginals 𝐹 1 , 𝐹 2 , . . . , 𝐹 𝑚 , respectively. Then, there exists a copula 𝐶 such that
If the marginals are continuous, then, the copula 𝐶 is unique; Otherwise, it is uniquely determined on Ran(𝐹 1 ) × Ran(𝐹 2 ) × ⋯ × Ran(𝐹 𝑚 ), where Ran(𝐹 ) denotes the range of 𝐹 [41].
Based on the Sklar's Theorem, we let 𝑤 1 (𝑡), 𝑤 2 (𝑡), … , 𝑤 𝑚 (𝑡) be 𝑚 random variables with joint distribution function 𝐹 (joint) and continuous marginals 𝐹 𝑤 1 (𝑡) , 𝐹 𝑤 2 (𝑡) , … , 𝐹 𝑤 𝑚 (𝑡) respectively. Then, there exists a unique copula 𝐶 such that
As a tool for statistical analysis, copulas allow for the modeling of marginals to be handled separately from the dependence structure characterized by the copula, and represent a flexible alternative in which one can bypass the complex specification and validation of multivariate distributions. Although there exist many candidate copulas in the literature, the choice of a particular parametric copula function for a particular problem is still challenging. For potentially complex dependence structures, a specific type of parametric copula may not be adequate.
Hence, to make our model robust and general, we consider a large family of copulas known as Archimedean [41]. Archimedean copulas are an associative class of copulas that model the dependence in arbitrarily high dimensions with only one parameter, which governs the strength of statistical dependence. The most prominent bivariate Archimedean copulas include Clayton, Gumbel, Frank, Joe, etc.
From the modeling point of view, one main advantage of Archimedean copulas is that any Archimedean copula 𝐶 admits the following representation:
where 𝜑(𝑢) is known as the generator function, which is strictly decreasing and convex on (0, 1) such that 𝜑(1) = 0. Suppose that it is possible to approximate the generator 𝜑(𝑢) by some function φ(𝑢). Then, we may restrict our attention to the inference on the approximate function without choosing a specific parametric form of the copula function. Instead of directly approximating the generator function which is unbounded at 0 + , we adopt the idea proposed in [47] which approximates the following function
where
using a cubic B-splines given by
where 𝐁 𝑢 contains the values of the B-splines at 𝑢 based on 𝑘 equidistant internal knots on [0, 1], and 𝜼 = (𝜂 0 , 𝜂 1 , … , 𝜂 𝑘+1 ) 𝑇 ∈ R 𝑘+2 . Given a knot sequence, 𝑢 0 = 𝑢 1 < ⋯ < 𝑢 𝑘 = 𝑢 𝑘+1 , the 𝑖th B-splines of order 𝑛 𝐵 at time 𝑢 is obtained using the standard recurrence:
)𝐵 𝑢 (𝑖 + 1, 𝑛 -1) for 𝑖 = 0, … , 𝑘 + 1 and 𝑛 = 1, … , 𝑛 𝐵 [48]. In particular, [47] showed that the elements of 𝜼 must satisfy the following conditions for the approximation (21) to be valid:
Introducing statistical dependency among 𝛼 1 (𝑡), 𝛼 2 (𝑡), … , 𝛼 𝑚 (𝑡) does not alter the parametric structure of the dynamic model (12). However, it does change the covariance matrix of q(𝑡).
Let
where L = diag(𝑳, 𝑳, … , 𝑳) and W (𝑡) = (𝑾 𝑇 1 (𝑡), 𝑾 𝑇 2 (𝑡), … , 𝑾 𝑇 𝑚 (𝑡)) 𝑇 . Solving the linear stochastic differential Eq. ( 23) at discrete times, we obtain
which maintains the same form of ( 12), but the covariance matrix of q(𝑡) is given by:
with
Given the special structure of 𝑳 defined under (5), L𝜮 W L𝑇 is a 𝑚𝑛 × 𝑚𝑛 sparse matrix which has non-zero entries only at its (𝑖 1 𝑛, 𝑖 2 𝑛)th positions for 𝑖 1 , 𝑖 2 = 1, 2, … , 𝑚. Hence, given the observation 𝒚(𝑡), it is possible to obtain the posterior distribution (i.e., filtering distribution) of the hidden system health variables [ ã(𝑡), b(𝑡)|𝒚(𝑡)] from the dynamic model (24). The obtained posterior distribution enables one to monitor the temporal dynamics of the critical system health conditions.
In this case study, we re-visit the motivating example presented in Section 1.1. The data used in this case study are shown in Fig. 2, including the observed daily COP, condenser coolant inlet temperature 𝑇 𝑐 , evaporator coolant outlet temperature 𝑇 𝑒 , and cooling capacity 𝑄 𝑒 over a 57-day study period (see Section 1.1 for more detailed descriptions). There exists a strong thermodynamics law that governs the relationship among these critical operating parameters,
; also see (1) for more details. Here, the true system health state, i.e., the internal irreversibility states (𝛾 1 , 𝛾 2 , 𝛾 3 ), are not directly observed and may gradually drift away from their nominal values. Hence, the goal is to estimate the (statistically correlated) degradation of the hidden system state variables from the monitored operating parameters.
The physics layer is constructed from the governing physics (1). Firstly, we note from Fig. 2 that the daily evaporator coolant outlet temperature, 𝑇 𝑒 , presents very small variability over the 57-day observation period (less than 0.5𝐾). In fact, 𝑇 𝑒 is pre-set and should remain at a certain level. From the physics point of view, this observation suggests that the chiller can still provide chilled water with the preset temperature, although the energy efficiency in cooling the water might have already decreased as indicated by the drop of COP. Hence, by treating 𝑇 𝑒 as a constant, we obtain from the governing physics (1) a linear model:
)
where 𝛼 1 = -𝛾 1 and 𝛼 2 = 𝛾 2 -𝛾 3 𝑇 𝑒 (𝑡) -1 are the two hidden system states, 𝑥 𝑡 = 𝑇 𝑐 (𝑡), and 𝑣(𝑡) ∼ 𝑁(0, 𝜎 2 𝑣 ) is the observation error; see (10). Following (3), the (marginal) temporal dynamics of the two system states are modeled by a first-order LTI stochastic differential equation such that α1 (𝑡) = 𝛽 1 (𝑡) + 𝑤 1 (𝑡) and α2 (𝑡) = 𝛽 2 (𝑡) + 𝑤 2 (𝑡) where 𝛽 1 and 𝛽 2 capture the drift of the two state variables from their nominal values. In addition, we establish the statistical dependency between 𝑤 1 (𝑡) and 𝑤 2 (𝑡) by ( 17) through a copula 𝐶 such that
where 𝐹 (joint) and 𝐹 respectively denote the joint and marginal distributions. The function 𝐶 corresponds to an Archimedean copula and its generator function is modeled using cubic B-splines; see ( 18)- (21). Note that, when the temporal dynamics of the two system states are modeled by a first-order LTI stochastic differential equation, ã(𝑡) and b(𝑡) in (24) are respectively (𝛼 1 (𝑡), 𝛼 2 (𝑡)) 𝑇 and (𝛽 1 (𝑡), 𝛽 2 (𝑡)) 𝑇 . Let 𝜽(𝑡) = (𝛼 1 (𝑡), 𝛼 2 (𝑡), 𝛽 1 (𝑡), 𝛽 2 (𝑡)) 𝑇 , and let 𝛽 1 (𝑡) and 𝛽 2 (𝑡) be two AR(1) processes, we obtain a dynamic model as follows:
where
,
In (28), the observed response 𝑦(𝑡) is determined by the latent process 𝜽(𝑡) up to a Gaussian error. Note that, the drifts 𝛽 1 and 𝛽 2 are also treated as auxiliary state variables, although they are not the actual system health state variables [49]. The augmented state variables evolve over time following a Markovian structure, and the statistical dependency among the temporal dynamics of 𝛼 1 and 𝛼 2 is captured by 𝜎 𝛽 1 ,𝛽 2 . Based on the Hoeffding's Lemma [41]
Let 𝝍 = (𝜎 𝑣 , 𝜎 𝛼 1 , 𝜎 𝛼 2 , 𝜎 𝛽 1 , 𝜎 𝛽 2 , 𝜎 𝛽 1 ,𝛽 2 ) be a collection of the unknown parameters. Note that, the last parameter 𝜎 𝛽 1 ,𝛽 2 depends on a set of unknown B-splines coefficient 𝜼 that defines the copula function. Given the observations 𝑦 1∶𝑇 = (𝑦(1), 𝑦(2), … , 𝑦(𝑇 )), the posterior distribution of the parameter and unobservable states is 𝜋(𝜽 0∶𝑇 , 𝝍; 𝑦 1∶𝑇 ) = 𝜋(𝜽 0∶𝑇 |𝑦 1∶𝑇 , 𝝍)𝜋(𝝍|𝑦 1∶𝑇 ), and the Gibbs sampling from 𝜋(𝜽 0∶𝑇 , 𝝍; 𝑦 1∶𝑇 ) requires one to simulate from the full conditional densities 𝜋(𝜽 0∶𝑇 ; 𝑦 1∶𝑇 , 𝝍) and 𝜋(𝝍; 𝑦 1∶𝑇 ). Details are provided as follows.
It is noted that, although the parameters 𝝍 1 = (𝜎 𝑣 , 𝜎 𝛼 1 , 𝜎 𝛼 2 , 𝜎 𝛽 1 , 𝜎 𝛽 2 ) can be efficiently sampled leveraging the well-known conjugate inverse Gamma priors, the sampling of 𝜎 𝛽 1 ,𝛽 2 requires the drawing 𝜼 that defines the copula function; see (30). Because 𝜂 0 = 𝜂 𝑘+1 = 0, one may sample 𝑘 -2 B-splines coefficients, η = (𝜂 1 , … , 𝜂 𝑘 ), using a Bayesian framework described in [50].
Let 𝑢 1 (𝑡) = 𝛽 1 (𝑡 + 𝛥) -𝛽 1 (𝑡) and 𝑢 2 (𝑡) = 𝛽 2 (𝑡 + 𝛥) -𝛽 2 (𝑡), the likelihood function of η is
where 𝐶 = 𝐶(𝑢 1 (𝑡), 𝑢 2 (𝑡); η) = 𝜑 -1 (𝜑(𝑢 1 (𝑡); η) + 𝜑(𝑢 2 (𝑡); η); 𝜼) and 𝜑 is the generator function defined in (18). Adopting a non-informative prior for η
the posterior distribution of η is
Once the unknown parameters have been sampled, the state variables can be efficiently drawn from 𝜋(𝜽 0∶𝑇 |𝑦 1∶𝑇 , 𝝍) using the wellknown Forward Filtering Backward Sampling (FFBS) with linear complexity in time and the number of state variables [51]. Algorithm 1 summarizes the steps that sample from the full conditional densities 𝜋(𝜽 0∶𝑇 ; 𝑦 1∶𝑇 , 𝝍) and 𝜋(𝝍; 𝑦 1∶𝑇 ).
Applying the algorithm above to the dataset, Fig. 4 shows the posterior distributions of the model parameters 𝝍 = (𝜎 𝒗 , 𝜎 𝛼 1 , 𝜎 𝛼 2 , 𝜎 𝛽 1 , 𝜎 𝛽 2 , 𝜎 𝛽 1 ,𝛽 2 ). Fig. 5 shows the posterior means as well as the 95% bootstrap confidence intervals of the state variable, 𝜽(𝑡). It is immediately seen that both 𝛼 1 (𝑡) and 𝛼 2 (𝑡) gradually shift away from their initial values over the 57-day monitoring period, indicating deteriorating system internal health. In particular, the amount of daily shift 𝛼 1 (𝑡) is captured by 𝛽 1 (𝑡), while the amount of daily shift 𝛼 2 (𝑡) is captured by 𝛽 2 (𝑡).
In our model, the statistical dependency of system state variables is established through a copula function with its generator function being modeled by non-parametric cubic B-splines. The purpose is to bypass the difficulty of specifying a parametric copula function, and enhance the modeling flexibility of our model. Fig. 6 shows the posterior mean of the function, 𝜆 in (19), estimated from the proposed model (black thick line). The posterior mean is obtained by averaging the samples of 𝜆 (the gray lines show 50 selected samples for visualization purposes). For comparison purposes, we also re-fit our model using parametric copula functions, Clayton, AMH (Ali-Mikhail-Haq) and Frank, as well as assuming independent system health state. The idea is that, if the shape of 𝜆 obtained using the non-parametric approach is similar to that obtained from a parametric approach, then, a parametric copula Algorithm 1: FFBS in a hybrid sampler Initialize 𝝍 0 and 𝜼 0 such that 𝜼 0 satisfies the convexity and sign conditions in (22). for 𝑙 = 1 ∶ 𝑁 do Draw 𝜽 (𝑙) 0∶𝑇 from 𝜋(𝜽 (𝑙) 0∶𝑇 |𝑦 1∶𝑇 , 𝝍 (𝑙-1) ) using FFBS; Draw 𝝍 (𝑙) 1 from the posterior distribution assuming conjugate inverse-Gamma priors; Obtain (𝑢 1 (𝑡), 𝑢 2 (𝑡)) 𝑇 𝑡=1 , and set . else set 𝝃 𝑗 = 𝝃 𝑗-1 end end end Set 𝜼 (𝑙) = 𝝃 𝑘 and compute 𝜎 𝛽 1 ,𝛽 2 from (30); end should be used instead of the non-parametric one that increases the computational complexity. However, this is not the case as shown in Fig. 6. The comparison in Fig. 6 shows that the shape of 𝜆 obtained from the non-parametric approach is more complex and cannot be well captured by any of these commonly used parametric copula function, showing an improved modeling capability using the Archimedean family of copulas.
Finally, the normality assumption of the model is validated. The first column of Fig. 7 shows the Q-Q plot of the residuals from the observation equation in (28). The five Q-Q plots in this columns are respectively based on five samples of 𝜽 drawn from the FFBS. Columns 2 to 4 of Fig. 7 shows the Q-Q plots of the residuals from the state transition equation in (28) for the four state variables within 𝜽. Similarly, the five Q-Q plots in this columns are respectively based on five samples of 𝜽 drawn from the algorithm. The normality assumption of model ( 28) is well justified.
As discussed in Section 1, the decrease of COP does not necessarily imply that the system internal health states (hidden) have degraded. COP depends on various external environmental factors and internal system health state. The change of working load and outdoor temperature may cause the drop of the observed COP, while the system is functioning properly. Hence, before DC engineers can stop the operation of the cooling system, it is necessary to understand if the deteriorated COP is indeed due to the degradation of system internal health state, rather than the variation of uncontrollable external factors (note that, it is often costly and risky to stop the normal operation of DC cooling systems without convincing evidence). The proposed model successfully addresses this question by revealing the degradation of hidden system health states; see Fig. 5. It is worth noting that, the estimated system health states possess well-defined physical interpretation and are related to system irreversibility states in (1). This finding provides interpretable justifications that DC engineers could stop the operation of the DC cooling system, and investigate the root causes behind system state degradation. In our case study, DC engineers eventually located the root cause behind the observed cooling performance deterioration: the debris in the water pipe (see Fig. 8). Note that, in air conditioning systems, chillers are utilized to provide cooling water which is distributed to cool the air in server rooms of a DC. Used water is re-circulated back to the chiller to be cooled again. When the debris blocks the water flow, the heat exchange rate drops, causing the drift of system internal irreversibility states. The debris was brought to the water pipe due to a design deficiency during the DC construction phase and was later fixed. The discovery of such actionable insights can be facilitated when system physics is incorporated into the proposed health prognostics of complex engineered systems with multiple hidden health states.
This paper proposed a physics-regularized data-driven approach for health prognostics of complex engineered systems with multiple hidden health states. The proposed methodologies enabled the integration of critical system working principles with streaming sensor observations. The proposed framework consists of a data layer and a physics layer. The data layer captures the statistically correlated temporal dynamics of hidden system states, while the physics layer imposes regularization on the system health states by invoking the physical relationship between multiple observed system operating parameters. The integration of physics and data-driven approaches is thus achieved in a non-intrusive manner. The non-parametric cubic B-splines has been successfully employed to describe the complex statistical correlations among system state variables. The application and the effectiveness of the proposed approach have been demonstrated by a case study based on a real data set.