Distributed Estimation Fusion and Control with Packet Losses for Industrial Wireless Sensor and Actuator Networks

Abstract

This paper is concerned with a reliable solution to the distributed estimation and control issues over hybrid two-tier industrial wireless sensor and actuator networks (IWSANs). It aims at applying wireless technology to industrial automation control domains in hash industrial environment. A main challenge in these application domains is that wireless communication channels must satisfy strict requirements of real time and reliability. Therefore, first, we propose a new interacting dual model (IDM) adaptive estimation algorithm to identify the unreliable wireless communication channels, and, second, to further enhance the accuracy and fault tolerance of state estimation for controlled plant, a federated multisensor estimation fusion approach is presented. Furthermore, considering that the IWSANs need commonly to be deployed over a vast geographical area, a distributed collaborative control strategy is adopted. Finally, an example of temperature control illustrates the effectiveness of the proposed method.

1. Introduction

Wireless communication represents a major industrial stake in the coming years following the Fieldbus and Industrial Ethernet [1]. With the increasing demand in spatially distributed and flexible industrial applications, automation enterprises have drawn great attention to the communication and control architectures based on industrial wireless sensor and actuator networks (IWSANs). The essential attribute of IWSANs brings several significant advantages compared with traditional wired counterparts from the perspective of flexibility and convenience of the deployment for sensors and actuators (in remote, dangerous, and hard-to-reach areas), low cost (due to the less demands for cables, materials, preinstalled infrastructure, maintenance, and labor), compositionality (where scalability of network can be carried out more easily), and so forth [2, 3]. The benefits of integrating wireless communication technology with industrial applications are compelling. Some exploratory work has been done such as drip irrigation control for agriculture using wireless sensor and actuator network (WSAN) [4], environment control system in WSAN [5, 6], wireless ventilation systems in mining industrial [7], wireless fire detection and control [8], towards usage of wireless MEMS in industrial context [9], and monitoring system for robots and inventories in assembly line based on wireless sensor networks (WSNs) [10].

However, a common outstanding problem in these aforementioned applications is that the unreliable wireless communication channels suffering from delays and losses of data are difficult to ensure real-time and reliability requirements of control systems. The problem hinders that wireless technology is accepted generally and permeated in many industrial application domains. Therefore, designing the communication protocol standards, network topology architecture, and appropriate estimation and control strategies over IWSANs for networked control systems (NCSs) in industrial applications are still open issues.

In current literatures, most commonly adopted models of unreliable wireless channels are I.I.D (independent and identically distributed) model and Gilbert-Elliott model. The I.I.D model describes packet losses by a I.I.D Bernoulli random variable (taking value 1 or 0), such as in [5]. Gilbert-Elliott model considers packet losses as a two-state Markov chain to describe packet arrival and drop behaviors in [11]. Based on the above-mentioned two models, many researches have got some significant results for the problems of estimation and control [12–20].

The fault-tolerant mechanism is also an important issue which is used to improve the reliability of wireless data transmission in hash industrial context. As the current research hot spot issue of signal processing area, multisensor information fusion [21] can improve the fault tolerance and accuracy of real-time field data by utilizing multiple spatially distributed sensors to provide complementary and overlapping coverage on plants. It has received significant attention. Xia et al. [22] consider the networked data-fusion system with packet losses. In [23–26], fault-tolerant fusion, identification, and control for wireless sensor networks are presented. Hierarchical sensor data aggregation/fusion is studied by [27, 28].

Currently, high-reliable estimation and control strategies for industrial NCSs are still challenging issues. Since the quality of state estimation of plant will have a significant impact on control performance, communication and control became closely related and thus cannot be considered independently. In this paper, we are committed to developing a distributed collaborative estimation and control strategy over IWSANs. The main contributions are summarized as follows.

• (1)

We use interacting dual model (IDM) adaptive estimation algorithm to identify online the unreliable communication channels within IWSANs. To the best of our knowledge, the problem formulation is novel. Further, a federated estimation fusion algorithm is presented to improve the performance of accuracy and fault tolerance of state estimation.

(2)

A distributed collaborative scheme is presented which is similar to [20] but is a direct coordination method among actuators without going through the sensors.

(3)

A hybrid two-tier IWSAN based on OCARI technology [1] is developed. Based on the network topology architecture, the presented distributed estimation and control methods in (1) and (2) can be readily “piggy-backed” into industrial wireless networks such as WIA-PA [29] and OCARI [30].

The rest of this paper is organized as follows. The control system model is introduced in Section 2. A new IDM estimation and federated fusion algorithm are stated in Section 3. We design the distributed collaborative control strategy in Section 4. A simulation example of temperature control is provided in Section 5. Finally, Section 6 concludes this paper.

Notations. ${ℝ}^n$ is the n-dimensional Euclidean space. ${\mathbf{I}}_n$ denotes the $n\hspace{0.17em}\times \hspace{0.17em}n$ identity matrix. $\mathbf{P}\mathbf{r}\{·\}$ , $\mathbf{p}(·)$ , and $𝔼\{·\}$ denote the operators for probability, probability density, and estimation, respectively. $\left|·\right|$ denotes cardinality of a set. $A\to B$ denotes the map of set A to B. $\mathbf{T}\mathbf{r}(·)$ denotes the trace of a matrix. $\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(·)$ denotes a diagonal matrix.

2. Problem Formulation

In many industry fields such as printing, textile, refrigeration, and semiconductor, industrial workplaces put forwards strict requirements for some distinct parameters such as temperature, humidity, air cleanliness, pressure balance, and ventilation [5, 7]. In this paper, we consider an example of temperature supervision and control system in an industrial workshop. Based on OCARI technology, we place hybrid two-tier IWSANs organized in a cluster of cells in a workshop as shown in Figure 1. Each cell has a star topology and consists of an actuator (cell coordinator) with its group of sensors (end device nodes). The actuator is in charge of coordinating its sensors within cell and routing data packets to gateway. Furthermore, every actuator can exchange information with its neighbor actuators via the upper layer actuator network with mesh structure, so actuators can take decisions and collaboratively achieve a global target upon the industrial environment. In addition, the cooperation mechanism among actuators can improve the reliability of system. For example, when an actuator malfunctions, its neighbor actuators can be complements or alternate devices. In this paper, the actuator, also regarded as a processing center, can perform a multisource estimate fusion with its “radiofixed” intracell sensors through only single-hop communication channels. The single-hop star topology can overcome the high packet loss and long-time delays induced by multihop transmissions. The design of “radiofixed,” where the sensors need not choose their transmission paths to actuator, decreases the communication complexity and improves the real-time and reliability of data transmission.

Figure 1

Topology of hybrid two-tier IWSANs.

Our object is to control the temperature within every cell domain in a workshop to meet set value, respectively, according to the actual production requirements.

2.1. Graph of Hybrid Two-Tier IWSANs

In order to facilitate the following work, the communication relationship between the actuators and sensors and among actuators themselves over the hybrid two-tier IWSANs can be described by a graph 𝒢. Assume that there are n cells in a workshop. An actuator (i.e., central air conditioning) and a group of m temperature sensors with wireless communication function are deployed within every cell. At first, we focus on the graph of actuator network ${𝒢}_a\triangleq \{𝒜,ℰ\}$ (i.e., a subgraph of 𝒢), where $𝒜=\{a_{\mathrm{1}},\dots ,a_n\}$ is the set of n actuators and

$$\begin{array}{cc}& ℰ=\left\{(a_i,a_l)\mid a_i,a_l\in 𝒜,\hspace{0.17em}\hspace{0.17em}{a}_l\in {𝒩}_{a_i},\forall i,l\in \{\mathrm{1},\dots ,n\}\right\}\end{array}$$

(1)

is the edge set. ${𝒩}_{a_i}$ is a set of the neighbor actuators of $a_i$ . Each actuator is only allowed to communicate with neighbor actuators over actuator network, that is, $\exists (a_i,a_l)\in ℰ$ , if and only if $a_l\in {𝒩}_{a_i}$ .

