Sensor Data Reconstruction in Industrial Environments with Cellular Connectivity

The reliable acquisition of monitoring information is critical for several industrial use cases relying on wireless sensor network deployments. However, missing sensor measurements are typical in industrial systems empowered by cellular connectivity due to the stochastic nature of the wireless channel. In this paper, we propose a sensor data reconstruction scheme that exploits the hidden data dynamics to accurately estimate the missing measurements. Based on an analytical framework for the network model and a closed-form expression for the outage probability, the impact on the reconstruction error performance is thoroughly explored. Considering a dataset with high spatiotemporal correlation in the sensor observations, our proposed scheme is shown to outperform two baseline data recovery methods in terms of reconstruction error for various network configurations. In addition, despite the presence of imperfect cellular connectivity, our proposed scheme exhibits high event-detection accuracy.


I. INTRODUCTION
The ongoing modernization of the aging industrial systems aims at the transformation of the entire production chain, composed of processes, machines, workers, and production lines, into a fully integrated, automated and interconnected paradigm to improve productivity and operational efficiency. This evolution resides at the core of the upcoming Industry 4.0 era and heavily relies on the digitalization and ubiquitous connectivity provided by 5G mobile systems and Internet of Things (IoT) technologies [1]- [3]. The integration of 5G ultra-reliable low-latency communication (URLLC) as well as the massive deployment of intelligent wireless sensors are expected to: i) unlock unprecedented industrial use cases, e.g., non-stationary machines, scalable line configuration, etc.; ii) allow the expansion of digital operations, e.g., for remote plant condition monitoring; and iii) contribute to the extension of machine lifetime, e.g., via predictive maintenance and timely anomaly detection.
Instrumental to the ongoing industrial transformation is the availability and accessibility of well-structured datasets related with the monitoring and condition maintenance of industrial equipment [4]. However, in industrial IoT systems empowered by wireless connectivity, the unreliable nature of the wireless medium constitutes one of the fundamental challenges introducing uncertainties in various state-monitoring operations. The complex fading conditions in factory plants, usually rich in metallic surfaces and physical obstructions, result in high network dynamics and harsh radio propagation environments which significantly affect the accuracy of the transmitted sensor measurements and often result in transmission failures. In addition, industrial datasets often exhibit high spatiotemporal correlation with context dependencies that render the data characterization and modeling tasks non-trivial.
During the recent years, time series analysis has attracted considerable research attention for the accurate detection, identification and diagnosis of abnormal behaviours, e.g., faulty conditions, disturbances, etc., in industrial environments [5]- [12]. In the majority of existing works, however, the stochastic nature of the wireless channel and its impact on the sensor measurements is either neglected [5]- [9] or limitedly considered [10]- [12]. In fact, it is a common approach in the literature to consider simplified models for the wireless channel reliability; either with the use of a generic expression for the probability of failure, e.g., as in [10], [11], or with the aid of indicator functions for the unsuccessful transmissions, e.g., as in [12], and without explicit network topology considerations. Instead, a more rigorous analysis of the communication impairments is necessary for a thorough assessment of the impact of various network parameters on the quality of the sensor measurements. The authors in [13] employ a Kalman filter based state estimation model for industrial automation applications and propose an optimized small-packet transmission scheme to improve the networkwide revenue. Even though erroneous sensor transmissions are incorporated in the system model, missing data reconstruction is not explicitly addressed and the focus is solely on the optimization of the small-packet transmission scheme. In addition, to the best of the authors' knowledge, there is no prior work studying the recovery of missing sensor measurements transmitted over a wireless channel under the presence of spatiotemporal dynamics in the sensors' observations. Contribution: Motivated by these literature gaps, our contribution in this paper is twofold: representation of an incomplete sensor data sequence with missing measurements due to the unreliable nature of the wireless channel. • We develop a data reconstruction scheme that takes into account the hidden dynamics, i.e., spatiotemporal context dependencies, in the dataset to estimate the missing sensor measurements. The performance assessment in terms of reconstruction error reveals significant reconstruction accuracy gains against baseline methods and offers useful insights for the design of industrial IoT systems empowered by cellular connectivity. Organization: The rest of the paper is organized as follows. The system model, assumptions and the methodology for the connectivity analysis are presented in Section II. Section III describes in detail the proposed missing data reconstruction scheme along with all the involved steps of an iterative expectation-maximization (EM) algorithm employed for the estimation of the missing sensor measurements. Model validation and performance comparison of our proposed scheme against baseline methods for missing data recovery are presented in Section IV. Finally, Section V concludes the paper.