Consider that $a_i$ can communicate with a group of m intracell sensors (i.e., ${𝒮}_i=\{s_i^{\mathrm{1}},\dots ,s_i^m\}$ ) over sensor network. Thus, the overall graph 𝒢 can be obtained via taking the subgraph ${𝒢}_a$ and adding $m\times n$ new vertices 𝒮, where $𝒮=S_{\mathrm{1}}\cup S_{\mathrm{2}},\dots ,\cup \hspace{0.17em}\hspace{0.17em}{S}_n$ corresponds to the set of all sensors. Define the edge set

$$\begin{array}{l}{ℰ}_{𝒪}=\{\left(s_i^j,a_i\right)\mid a_i\in 𝒜,s_i^j\in {𝒮}_i,\forall i\in \{\mathrm{1},\dots ,n\},\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\forall j\in \left\{\mathrm{1},\dots ,m\right\}\phantom{\left(a_i,s_i^j\right)}\}.\end{array}$$

(2)

Then, we obtain $𝒢=\{𝒜\cup 𝒮,ℰ\cup {ℰ}_{𝒪}\}$ .

Remark 1.

In this paper, ${ℰ}_{𝒪}$ represents the radioconnectivity between actuators and their intracell sensors. However, ℰ represents not only the radio connectivity but also interactive acting among neighboring actuators through some ventilation pipes.

2.2. Modeling

Consider the above-mentioned temperature control system, where the temperature is controlled in an industrial workshop using hybrid two-tier IWSANs as presented in Figure 1. We focus on the following discrete-time plant model:

$$\begin{array}{cc}& x_{i,k+\mathrm{1}}=g_i{x}_{i,k}+w_{ii}{u}_{i,k}^a+{\displaystyle \sum }_{a_l\in {𝒩}_{a_i}}‍w_{il}{u}_{l,k}^a,\\ y_{i,k}=h_i{x}_{i,k},\end{array}$$

(3)

where $x_{i,k}$ , $i\in \{\mathrm{1},\dots ,n\}$ is the temperature state of ith cell domain of plant at time k and is measured by a group of m intracell sensors, that is, ${𝒮}_i=\{s_i^{\mathrm{1}},\dots ,s_i^m\}$ , $y_{i,k}={[y_{i,k}^{\mathrm{1}},\dots ,y_{i,k}^m]}^T$ is the generalized measurement value of $x_{i,k}$ collected by ${𝒮}_i$ and can be sent to actuator $a_i$ , and $u_{i,k}^a$ and $u_{l,k}^a$ , $l\in \{\mathrm{1},\dots ,n\}$ , are the input (output wind temperature) corresponding to the acting applied to the state $x_{i,k}$ by local actuator $a_i$ and neighbor actuator $a_l$ , respectively.

Due to the unreliable wireless communication channels and the interference from the harsh industrial environment, it is necessary to take into account a more “realistic” plant model. First, consider the unreliable wireless communication channels within actuator network. Let I.I.D Bernoulli stochastic variable ${\gamma }_{il,k}$ with $𝔼\{{\gamma }_{il,k}\}={\gamma }_{il}$ indicate whether data packet (including state information) from $a_i$ is received ( ${\gamma }_{il,k}=\mathrm{1}$ ) or lost ( ${\gamma }_{il,k}=\mathrm{0}$ ) by neighbor $a_l\in {𝒩}_{a_i}$ . Apparently, ${\gamma }_{ii,k}\equiv \mathrm{1}$ , for all $i\in \{\mathrm{1},\dots ,n\}$ . After any neighbor $a_l$ of $a_i$ obtains the state information of $a_i$ , its control unit can compute an appropriate input to $a_i$ . Therefore, when $a_l\in {𝒩}_{a_i}$ and ${\gamma }_{il,k}=\mathrm{1}$ , $x_{i,k}$ will be actuated cooperatively by local actuator $a_i$ and its neighbor actuators $\{a_l\}$ . For instance, in Figure 1, except local $a_{\mathrm{1}}$ , cell 1 domain can also receive the other inputs sent from $a_{\mathrm{2}}$ , $a_{\mathrm{3}}$ , and $a_{\mathrm{4}}$ through some ventilation pipes, respectively.

For a concise derivation, according to the zero-input strategy in the literature [31], we can write

$$\begin{array}{cc}& u_{i,k}^a=u_{i,k},\\ u_{l,k}^a=u_{l,k},\end{array}$$

(4)

where $u_{i,k}$ and $u_{l,k}$ are the desired control laws computed by embedded control units of $a_i$ and $a_l$ , respectively.

Secondly, we use the Gilbert-Elliott model to simulate the unreliable wireless channels within sensor network. At time k, the stochastic variable ${\theta }_{i,k}$ is used to describe the packet loss model ${\theta }_{i,k}^{\mathrm{0}}$ and the packet arrival model ${\theta }_{i,k}^{\mathrm{1}}$ of communication channel from ${𝒮}_i$ to $a_i$ , for all $i\in \{\mathrm{1},\dots ,n\}$ , where ${\theta }_{i,k}^{\alpha }\triangleq \{{\theta }_{i,k}=\alpha \}$ , $\alpha =\mathrm{0,1}$ . To capture the temporal correlation of channel variation, ${\theta }_{i,k}$ can be modeled as a two-state Markovian chain with the initial model probability vector ${\pi }_{i,\mathrm{0}}$ and the transition probability matrix (TPM) ${\mathrm{\Xi }}_i$ as

$$\begin{array}{cc}& {\pi }_{i,\mathrm{0}}=\left[{\pi }_{i,\mathrm{0}}^{\mathrm{0}}\hspace{0.5em}{\pi }_{i,\mathrm{0}}^{\mathrm{1}}\right],\end{array}$$

(5)

$$\begin{array}{cc}& {\mathrm{\Xi }}_i=\left[\begin{array}{cc}{\xi }_i^{\mathrm{00}}& {\xi }_i^{\mathrm{01}}\\ {\xi }_i^{\mathrm{10}}& {\xi }_i^{\mathrm{11}}\end{array}\right],\end{array}$$

(6)

where ${\pi }_{i,\mathrm{0}}^{\alpha }\triangleq \mathbf{P}\mathbf{r}\{{\theta }_{i,\mathrm{0}}^{\alpha }\}\hspace{0.17em}\hspace{0.17em}(\alpha =\mathrm{0,1})$ is the initial model probability, ${\xi }_i^{\alpha \beta }\triangleq \mathbf{P}\mathbf{r}\{{\theta }_{i,k}^{\beta }\mid {\theta }_{i,k-\mathrm{1}}^{\alpha }\}\hspace{0.17em}\hspace{0.17em}(\alpha ,\beta =\mathrm{0,1})$ is the transition probability of switching from model ${\theta }_{i,k-\mathrm{1}}^{\alpha }$ to model ${\theta }_{i,k}^{\beta }$ , and ${\xi }_i^{\mathrm{00}}$ , ${\xi }_i^{\mathrm{11}}$ , ${\xi }_i^{\mathrm{01}}=\mathrm{1}-{\xi }_i^{\mathrm{00}}$ , and ${\xi }_i^{\mathrm{10}}=\mathrm{1}-{\xi }_i^{\mathrm{11}}$ are called packet remaining loss rate, packet remaining arrival rate, packet recovery rate, and packet failure rate, respectively. Apparently, ${\mathrm{\Xi }}_i$ , the TPM, is characterized by the ${\xi }_i^{\mathrm{00}}$ and ${\xi }_i^{\mathrm{11}}$ .

When aggregating all states at time k into the $x_k={[x_{\mathrm{1},k}^T,\dots ,x_{n,k}^T]}^T$ and considering packet losses, noises, and (4), the above plant model (3) can be written in vector form as

$$\begin{array}{cc}& x_{k+\mathrm{1}}=G{x}_k+W\left({\gamma }_k\right)u_k+{\omega }_k,\\ y_k=H\left({\theta }_k\right)x_k+{\nu }_k,\end{array}$$

(7)