A. Network Model
We consider an industrial process monitoring scenario composed of a data fusion center and a number of deployed sensors able to monitor the operational state and condition of the various components. The sensors are equipped with radio interfaces and transmit their state measurements to the fusion center which is co-located with their associated base station 1 (BS). A stochastic geometry framework is considered for the positions of the BSs and the sensors. In particular, we assume that the BSs are distributed on the plane according to a homogeneous Poisson point process (PPP) Ψ b of intensity λ b , while the sensors are spatially distributed according to a homogeneous PPP Ψ u = {u k ; i = 1, 2, 3, . . . } with intensity λ u . We further assume saturated uplink traffic conditions where all sensors always have data to transmit while λ u is considered high enough such that each BS will serve at least one sensor per channel.

B. Channel Model and Power Control
Regarding the channel model, a generic power-law pathloss model is considered where the signal power decays with the propagation distance r at the rate r −α , where α denotes the path-loss exponent. In addition to the path-loss attenuation, the channel power gains are assumed to be independent of each other and of the spatial locations, and identically distributed (i.i.d). We assume Rayleigh fast fading where the channel power gains, denoted by h, are exponentially distributed with unit mean. Let p k denote the transmit power of sensor k. All sensors are considered to have equal maximum transmit power denoted by p max and employ a truncated channel inversion power control policy to compensate for the path loss and 1 A closest BS association scheme is adopted in this paper. keep the average received signal power at the BS equal to a threshold ρ which is greater than the receiver sensitivity ρ min .

C. Outage Probability
If the transmit power required for the path-loss inversion is higher than p max , the sensor is considered to be in outage due to insufficient transmit power and does not transmit its measurement. Based on the transmit power analysis conducted in [14], the probability that a sensor experiences outage due to insufficient transmit power for α = 4 is given by From Eq. (1), it can be observed that P o1 depends on the BS intensity λ b , the maximum sensor transmit power p max , and the threshold ρ which constitutes a network design parameter. A measurement from an active sensor (i.e., a sensor not in power outage) is considered to be successfully decoded at the serving BS if and only if the signal-to-interference-plusnoise ratio (SINR) of the useful signal is greater than a certain threshold γ th . Based on the stationarity of the PPP, the SINR experienced at the origin BS for a target active sensor can be expressed as where the numerator corresponds to the useful signal power, σ 2 is the noise power while the second term in the denominator denotes the aggregate interference from other uplink sensor transmissions on the same channel 2 . Then, according to the SINR analysis conducted in [14] and for α = 4, the SINR outage probability P o2 = P(SINR < γ th ) for an active sensor can be calculated as where γ(a, b) = b 0 t a−1 e −t dt denotes the lower incomplete gamma function. The total outage probability can then be expressed as where P o1 and P o2 can be calculated using Eqs. (1) and (3) respectively.

III. SENSOR DATA RECONSTRUCTION
At the fusion center located at each serving BS, the received measurements can be represented by a partially observable time sequence Y = [y 1 , y 2 , . . . , y T ], where each vector y t contains the received measurements at time-step t from the deployed sensors. In particular, by incorporating the outage probability P o (given by Eq. (4)), the entries y t,k of y t can be expressed as where x t,k corresponds to the transmitted measurement from sensor k at time-step t while υ t,k is a zero-mean white Gaussian noise with covariance matrix Ω. It is noted that Eq. (5) takes into account the stochastic nature of the wireless channel which may result in a received vector y t with intermittent measurements. We consider a time sequence of latent variables (i.e., hidden states) Z = [z 1 , z 2 , . . . , z T ] to model the dynamics and the hidden patterns of the received measurements. We also introduce an indicator matrix, Φ, for the missing measurements; i.e., φ t,k = 0 whenever the k-th sensor measurement in y t is missing at time t; otherwise, φ t,k = 1. Let us also denote the observed part of Y as Y r and the missing part as Y m . Following the rationale of linear dynamical systems [15], our model for the received measurements at the fusion center can be described by the following two equations: To capture temporal correlation, we assume that the latent variables at each time tick depend linearly on the previous values via the linear state transition matrix A. At each time tick, the received vector y t , including both observed and missing sensor measurements, is assumed to be a linear function of the latent variables z t via the linear projection matrix C. This mapping implicitly captures the spatial correlation among the different sensor measurements [16]. Both hidden state evolution and received measurement processes are corrupted by zero-mean white Gaussian noise, w t and v t , with covariance matrices, Q and R, respectively. Further, w t and v t are assumed to be independent. The initial state z 0 of the latent variables is also a Gaussian random variable with mean π 1 and covariance V 1 . Therefore, the parameter vector of our model is θ = (A, C, Q, R, π 1 , V 1 ).
Based on Eqs. (6) and (7), we can express the conditional probabilities for the hidden state and the received sequence, respectively, as follows: where D (ω t , µ t , Ξ) = (ω t − µ t ) Ξ −1 (ω t − µ t ) denotes the square of the Mahalanobis distance of a vector ω t with mean vector µ t and covariance matrix Ξ. Based on the Markov property implicit in the model, the factored representation of the joint probability distribution of Z and Y is given by and the joint log-likelihood can be written as Given that the received sequence Y is characterized by intermittent measurements due to the imperfect cellular connectivity, our goal is to maximize the conditional expectation of the received data log-likelihood, i.e., To that aim, we apply an iterative EM algorithm following a coordinate descent procedure [17]. We provide the details in the following subsection.

A. The EM algorithm
The EM algorithm provides an iterative method for finding the maximum likelihood estimates of θ based on the observed measurements, Y r , by successively maximizing Eq. (12). The E-step of EM algorithm requires computing L (θ) in Eq. (12). Based on Eq. (11), this computation amounts to deriving the following three expectations: Let z τ t and V τ t denote E (z t |Y τ 1 ) and Var (z t |Y τ 1 ), respectively, for the subsequence of received measurements until time τ . Note that z 1 0 = π 1 and V 1 0 = V 1 . Let also θ be an initialization of the parameter vector. The conditional expectations in Eqs. (13)-(15) can be expressed aŝ and their computation can be decomposed into the following sets of forward and backward recursion: i) Forward recursion: ii) Backward recursion: where Eq. (27) is initialized as V T T,T −1 = (I − K T C) AV T −1 T −1 . After calculating the conditional expectations of the latent variables (E-step), the M-step re-estimates the parameter vector θ to be used in the E-step. To estimate θ = (A, C, Q, R, π 1 , V 1 ), we take the respective partial derivative of Eq. (12), set to zero, and solve for the value of each respective parameter. In particular, the updated parameters are computed as follows: i) Projection matrix: ii) Measurement noise covariance: iii) State transition matrix: iv) State noise covariance: v) Initial state mean: vi) Initial state covariance: Finally, using the Markov property, the missing sensor measurements Y m can be computed from the estimation of the latent variables as The Eqs.