where $x_k\in {ℝ}^n$ , $u_k\in {ℝ}^n$ , $y_k\in {ℝ}^{nm}$ , ${\omega }_k\in {ℝ}^n$ , ${\nu }_k\in {ℝ}^{nm}$ , ${\omega }_k$ and ${\nu }_k$ are mutually uncorrelated Gaussian, white, and zero-mean noises with covariance $Q>\mathrm{0}$ and covariance $R>\mathrm{0}$ , respectively, $G=\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(g_{\mathrm{1}},\dots ,g_n)$ , $W({\gamma }_k)\in {ℝ}^{n\times n}$ , $w_{il}({\gamma }_{il,k})=\mathrm{0}$ if $a_l\notin {𝒩}_{a_i}$ or ${\gamma }_{il,k}=\mathrm{0}$ , and $H({\theta }_k)=\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}(h_{\mathrm{1}}({\theta }_k),\dots ,h_n({\theta }_k))$ with $h_i({\theta }_k)=[h_i^{\mathrm{1}}({\theta }_{i,k}),\dots ,h_i^m({\theta }_{i,k}){]}^T$ .

In this paper, the main objective is to control temperature within every cell domain to a required set value. For the aforementioned stochastic system (7), firstly, let us give the following quadratic cost function:

$$\begin{array}{cc}& J=\frac{\mathrm{1}}{n}𝔼\left\{{\left[{\psi }_k-x_k\right)}^T\left({\psi }_k-x_k\right]\right\},\end{array}$$

(8)

where ${\psi }_k={[{\psi }_{\mathrm{1},k},\dots ,{\psi }_{n,k}]}^T$ is the temperature set value vector. Then, by estimating the temperature state, users wish to find an admissible optimal control strategy $u_k^{*}$ to minimize the cost function J as

$$\begin{array}{cc}& u_k^{*}\triangleq {𝔼}_{y^k}\left\{\underset{u_k}{\text{argmin}}\hspace{0.17em}J\right\},\end{array}$$

(9)

where $u_k\in u^k$ , $u^k=\{u_{\tau }\mid \tau =\mathrm{1},\dots ,k\}$ , and $y^k=\{y_{\tau }\mid \tau =\mathrm{1},\dots ,k\}$ .

3. Estimation and Fusion

In this section, we present a new adaptive estimation and fusion algorithm based on IDM and federated filter (IDM-FF). The architecture of the IDM-FF is shown in Figure 2. The IDM-FF method can be divided into two phases: first, the single-sensor processing phase, where temperature of every cell domain is measured by m sensors, so the actuator provides m LFs (local filters) to preprocess the data received from the m sensors, respectively, that is, every LF runs an adaptive IDM filtering algorithm to identify online the unreliable communication channel in order to estimate the measurement data of sensors, and, second, the multisensor data fusion phase, where the MF (main filter) is only to fuse the data obtained from its m preprocessing LFs so as to improve the accuracy and fault tolerance against the harsh environments and electromagnetic interference issues in industrial workshop.

Figure 2

Architecture of the IDM-FF.

3.1. IDM Adaptive Filtering

As presented in Section 2.2, the unreliable wireless channels within sensor network can be modeled as a two-state Markovian chain. In many practical situations, however, the prior knowledge of TPM (6) about the channel may be inadequate and lacking especially for large-scale wireless networks. Therefore, we propose a new adaptive channel-aware method based on IDM filter to not only identify the channel status but also learn online a posterior estimation of TPM from collected data.

First, we give the adaptive learning process online to obtain the TPM without the prior knowledge; then the updated TPM can be stored and utilized in IDM filtering algorithm by LF for channel decoding. According to the Bayesian framework, the estimation ${\overline{\mathrm{\Xi }}}_i$ of ${\mathrm{\Xi }}_i$ , model probability ${\pi }_{i,k}$ , and model likelihood ${\mathrm{\Lambda }}_{i,k}$ are defined, respectively, as follows:

$$\begin{array}{l}{\overline{\mathrm{\Xi }}}_{i,k}\triangleq 𝔼\left[{\mathrm{\Xi }}_i\mid {\left(y_i^j\right)}^k\right],\\ {\pi }_{i,k}={\left[{\pi }_{i,k}^{\mathrm{0}},{\pi }_{i,k}^{\mathrm{1}}\right]}^T\hspace{0.17em}\hspace{0.17em}\text{with}\\ {\pi }_{i,k}^{\beta }\triangleq \mathbf{P}\mathbf{r}\left\{{\theta }_{i,k}^{\beta }\mid {\overline{\mathrm{\Xi }}}_{i,k-\mathrm{1}},{\left(y_i^j\right)}^{k-\mathrm{1}}\right\},\hspace{1em}\beta =\mathrm{0,1},\\ \begin{array}{l}{\mathrm{\Lambda }}_{i,k}={\left[{\mathrm{\Lambda }}_{i,k}^{\mathrm{0}},{\mathrm{\Lambda }}_{i,k}^{\mathrm{1}}\right]}^T\hspace{0.17em}\hspace{0.17em}\text{with}\end{array}\\ {\mathrm{\Lambda }}_{i,k}^{\beta }\triangleq \mathbf{p}\left[y_{i,k}^j\mid {\theta }_{i,k}^{\beta },{\overline{\mathrm{\Xi }}}_{i,k-\mathrm{1}},{\left(y_i^j\right)}^{k-\mathrm{1}}\right],\hspace{1em}\beta =\mathrm{0,1},\end{array}$$

(10)

where $i\in \{\mathrm{1},\dots ,n\}$ , ${(y_i^j)}^k=\{y_{i,\tau }^j\mid \tau =\mathrm{1},\dots ,k\}$ .

At time k, we consider an adaptive recursive process to obtain the TPM according to posterior model probability ${\pi }_{i,k-\mathrm{1}}$ and likelihood ${\mathrm{\Lambda }}_{i,k}$ .

By the whole probability formula, one can obtain

$$\begin{array}{l}\mathbf{p}[y_{i,k}^j\mid {\mathrm{\Xi }}_i,{\left(y_i^j\right)}^{k-\mathrm{1}}]\\ \hspace{1em}\hspace{1em}={\displaystyle \sum }_{\beta =\mathrm{0}}^{\mathrm{1}}‍\mathbf{p}[y_{i,k}^j\mid {\theta }_{i,k}^{\beta },{\mathrm{\Xi }}_i,{\left(y_i^j\right)}^{k-\mathrm{1}}]\mathbf{P}\mathbf{r}\{{\theta }_{i,k}^{\beta }\mid {\mathrm{\Xi }}_i,{\left(y_i^j\right)}^{k-\mathrm{1}}\}\\ \hspace{1em}\hspace{1em}={\displaystyle \sum }_{\beta =\mathrm{0}}^{\mathrm{1}}‍\mathbf{p}[y_{i,k}^j\mid {\theta }_{i,k}^{\beta },{\mathrm{\Xi }}_i,{\left(y_i^j\right)}^{k-\mathrm{1}}]\\ \hspace{1em}\hspace{1em}\hspace{1em}\times {\displaystyle \sum }_{\alpha =\mathrm{0}}^{\mathrm{1}}‍\mathbf{P}\mathbf{r}\{{\theta }_{i,k}^{\beta }\mid {\theta }_{i,k-\mathrm{1}}^{\alpha },{\mathrm{\Xi }}_i,{\left(y_i^j\right)}^{k-\mathrm{1}}\}\\ \hspace{1em}\hspace{1em}\hspace{1em}\times \mathbf{P}\mathbf{r}\{{\theta }_{i,k-\mathrm{1}}^{\alpha }\mid {\mathrm{\Xi }}_i,{\left(y_i^j\right)}^{k-\mathrm{1}}\}\\ \hspace{1em}\hspace{1em}\approx {\displaystyle \sum }_{\beta =\mathrm{0}}^{\mathrm{1}}‍{\mathrm{\Lambda }}_{i,k}^{\beta }{\displaystyle \sum }_{\alpha =\mathrm{0}}^{\mathrm{1}}‍{\xi }_i^{\alpha \beta }{\pi }_{i,k-\mathrm{1}}^{\alpha }\\ \hspace{1em}\hspace{1em}={\pi }_{i,k-\mathrm{1}}^T{\mathrm{\Xi }}_i{\mathrm{\Lambda }}_{i,k},\end{array}$$

(11)

where the following approximations

$$\begin{array}{cc}& \mathbf{p}[y_{i,k}^j\mid {\theta }_{i,k}^{\beta },{\mathrm{\Xi }}_i,{\left(y_i^j\right)}^{k-\mathrm{1}}]\approx {\mathrm{\Lambda }}_{i,k}^{\beta },\\ \mathbf{P}\mathbf{r}\{{\theta }_{i,k-\mathrm{1}}^{\alpha }\mid {\mathrm{\Xi }}_i,{\left(y_i^j\right)}^{k-\mathrm{1}}\}\approx {\pi }_{i,k-\mathrm{1}}^{\alpha }\end{array}$$

(12)

are used in (11). These local linearity approximations are equivalent to replacing the unknown ${\mathrm{\Xi }}_i$ in (11) with its best estimate ${\overline{\mathrm{\Xi }}}_{i,k-\mathrm{1}}$ for the current recursive cycle $(k-\mathrm{1})\to k$ . This approximation approach can provide high accuracy as stated in [32].

According to (11), $\mathbf{p}[y_{i,k}^j\mid {(y_i^j)}^{k-\mathrm{1}}]\mathrm{\hspace{0.17em}\hspace{0.17em}}$ becomes

$$\begin{array}{l}\mathbf{p}\left[y_{i,k}^j\mid {\left(y_i^j\right)}^{k-\mathrm{1}}\right]\\ \hspace{1em}\hspace{1em}={\displaystyle \int }\mathbf{p}[y_{i,k}^j\mid {\mathrm{\Xi }}_i,{\left(y_i^j\right)}^{k-\mathrm{1}}]\mathbf{p}\left[{\mathrm{\Xi }}_i\mid {\left(y_i^j\right)}^{k-\mathrm{1}}\right]d{\mathrm{\Xi }}_i\\ \hspace{1em}\hspace{1em}={\displaystyle \int }{\pi }_{i,k-\mathrm{1}}^T{\mathrm{\Xi }}_i{\mathrm{\Lambda }}_{i,k}\mathbf{p}\left[{\mathrm{\Xi }}_i\mid {\left(y_i^j\right)}^{k-\mathrm{1}}\right]d{\mathrm{\Xi }}_i\\ \hspace{1em}\hspace{1em}={\pi }_{i,k-\mathrm{1}}^T{\overline{\mathrm{\Xi }}}_{i,k-\mathrm{1}}{\mathrm{\Lambda }}_{i,k}.\end{array}$$

(13)

Further, the approximate posterior probability density function (PDF) of TPM will be updated as follows:

$$\begin{array}{cc}& \mathbf{p}\left[{\mathrm{\Xi }}_i\mid {\left(y_i^j\right)}^k\right]=\frac{\mathbf{p}[y_{i,k}^j\mid {\mathrm{\Xi }}_i,{\left(y_i^j\right)}^{k-\mathrm{1}}]}{\mathbf{p}\left[y_{i,k}^j\mid {\left(y_i^j\right)}^{k-\mathrm{1}}\right]}\mathbf{p}\left[{\mathrm{\Xi }}_i\mid {\left(y_i^j\right)}^{k-\mathrm{1}}\right].\end{array}$$

(14)

Substituting (11) and (13) into (14), we can obtain

$$\begin{array}{cc}& \mathbf{p}\left[{\mathrm{\Xi }}_i\mid {\left(y_i^j\right)}^k\right]=\frac{{\pi }_{i,k-\mathrm{1}}^T{\mathrm{\Xi }}_i{\mathrm{\Lambda }}_{i,k}}{{\pi }_{i,k-\mathrm{1}}^T{\overline{\mathrm{\Xi }}}_{i,k-\mathrm{1}}{\mathrm{\Lambda }}_{i,k}}\mathbf{p}\left[{\mathrm{\Xi }}_i\mid {\left(y_i^j\right)}^{k-\mathrm{1}}\right].\end{array}$$

(15)

Finally, we compute the estimation of TPM by a numerical integration as

$$\begin{array}{l}{\overline{\mathrm{\Xi }}}_{i,k}={\displaystyle \int }{\mathrm{\Xi }}_i\mathbf{p}\left[{\mathrm{\Xi }}_i\mid {\left(y_i^j\right)}^k\right]d{\mathrm{\Xi }}_i\\ =\frac{\mathrm{1}}{N}{\displaystyle \sum }_{s=\mathrm{1}}^N‍{\mathrm{\Xi }}_i^s\mathbf{p}\left[{\mathrm{\Xi }}_i^s\mid {\left(y_i^j\right)}^k\right],\end{array}$$

(16)

which means the posterior PDF distribution in N grid TPMs ${\mathrm{\Xi }}_i^s$ points over $[\mathrm{0,1}]$ .

Next, we use the IDM filter algorithm with ${\overline{\mathrm{\Xi }}}_{i,k-\mathrm{1}}$ and ${\pi }_{i,k-\mathrm{1}}$ to estimate the measurement data. The basic idea of the presented IDM filtering approach requires running two same parallel optimal filters for the two possible channel model sequences (i.e., ${\theta }_{i,k}^{\mathrm{0}}$ and ${\theta }_{i,k}^{\mathrm{1}}$ ), and then the measurement data can be calculated using a combined estimate. We define the local mean and variance of $x_{i,k}^j$ measured and sent by sensor $s_i^j$ with the model ${\theta }_{i,k}^{\alpha }$ , $\alpha =\mathrm{0,1}$ , at time k as follows:

$$\begin{array}{cc}& {\widehat{x}}_{i,k\mid k}^{j,\alpha }\triangleq 𝔼\left[x_{i,k}^{j,\alpha }\mid {\left(y_i^j\right)}^k\right],\\ P_{i,k\mid k}^{j,\alpha }\triangleq 𝔼\left[\left(x_{i,k}^{j,\alpha }-{\widehat{x}}_{i,k\mid k}^{j,\alpha }\right){\left(x_{i,k}^{j,\alpha }-{\widehat{x}}_{i,k\mid k}^{j,\alpha }\right)}^T\mid {\left(y_i^j\right)}^k\right].\end{array}$$

(17)

For brevity, we provide a summary of IDM filter algorithm as follows.

Step 1 (interacting for two models of varying channel).

The mixed probabilities ${\xi }_{i,k}^{\beta \mid \alpha }$ for channel models ${\theta }_{i,k}^{\alpha }$ and ${\theta }_{i,k}^{\beta }$ at time k are calculated as

$$\begin{array}{cc}& {\xi }_{i,k}^{\beta \mid \alpha }=\frac{{\overline{\xi }}_{i,k-\mathrm{1}}^{\beta \alpha }{\pi }_{i,k-\mathrm{1}}^{\beta }}{{{\displaystyle \sum }}_{\beta =\mathrm{0}}^{\mathrm{1}}‍{\overline{\xi }}_{i,k-\mathrm{1}}^{\beta \alpha }{\pi }_{i,k-\mathrm{1}}^{\beta }},\end{array}$$

(18)

where ${\pi }_{i,k-\mathrm{1}}^{\beta }$ is probabilities for model ${\theta }_{i,k}^{\beta }$ at time $k-\mathrm{1}$ . ${\overline{\xi }}_{k-\mathrm{1}}^{\beta \alpha }$ can be calculated and updated using (16) and is the recursively Bayesian posterior estimation of transition probabilities. Then, the mixed inputs (mean and covariance) for every filter can be computed as