IV. RESULTS AND DISCUSSION
In this section, we aim to validate our proposed analytical framework against simulation results and provide a performance comparison in terms of reconstruction error of our proposed scheme with respect to two baseline methods. For the simulation setup, we consider a PPP for the BSs with intensity λ b and sensors are randomly spread over the area until all BSs are serving at least one sensor, i.e., to achieve traffic saturation conditions. Table I summarizes the parameter values used in our simulations, unless otherwise stated.
The considered dataset consists of IEC-61850 network traffic traces generated by intelligent electronic devices based on open-source libraries for protocol implementation 3 and a discrete-event network simulator [18]. The traffic streams correspond to both intra-and inter-substation communication in power systems and represent both normal and disturbance scenarios. The dataset exhibits high levels of spatial and temporal correlation across the device measurements to model the real behaviour of substation automation systems, e.g., local cascade failures result in inter-dependent transmissions of measurements from neighboring devices.
For performance comparison, we consider the following two baseline data reconstruction schemes: i) a linear interpolation method which exploits temporal continuity to compute the missing value estimates using observed sensor measurements in neighboring time ticks; and ii) a method based on Singular Value Decomposition (SVD) which transforms the weighted low-rank approximation problem to a maximum-likelihood problem with missing values and approximates them by computing the SVD iteratively [19]. In our proposed scheme, we initialize our estimated received time sequenceŶ with Y r while the missing sensor measurements are initially reconstructed by means of linear interpolation and then iteratively imputed as in Eq. (34). The process continues by updating the expectations of the latent variables based on the newly imputed values of the missing measurements until convergence. We further assume that the noise covariances in θ constitute diagonal matrices. The effectiveness of reconstruction is evaluated in terms of the mean squared error (MSE), defined as the average of the squared differences between the real and reconstructed missing measurements, i.e., (35) To reduce random effects, we repeat each simulation 10000 times and we report the average of the MSE. Fig. 1 illustrates the MSE performance with respect to the SINR threshold γ th for different reconstruction schemes. It can be observed that the MSE increases with increasing SINR threshold γ th , as the aggregate interference leads to significant signal degradation and results in a high number of sensor links in SINR outage. In addition, the increased outage probability in the high γ th regime, impacts the achieved MSE and may lead to lower performance of the data reconstruction process. However, in stark contrast with the linear interpolation and SVD-based methods, our proposed scheme keeps the MSE in relatively lower levels which is especially important for high values of γ th . We can also observe a tight match between the analytical curves and the simulation results which validates the accuracy of the outage probability expression in Eq. (4). Fig. 2 depicts the evolution of the MSE with respect to various network configurations for the three missing data reconstruction schemes. It can be observed that MSE increases with increasing ρ, as a higher received power threshold at the BS leads to an increased number of lost sensor measurements. This can be intuitively explained as follows. When ρ increases, the sensors are required to transmit at higher power levels; this, in turn, leads to an increasing number of sensors failing to establish an uplink connection with their serving BS due to  insufficient transmit power to invert their path loss, given that p max is considered constant (refer also to Eq. (1)). It is worth noting that the achieved MSE levels of our proposed scheme remain significantly lower with respect to the linear interpolation and SVD-based methods. In addition, data reconstruction performance improves with increasing BS intensity λ b . When the BSs are dense enough, the distance between a generic sensor and its corresponding serving BS decreases and the required transmit power for path-loss inversion is lower. It is worth noting that the data reconstruction gains of our scheme can be capitalized even for low values of λ b , demonstrating the robustness of our approach even for a sparse deployment of BSs. Finally, Fig. 3 illustrates the event-detection performance of our proposed scheme. In particular, as shown in the top plot of Fig. 3, we have artificially "injected" various events in the dataset, corresponding to abnormal conditions, i.e., voltage drop beyond a maximum allowable level of 3%. The superposition of multiple events in the dataset generates different patterns which translate to different model parameters θ and latent variables Z. Based on its property to compute the posterior distribution of the latent variables in each time tick, our proposed scheme is able to accurately capture the dataset dynamics, i.e., events, by tracking the peak values of the MSE shown in the bottom plot of Fig. 3. It is worth noting that the MSE peaks nearly coincide with the time ticks of the events in the top plot of Fig. 3, despite the presence of imperfect connectivity.

V. CONCLUSIONS
Industrial systems empowered by cellular connectivity heavily rely on the reliable acquisition of sensor measurements. However, due to the wireless channel impairments, the sensor measurements received at the fusion center are often highly intermittent. In this paper, a reconstruction method for the estimation of missing sensor measurements has been proposed, considering a stochastic geometry framework for the network topology in industrial environments. Based on a closed-form expression for the outage probability, the impact on the missing values' estimation accuracy is explored under different network configurations with the aid of a dataset which consists of spatiotemporal sensor measurements. The simulation results demonstrate the superior performance of our proposed scheme against two baseline data recovery approaches, while the achieved event-detection performance is shown to be highly accurate despite the incomplete received data.
Future work will aim at extending the proposed framework to capture non-linear and multi-scale spatiotemporal behavior of the sensor measurements by resorting to the utilization of switching linear Gaussian dynamical systems. In addition, we plan to further investigate the applicability of our proposed scheme for other data mining tasks, e.g., for data compression when the storing capabilities of the fusion center impose limitations on the amount of stored sensor data.