$$\begin{array}{l}{\tilde{x}}_{i,k-\mathrm{1}}^{j,\alpha }={\displaystyle \sum }_{\beta =\mathrm{0}}^{\mathrm{1}}‍{\xi }_{i,k}^{\beta \mid \alpha }{\widehat{x}}_{i,k-\mathrm{1}\mid k-\mathrm{1}}^{j,\beta },\\ {\tilde{P}}_{i,k-\mathrm{1}}^{j,\alpha }={\displaystyle \sum }_{\beta =\mathrm{0}}^{\mathrm{1}}‍{\xi }_{i,k}^{\beta \mid \alpha }\\ \hspace{1em}\times \{\phantom{{\left[{\widehat{x}}_{i,k-\mathrm{1}\mid k-\mathrm{1}}^{j,\beta }-{\tilde{x}}_{i,k-\mathrm{1}}^{j,\alpha }\right]}^T}{P}_{i,k-\mathrm{1}\mid k-\mathrm{1}}^{j,\beta }+\left[{\widehat{x}}_{i,k-\mathrm{1}\mid k-\mathrm{1}}^{j,\beta }-{\tilde{x}}_{i,k-\mathrm{1}}^{j,\alpha }\right]\\ \phantom{mm}\hspace{1em}\times {\left[{\widehat{x}}_{i,k-\mathrm{1}\mid k-\mathrm{1}}^{j,\beta }-{\tilde{x}}_{i,k-\mathrm{1}}^{j,\alpha }\right]}^T\},\end{array}$$

(19)

where ${\widehat{x}}_{i,k-\mathrm{1}\mid k-\mathrm{1}}^{j,\beta }$ and ${\widehat{P}}_{i,k-\mathrm{1}\mid k-\mathrm{1}}^{j,\beta }$ are the updated state mean and covariance for model ${\theta }_{i,k}^{\beta }$ , respectively, at time $k-\mathrm{1}$ .

Step 2 (information Kalman filler (IKF)).

For each model ${\theta }_{i,k}^{\alpha }$ , the IKF is done as

• $\hspace{0.17em}\text{IK}{\text{F}}_u$

:

$$\begin{array}{l}{\widehat{x}}_{i,k\mid k}^{j,\alpha }=P_{i,k\mid k}^{j,\alpha }{\left(P_{i,k\mid k-\mathrm{1}}^{j,\alpha }\right)}^{-\mathrm{1}}{\widehat{x}}_{i,k\mid k-\mathrm{1}}^{j,\alpha }\\ \hspace{1em}+P_{i,k\mid k}^{j,\alpha }{\left(h_i^j\left({\theta }_{i,k}^{\alpha }\right)\right)}^T{\left(R_i^j\right)}^{-\mathrm{1}}{y}_{i,k}^{j,\alpha },\\ {\left(P_{i,k\mid k}^{j,\alpha }\right)}^{-\mathrm{1}}={\left(P_{i,k\mid k-\mathrm{1}}^{j,\alpha }\right)}^{-\mathrm{1}}\\ \hspace{1em}+{\left(h_i^j\left({\theta }_{i,k}^{\alpha }\right)\right)}^T{\left(R_i^j\right)}^{-\mathrm{1}}{h}_i^j\left({\theta }_{i,k}^{\alpha }\right),\end{array}$$

(20)

 

$\text{IK}{\text{F}}_p$ :

$$\begin{array}{l}{\widehat{x}}_{i,k\mid k-\mathrm{1}}^{j,\alpha }=g_{ii}{\tilde{x}}_{i,k-\mathrm{1}}^{j,\alpha }+w_{ii}{u}_{i,k-\mathrm{1}}\\ +{\displaystyle \sum }_{a_l\in {𝒩}_{a_i}}‍w_{il}\left({\gamma }_{il,k-\mathrm{1}}\right)u_{l,k-\mathrm{1}},\\ {\left(P_{i,k\mid k-\mathrm{1}}^{j,\alpha }\right)}^{-\mathrm{1}}={\left(g_{ii}{\tilde{P}}_{i,k-\mathrm{1}}^{j,\alpha }{g}_{ii}^T+Q_i^j\right)}^{-\mathrm{1}},\end{array}$$

(21)

 

where the prediction and update steps of IKF are denoted by $\text{IK}{\text{F}}_p$ and $\text{IK}{\text{F}}_u$ , respectively, $Q_i^j={\beta }_i^j{Q}_i$ , according to the rule of information-distribution technology of federated filter stated in Section 3.2.

Step 3 (model probabilities update).

By Bayes’ rule, the probabilities of each model ${\theta }_{i,k}^{\alpha }$ at time k are updated as

$$\begin{array}{cc}& {\pi }_{i,k}^{\alpha }=\frac{\left({{\displaystyle \sum }}_{\beta =\mathrm{0}}^{\mathrm{1}}‍{\overline{\xi }}_{i,k-\mathrm{1}}^{\beta \alpha }{\xi }_{i,k-\mathrm{1}}^{\beta }\right){\mathrm{\Lambda }}_{i,k}^{\alpha }}{{{\displaystyle \sum }}_{\alpha =\mathrm{0}}^{\mathrm{1}}‍\left({{\displaystyle \sum }}_{\beta =\mathrm{0}}^{\mathrm{1}}‍\left({\overline{\xi }}_{i,k-\mathrm{1}}^{\beta \alpha }{\xi }_{i,k-\mathrm{1}}^{\beta }\right){\mathrm{\Lambda }}_{i,k}^{\alpha }\right)}.\end{array}$$

(22)

Step 4 (combined estimation).

The combined estimate for the state mean and covariance at time k are computed as

$$\begin{array}{l}{\widehat{x}}_{i,k\mid k}^j={\displaystyle \sum }_{\alpha =\mathrm{0}}^{\mathrm{1}}‍{\pi }_{i,k}^{\alpha }{\widehat{x}}_{i,k\mid k}^{j,\alpha },\\ P_{i,k\mid k}^j={\displaystyle \sum }_{\alpha =\mathrm{0}}^{\mathrm{1}}‍{\pi }_{i,k}^{\alpha }\times \{\phantom{{\left[{\widehat{x}}_{i,k\mid k}^{j,\alpha }-{\widehat{x}}_{i,k\mid k}^j\right]}^T}{P}_{i,k\mid k}^{i,\alpha }+\left[{\widehat{x}}_{i,k\mid k}^{j,\alpha }-{\widehat{x}}_{i,k\mid k}^j\right]\\ \phantom{\text{mmmm}}\hspace{1em}\times {\left[{\widehat{x}}_{i,k\mid k}^{j,\alpha }-{\widehat{x}}_{i,k\mid k}^j\right]}^T\}.\end{array}$$

(23)

From the above results, we can know that the purpose of using IDM filter in every LF is to estimate the probability of channel model switching between packet loss model ${\theta }_{i,k}^{\mathrm{0}}$ and packet arrival model ${\theta }_{i,k}^{\mathrm{1}}$ to reduce the channel uncertainty and then obtain the combined estimation of state mean and covariance.

3.2. Federated Multisensor Estimation Fusion

In order to enhance the accuracy and fault tolerance of the whole measurement system in harsh industrial environment, the federated filter (FF) [33, 34], a decentralized filtering technology, is adopted to fuse the data of LFs (i.e., ${\widehat{x}}_{i,k\mid k}^j$ , $P_{i,k\mid k}^j$ , $j=\mathrm{1},\dots ,m$ ) for every $x_{i,k}$ in MF.

To remove the correlation among the states in different LFs, according to the rule of information distribution, we have

$$\begin{array}{cc}& Q_i^j={\beta }_i^j{Q}_i,\end{array}$$

(24)

where ${\beta }_i^j>\mathrm{0}$ is the information-distribution coefficient of the jth LF and is subjected to

$$\begin{array}{cc}& {\displaystyle \sum }_{j=\mathrm{1}}^m‍{\beta }_i^j=\mathrm{1}.\end{array}$$

(25)

Because $(h_i^j{)}^T(R_i^j{)}^{-\mathrm{1}}{h}_i^j$ can reflect the measurement accuracy of local measurement, ${\sigma }_i^j=(h_i^j{)}^T(R_i^j{)}^{-\mathrm{1}}{h}_i^j$ is defined as an information variable. Thus, the optimal information-distribution method [35] is given as

$$\begin{array}{cc}& {\beta }_i^j=\frac{\mathbf{T}\mathbf{r}\left({\sigma }_i^j\right)}{{{\displaystyle \sum }}_{j=\mathrm{1}}^m‍\mathbf{T}\mathbf{r}\left({\sigma }_i^j\right)}.\end{array}$$

(26)

Then, the global state estimation ${\widehat{x}}_{i,k\mid k}^g$ and its covariance $P_{i,k\mid k}^g$ obtained from the MF can be given as follows:

$$\begin{array}{l}{P}_{i,k\mid k}^g={\left[{\left(P_{i,k\mid k}^{\mathrm{1}}\right)}^{-\mathrm{1}}+{\left(P_{i,k\mid k}^{\mathrm{2}}\right)}^{-\mathrm{1}}+\cdots +{\left(P_{i,k\mid k}^m\right)}^{-\mathrm{1}}\right]}^{-\mathrm{1}},\\ {\widehat{x}}_{i,k\mid k}^g=P_{i,k\mid k}^g[{\left(P_{i,k\mid k}^{\mathrm{1}}\right)}^{-\mathrm{1}}{\widehat{x}}_{i,k\mid k}^{\mathrm{1}}+{\left(P_{i,k\mid k}^{\mathrm{2}}\right)}^{-\mathrm{1}}{\widehat{x}}_{i,k\mid k}^{\mathrm{2}}\\ \phantom{\text{mm}}\hspace{1em}+\cdots +{\left(P_{i,k\mid k}^m\right)}^{-\mathrm{1}}{\widehat{x}}_{i,k\mid k}^m],\end{array}$$

(27)

which indicates that the federated global estimation fusion ${\widehat{x}}_{i,k\mid k}^g$ for $x_{i,k}$ is a linear weighted combination of the local estimates ${\widehat{x}}_{i,k\mid k}^j$ with weighting matrices $P_{i,k\mid k}^g{(P_{i,k\mid k}^j)}^{-\mathrm{1}}$ , $j=\mathrm{1},\dots ,m$ . To a better fault tolerance of sensor measurements, we select a no-reset mode of federated filter which will not cross infect between LFs, as shown in Figure 2.

In summary, the main motivation for the above-mentioned results stems from the need to obtain reliable measurement data from sensors over wireless sensor network in the harsh industrial environment. In the following, we consider the cooperative control problem among actuators.

4. Distributed Control

In every cell domain within workshop, actuator can decide its control input according to the estimated local temperature state provided by its IDM-FF presented in Section 3. However, the stability of the whole system cannot be guaranteed due to overlapping actuation between actuators. Fortunately, the upper layer actuator network of two-tier IWSANs provides a foundation for coordinating the control inputs of actuators to improve the performance of control system. This implies that every actuator can obtain more complete state information of plant by exchanging information with neighbor actuators over the actuator network. Therefore, in this paper, we consider that every actuator can firstly perform a locally optimized control law without accounting for the effect of neighbor actuators, and then a cooperation mechanism is adopted to eliminate the overshoot caused by the interaction with neighbor actuators. Comparing the cooperation mechanism with the one presented in [20], the main difference between them is that any actuator can directly communicate with its neighbor actuators in our method, unlike the indirect communication between actuators going through the sensors in [20].

In this section, a locally optimal control method is considered firstly for actuator $a_i$ without collaborations with neighbor actuators: first, every $a_i$ can exchange information with neighbor $a_l$ if and only if $\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}{\gamma }_k=\mathrm{1}$ via the upper actuator network. Second, every $a_i$ can decide its control input to its own ith cell domain and all neighbor l cell domains, respectively, based on the received information packets at the first step. Next, we describe the information flow among actuators by ${𝒢}_a$ . The input matrix $W({\gamma }_k)$ in (7) can be degenerated as a binary matrix $W^{\prime }({\gamma }_k)$ to indicate the connectivity of ${𝒢}_a$ by replacing any nonzero elements with one. Thus, we can use $W^{\prime }({\gamma }_k)=[w_{il}^{\prime }({\gamma }_{il,k})]\in {ℝ}^{n\times n}$ as the adjacency matrix of ${𝒢}_a$ :

$$\begin{array}{cc}& w_{il}^{\prime }\left({\gamma }_{il,k}\right)=\{\begin{array}{cc}1\hfill & \text{if}\hspace{0.17em}\hspace{0.17em}{w}_{il}\left({\gamma }_{il,k}\right)\ne \mathrm{0}\hfill \\ 0\hfill & \text{otherwrise}.\hfill \end{array}\end{array}$$

(28)

Remark 2.

Consider that $a_i$ is a neighborhood itself; that is, $w_{ii}^{\prime }({\gamma }_k)\equiv \mathrm{1}$ in (28).

For $a_i$ , we define the number of its neighbors ${𝒩}_{a_i}$ as

$$\begin{array}{cc}& {\rho }_{i,k}\triangleq {\left|{𝒩}_{a_i}\right|}_{{\gamma }_k=\mathrm{1}}={\displaystyle \sum }_{a_l\in {𝒩}_{a_i}}‍w_{il}^{\prime }\left({\gamma }_k=\mathrm{1}\right),\end{array}$$

(29)

where $\mathrm{1}⩽{\rho }_{i,k}⩽\left|{𝒩}_{a_i}\right|⩽n$ .

According to the aforementioned distributed collaborative control mechanism, every actuator needs to obtain the information of neighbor actuators before deciding a control law. So, $a_i$ can use a row selection matrix ${\mathrm{\Gamma }}_{i,k}\in {ℝ}^{{\rho }_{i,k}\times n}$ to collect information from ${\rho }_{i,k}$ neighbor actuators. ${\mathrm{\Gamma }}_{i,k}$ is defined as

$$\begin{array}{cc}& {\mathrm{\Gamma }}_{i,k}={\left[e_l^T\right]}_{l\mid w_{il}^{\prime }\left({\gamma }_k\right)=\mathrm{1}},\hspace{1em}l=\mathrm{1},\dots ,n,\end{array}$$

(30)

where $e_l$ denotes the lth vector of the standard basis of ${ℝ}^n$ . Apparently, ${\mathrm{\Gamma }}_{i,k}^T$ is a column selection matrix. Thus, actuator $a_i$ can receive available information from neighbor actuators (including itself) as

$$\begin{array}{cc}& {\mathrm{\Psi }}_{i,k}={\mathrm{\Gamma }}_{i,k}{\psi }_k,\end{array}$$

(31)

$$\begin{array}{cc}& X_{i,k+\mathrm{1}}=G_{i,k}{X}_{i,k}+{\overline{W}}_k{u}_k+{\mathrm{\Omega }}_{i,k},\end{array}$$

(32)

where $X_{i,k}={\mathrm{\Gamma }}_{i,k}{x}_k$ , ${\mathrm{\Omega }}_{i,k}={\mathrm{\Gamma }}_{i,k}{\omega }_k$ with covariance $Q_i^{\prime }>\mathrm{0}$ , $G_{i,k}={\mathrm{\Gamma }}_{i,k}G{\mathrm{\Gamma }}_{i,k}^T$ , and ${\overline{W}}_k={\mathrm{\Gamma }}_{i,k}W({\gamma }_k)$ .

Remark 3.

Note that the process of selecting neighbor actuators of $a_i$ via ${\mathrm{\Gamma }}_{i,k}$ is a dynamic process based on ${\gamma }_k$ due to the uncertainty of wireless actuator network. For example, $a_l\in {𝒩}_{a_i}$ is a neighbor actuator of $a_i$ if ${\gamma }_{k-\mathrm{1}}=\mathrm{1}$ (packet reception) at time $k-\mathrm{1}$ but is not a neighbor actuator of $a_i$ if ${\gamma }_k=\mathrm{0}$ (packet loss) at time k.

From the aforementioned results, we can obtain the cost function only computed by $a_i$ as

$$\begin{array}{cc}& J_{i,k}=\frac{\mathrm{1}}{{\rho }_{i,k}}𝔼\left\{{\left({\mathrm{\Psi }}_{i,k}-X_{i,k}\right)}^T\left({\mathrm{\Psi }}_{i,k}-X_{i,k}\right)\right\}.\end{array}$$

(33)

Theorem 4.

At time $k-\mathrm{1}$ , one can find the locally optimal control law $u_{i,k}^{*}$ of actuator $a_i$ without accounting for the effects of the neighbor actuators ${𝒩}_{a_i}$ as follows:

$$\begin{array}{l}{u}_{i,k}^{*}=𝔼\left\{\underset{u_{i,k}}{\text{argmin}}{J}_{i,k}\right\}\\ ={\left({\overline{W}}_{i,k}^T{\overline{W}}_{i,k}\right)}^{-\mathrm{1}}{\overline{W}}_{i,k}^T\left({\mathrm{\Psi }}_{i,k+\mathrm{1}}-G_{i,k}{\widehat{X}}_{i,k\mid k}\right),\end{array}$$

(34)

where ${\widehat{X}}_{i,k\mid k}={\mathrm{\Gamma }}_{i,k}[{\widehat{x}}_{\mathrm{1},k\mid k}^g,\dots ,{\widehat{x}}_{n,k\mid k}^g{]}^T$ , and ${\overline{W}}_{i,k}$ is the ith column of matrix ${\overline{W}}_k$ .

Proof.

According to (33), one can obtain

$$\begin{array}{l}{J}_{i,k+\mathrm{1}}=\frac{\mathrm{1}}{{\rho }_{i,k+\mathrm{1}}}𝔼\left\{{\left({\mathrm{\Psi }}_{i,k+\mathrm{1}}-X_{i,k+\mathrm{1}}\right)}^T\left({\mathrm{\Psi }}_{i,k+\mathrm{1}}-X_{i,k+\mathrm{1}}\right)\right\}\\ =\frac{\mathrm{1}}{{\rho }_{i,k+\mathrm{1}}}(\phantom{{\mathrm{\Psi }}_{i,k+\mathrm{1}}^T}𝔼\left\{-\mathrm{2}{\mathrm{\Psi }}_{i,k+\mathrm{1}}^T{X}_{i,k+\mathrm{1}}+X_{i,k+\mathrm{1}}^T{X}_{i,k+\mathrm{1}}\right\}\\ \hspace{1em}\phantom{\text{mm}}\hspace{1em}+{\mathrm{\Psi }}_{i,k+\mathrm{1}}^T{\mathrm{\Psi }}_{i,k+\mathrm{1}}).\end{array}$$

(35)

Substituting (32) into (35), we have

$$\begin{array}{l}{J}_{i,k+\mathrm{1}}\\ \hspace{1em}=\frac{\mathrm{1}}{{\rho }_{i,k+\mathrm{1}}}(𝔼\{\phantom{{\overline{W}}_k^T}-\mathrm{2}{\mathrm{\Psi }}_{i,k+\mathrm{1}}^T(G_{i,k}{X}_{i,k}+{\overline{W}}_k{u}_k+{\mathrm{\Omega }}_{i,k})\\ \hspace{1em}\hspace{1em}\hspace{1em}\phantom{\text{mmm}}\hspace{1em}+{\left(G_{i,k}{X}_{i,k}+{\overline{W}}_k{u}_k+{\mathrm{\Omega }}_{i,k}\right)}^T\\ \hspace{1em}\hspace{1em}\hspace{1em}\phantom{\text{mmm}}\hspace{1em}\times (G_{i,k}{X}_{i,k}+{\overline{W}}_k{u}_k+{\mathrm{\Omega }}_{i,k})\phantom{{\overline{W}}_k^T}\}+{\mathrm{\Psi }}_{i,k+\mathrm{1}}^T{\mathrm{\Psi }}_{i,k+\mathrm{1}})\\ \hspace{1em}=\frac{\mathrm{1}}{{\rho }_{i,k+\mathrm{1}}}(𝔼\{u_k^T{\overline{W}}_k^T{\overline{W}}_k{u}_k+\mathrm{2}{u}_k^T{\overline{W}}_k^T\\ \hspace{1em}\hspace{1em}\hspace{1em}\phantom{\text{mmn}}\hspace{1em}\times (G_{i,k}{X}_{i,k}-{\mathrm{\Psi }}_{i,k+\mathrm{1}})\phantom{{\overline{W}}_k^T}\}+{\mathrm{\Psi }}_{i,k+\mathrm{1}}^T{\mathrm{\Psi }}_{i,k+\mathrm{1}}+\mathbf{T}\mathbf{r}\left(Q_i^{\prime }\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\phantom{\text{m}}\hspace{1em}+𝔼\left\{X_{i,k}^T{G}_{i,k}^T{G}_{i,k}{X}_{i,k}-\mathrm{2}{X}_{i,k}^T{G}_{i,k}^T{\mathrm{\Psi }}_{i,k+\mathrm{1}}\phantom{{\overline{W}}_k^T}\right\}).\end{array}$$

(36)

The minimum of $J_{i,k+\mathrm{1}}$ in (36) controlled by $u_k$ is

$$\begin{array}{l}{V}_{i,k+\mathrm{1}}=\frac{\mathrm{1}}{{\rho }_{i,k+\mathrm{1}}}(\underset{u_k}{\min }𝔼\{u_k^T{\overline{W}}_k^T{\overline{W}}_k{u}_k\\ \hspace{1em}\hspace{1em}\phantom{mmmmm}\hspace{1em}+\mathrm{2}{u}_k^T{\overline{W}}_k^T(G_{i,k}{X}_{i,k}-{\mathrm{\Psi }}_{i,k+\mathrm{1}})\}\phantom{{\overline{W}}_k^T}).\end{array}$$

(37)

By solving $\partial V_{i,k+\mathrm{1}}/\partial u_k=\mathrm{0}$ , that is,

$$\begin{array}{cc}& \frac{\partial V_{i,k+\mathrm{1}}}{\partial u_k}=\mathrm{2}{\overline{W}}_k^T{\overline{W}}_k{u}_k+\mathrm{2}{\overline{W}}_k^T(G_{i,k}{\widehat{X}}_{i,k\mid k}-{\mathrm{\Psi }}_{i,k+\mathrm{1}})=\mathrm{0},\end{array}$$

(38)

where ${\widehat{X}}_{i,k\mid k}={\mathrm{\Gamma }}_{i,k}{[{\widehat{x}}_{\mathrm{1},k\mid k}^g,\dots ,{\widehat{x}}_{n,k\mid k}^g]}^T$ , we can obtain the optimal control law $u_{i,k}^{*}$ as

$$\begin{array}{cc}& u_{i,k}^{*}={\left({\overline{W}}_k^T{\overline{W}}_k\right)}^{-\mathrm{1}}{\overline{W}}_k^T\left({\mathrm{\Psi }}_{i,k+\mathrm{1}}-G_{i,k}{\widehat{X}}_{i,k\mid k}\right).\end{array}$$

(39)

Since $u_{\mathrm{1},k},u_{\mathrm{2},k},\dots ,u_{n,k}$ are independent, the locally optimal control law $u_{i,k}^{*}$ of $a_i$ without accounting for the effects of the neighbor actuators ${𝒩}_{a_i}$ can be expressed in (34). This completes the proof.

The locally optimal control law $u_{i,k}^{*}$ stated in Theorem 4 does not consider the issues of cooperation among neighbor actuators for the stability of the entire system. Therefore, in order to eliminate the overlapping action between neighbor actuators, we adopt a distributed cooperation mechanism for every actuator as follows:

$$\begin{array}{cc}& u_{i,k}^{†}={\mu }_{i,k}{u}_{i,k}^{*},\hspace{1em}i=\mathrm{1},\dots ,n,\end{array}$$

(40)

where ${\mu }_{i,k}$ is a tuning parameter. Next, we introduce how to compute ${\mu }_{i,k}$ for $a_i$ . First, at time k, every actuator $a_l$ , $a_l\in {𝒩}_{a_i}$ , can use an information set ${ℐ}_{l,k}=\{w_{lj}({\gamma }_{lj,k})u_{j,k}^{*}\mid a_j\in {𝒩}_{a_l},{\gamma }_{lj,k}=\mathrm{1}\}$ , where $u_{j,k}^{*}$ is calculated by (34). Second, every $a_l$ computes the number of elements in set ${ℐ}_{l,k}$ that have the same sign as one another and then obtains a factor set $\mathrm{\Lambda }=\{{\lambda }_{jl,k}\}$ . For $a_l$ , ${\lambda }_{jl,k}$ means the number of its neighbor actuators that can produce overlapping actuation with $a_l$ at time k. Then, after receiving $\mathrm{\Lambda }$ , we can obtain the tuning parameter for $a_i$ as

$$\begin{array}{cc}& {\mu }_{i,k}=\frac{\mathrm{1}}{{\rho }_{i,k}}{\displaystyle \sum }_{a_l\in N_{a_i}}‍\frac{\mathrm{1}}{{\lambda }_{il,k}}.\end{array}$$

(41)

Note that the aforementioned cooperation mechanism in (41) is similar to [20]; however, our method runs a direct cooperation among actuators without going through the sensors.

5. Numerical Example

In this section, we consider an example of temperature control in industrial environment over the hybrid two-tier IWSANs. The simulation experiment is built on MATLAB software platform. The monitored workshop field is 100 m × 100 m and deployed 7 cells. The actuating radius of actuators within each cell is $r_a=\mathrm{17}$  m, and the detection radius of each sensor is $r_s=\mathrm{8}$  m. Assume that the environment temperature of every cell domain can be controlled by a central air conditioning (i.e., an actuator with estimate/control unit) and collected by 5 intracell sensors. According to the actual production requirements, assume that our object is to control the temperature within every cell domain to meet set value ${\mathrm{\Psi }}_i$ , $i\in \{\mathrm{1},\dots ,\mathrm{7}\}$ , chosen from $[{\mathrm{10}}^{°}\text{C},{\mathrm{40}}^{°}\text{C}]$ . Let the plant model be specified as follows. Consider $g_{ii}=\mathrm{0.9}$ , for all i, and weight coefficient $w_{il}({\gamma }_{il,k})$ is distributed randomly in $[\mathrm{0,1}]$ . The average loss packet rate of channel is $𝔼({\gamma }_{il,k})=\mathrm{75}\%$ in the actuator network. For the sensor network, assume that ${\theta }_{\mathrm{1},k}={\theta }_{\mathrm{2},k}=\cdots ={\theta }_{\mathrm{7},k}={\theta }_k$ , $h_i^j({\theta }_k^{\alpha })=\alpha$ , $\alpha =\mathrm{0,1}$ , the initial channel model probability is ${\pi }_{\mathrm{0}}=[\mathrm{0.5}\hspace{0.5em}\mathrm{0.5}]$ , the true TPM is chosen with ${\xi }^{\mathrm{00}}=\mathrm{0.3}$ and ${\xi }^{\mathrm{11}}=\mathrm{0.9}$ , and the IDM uses a numerical integration estimator with $\mathrm{50}$ grid points over $[\mathrm{0,1}]$ for posterior TPM.

In simulation running, let the number of the total simulation steps and the measurement packets from sensors successfully arriving at actuator be Ϝ and $\mathrm{\Delta }$ (i.e., the number of model ${\theta }_k^{\mathrm{1}}$ ), respectively. According to the TPM in (6), the times channel state ${\theta }_k$ switching from ${\theta }_k^{\mathrm{1}}$ to ${\theta }_k^{\mathrm{0}}$ is about $\mathrm{\Delta }{\xi }^{\mathrm{10}}$ . We introduce a variable $T_{\text{loss}}\in \mathbb{Z}$ that denotes the times of ${\theta }_k$ remaining at ${\theta }_k^{\mathrm{0}}$ . Therefore, the average times ${\theta }_k$ remaining at ${\theta }_k^{\mathrm{0}}$ can be estimated as

$$\begin{array}{cc}& 𝔼\left\{T_{\text{loss}}\right\}={\displaystyle \sum }_{t=\mathrm{1}}^{\infty }‍t{\left({\xi }^{\mathrm{00}}\right)}^{t-\mathrm{1}}=\frac{\mathrm{1}}{{\left(\mathrm{1}-{\xi }^{\mathrm{00}}\right)}^{\mathrm{2}}}=\frac{\mathrm{1}}{{\left({\xi }^{\mathrm{01}}\right)}^{\mathrm{2}}}.\end{array}$$

(42)

The total times ${\theta }_k^{\mathrm{0}}$ in sequence of ${\theta }_k$ can be described approximately as $\mathrm{\Delta }{\xi }^{\mathrm{10}}/{({\xi }^{\mathrm{01}})}^{\mathrm{2}}$ ; thus

$$\begin{array}{cc}& \mathrm{\Delta }+\mathrm{\Delta }\frac{{\xi }^{\mathrm{10}}}{{\left({\xi }^{\mathrm{01}}\right)}^{\mathrm{2}}}\approx Ϝ.\end{array}$$

(43)

Finally, the above equation can be rewritten as

$$\begin{array}{cc}& \mathrm{\Delta }\approx \frac{{\left({\xi }^{\mathrm{01}}\right)}^{\mathrm{2}}}{{\xi }^{\mathrm{10}}+{\left({\xi }^{\mathrm{01}}\right)}^{\mathrm{2}}}Ϝ.\end{array}$$

(44)

If we take $Ϝ=\mathrm{100}$ , then $\mathrm{\Delta }\approx \mathrm{83}$ , $\mathrm{\Delta }{\xi }^{\mathrm{10}}\approx \mathrm{9}$ , and the average arrival rates of channel state ${\theta }_k$ are about $\mathrm{83}\%$ . We report the approximate distribution of $T_{\text{loss}}$ in Table 1. A sample sequence of ${\theta }_k$ is generated as shown in Figure 3.

Table 1

Analysis of packet losses.

Distribution of $T_{\text{loss}}$
$T_{\text{loss}}$	$\mathbf{P}\mathbf{r}\{T_{\text{loss}}=t\}={\xi }^{\mathrm{01}}({\xi }^{\mathrm{00}}{)}^{t-\mathrm{1}}$	Times of $T_{\text{loss}}=t$
1	0.7000	5
2	0.2100	2
3	0.0630	1
4	0.0189	1

Figure 3

Example of sample sequence of ${\theta }_k$ .

The results of simulation are shown in Figures 4–8. Figure 4 illustrates the control performance of the system. It is clearly seen that the system with cooperation is stable, but the cost function curve without cooperation has an about $\mathrm{20}\%$ deviation due to overlapping effects among neighbor actuators. Figure 5 illustrates the convergence of the TPM estimators. Figure 6 shows that the adaptive IDM has a better accuracy of state estimate than the nonadaptive one with a mismatched TPM ( ${\xi }^{\mathrm{00}}={\xi }^{\mathrm{11}}=\mathrm{0.6}$ ). Because the prior knowledge of true TPM for wireless communication channel with random packet losses is unknown in most real situations, it is preferable to use an adaptive algorithm to estimate the TPM. Figure 7 compares the accuracy of mean trace for estimation covariance for LF ( $\mathrm{Tr}(p^l)$ ), SF ( $\mathrm{Tr}(p^s)$ ), and MF ( $\mathrm{Tr}(p^g)$ ). Note that SF (single filter) indicates this case in which the temperature state $x_i$ is measured by single sensor rather than a multisensor measurement of federated estimation fusion. It shows that, first, the accuracy of MF is improved much compared with LF, and, second, federated estimators perform better than estimators without fusion. Figure 8 shows that the adaptive IDM filtering can provide a better true model probability. As one can expect, the distributed collaborative control approach can ameliorate effectively the control performance of the whole system, and the adaptive IDM-FF shows substantially better overall performance than the nonadaptive one.

Figure 4

Control performance.

Figure 5

Convergence of TPM estimators.

Figure 6

Mean absolute estimate error.

Figure 7

Mean trace for estimate covariance.

Figure 8

IDM probability of the true model.

6. Conclusions

In this paper, we consider the problem of distributed control and estimate over hybrid two-tier IWASNs for industrial automation control applications. In order to identify the unreliable communication channel, a novel IDM with the adaptive channel-aware algorithm is designed. The federated estimation fusion algorithm is presented to further improve the performance of accuracy and fault tolerance of state estimation. A distributed control scheme with coordination is proposed. The system stability and effectiveness of the presented methods are shown by simulation results.

Future work includes finding a multi-index (such as temperature, humidity, and air cleanliness) fusion and control strategy for production environments in a workshop and considering a coordination mechanism to minimize the average energy consumption of system.