首页 关于本刊 编 委 会 期刊动态 作者中心 审者中心 读者中心 下载中心 联系我们 English
 自动化学报  2017, Vol. 43 Issue (7): 1234-1240 PDF

1. 中国科学技术大学自动化系 合肥 230027;
2. 火灾科学国家重点实验室 中国科学技术大学先进技术研究院合肥 230027;
3. 中国科学技术大学先进技术研究院 合肥 230027;
4. 中国科学院空间信息处理与应用系统技术重点实验室 北京 100190;
5. 安庆师范大学物理与电子信息学院 安庆 246011

Optimization Algorithms for Predictive Control Approach to Networked Bilinear Systems
Binglin Wang1, Yu Kang1,2,3,4, Jiahu Qin1, Yanmei Li5
1. Department of Automation, University of Science and Technology of China, Hefei 230027, China;
2. State Key Laboratory of Fire Science, Department of Automation, Institute of Advanced Technology, University of Science and Technology of China, Hefei 230027, China, and also with the Key Laboratory of Technology in Geo-Spatial Information Processing and Application System, Chinese Academy of Sciences, Beijing 100190, China;
3. Physics and Electronic Engineering, Anqing Normal University, Anqing 246011, China
Abstract: This paper is concerned with the networked predictive control of discrete-time bilinear systems. To deal with the network-induced communication delay that exists in both forward channel (controller to actuator) and feedback channel (sensor to controller), a bilinear networked predictive control scheme is proposed. Then a non-convex optimization problem of solving the predictive control sequence is presented, for which two gradually-optimized algorithms are proposed based on the special structure of bilinear system dynamics model. The numerical simulation indicates that the resulting predictive control sequence can compensate for the network-induced issues actively, which proves the effectiveness of the proposed predictive control strategy.
Key words: Bilinear systems     networked control systems (NCSs)     non-convex optimization     predictive control
1 Introduction

Networked control systems (NCSs) are spatially distributed control systems, where controllers, sensors and actuators are connected via a limited bandwidth communication network. There are many advantages of NCSs, such as flexible architectures, low installation costs and high reliability. However, the insertion of the communication network will inevitably lead to communication delay, data packets dropout and disorder [1], which usually result in performance degradation or even instability of systems [2], [3].

Focusing on these issues, various approaches have been developed, for example, stochastic control approach [4], switched system approach [5], queuing approach [6], robust control approach [7] and so on. However, these methods are proposed based on some restrictive assumptions on NCSs which ignore the features of NCSs. In recent years, a new approach called the networked predictive control is proposed by taking full advantage of the specialities of NCSs, such as timestamp technique, packet-based transmission and caching technology [8], [9]. As shown in Fig. 1, the networked predictive control system mainly consists of control prediction generator (CPG) and network delay compensator (NDC), where a set of future control predictions are generated by CPG and the appropriate control value in the latest control prediction sequence is chosen to compensate the network-induced issues actively by NDC. In addition, two data buffers are located at the controller and actuator to reorder and store the received data from networks. Based on the networked predictive control scheme, many results have been obtained for networked linear control systems. The design, analysis and implementation of the networked predictive control scheme have been widely studied [10]-[12]. The stability and robustness criteria have been derived when uncertainties exist in NCSs [13], [14].

 Figure 1 Networked predictive control system.

Following the predictive control based approach of linear systems in [8], we aim to study the networked predictive control for bilinear systems in this paper. As a kind of important nonlinear system, the bilinear systems are applied widely in engineering, economics, biology and other areas [15]. Further, the bilinear systems can approach other nonlinear systems with higher precision than any other traditional linear systems [16]. Therefore, it is of both theoretical and practical interest to explore the networked predictive control of bilinear systems. Actually, the networked predictive control scheme has been extended to some nonlinear systems [17]-[21]. However, the analytical methods are restricted to some specific nonlinear systems, for instance, hammerstein systems [17] and wiener systems [18], which are not applicable to bilinear systems. The networked predictive control schemes for general nonlinear systems are proposed in [19]-[21]. However, the linearization theory of nonlinear systems adopted in [19] is always accompanied by strong conservatism. The system states and control variables in [20] are subject to some constraints, which will reduce the versatility of system. The nonlinear dynamical controller is given in [21] in advance, which is always difficult to be obtained in reality.

Different from most of the existing literatures which state the form of controller firstly, a constrained optimization problem consisted of system states and inputs is considered in this paper to obtain the future control inputs. However, the system predictive states will present a highly coupled form because of the nonlinearity inherent in bilinear model, which makes the optimization problem to be non-convex. A general method for solving such an optimization problem is the sequential quadratic programming (SQP) [22], where a quadratic program is solved to determine the search direction and step size. However, if the approximation of SQP is poor, the search direction derived would also be poor. Another effective method is linearizing the nonlinear model about the input trajectory [23], [24]. However, the linearization method always impose strong constraints which subsequently leads to strong conservatism. Different from approximating the original problem to a quadratic program by SQP or linearization, the special structure of bilinear model will be fully exploited in this paper, and the original non-convex optimization problem will be transformed into a quadratic program subproblem. Based on the above analysis, two algorithms, called the stepwise iterative (SI) algorithm and approximate forward stagewise (AFS) algorithm, will be proposed to calculate the predictive control sequence of the networked bilinear systems. Subsequently, the future inputs will be transmitted to the actuator via the forward channel.

The remainder of this paper is organized as follows: In Section 2, some assumptions and problem formulations are introduced. Further, the networked predictive control scheme is implemented and the related issues are converted to a non-convex optimization problem. Two gradually-optimized algorithms are given to solve the proposed optimization problem in Section 3. The numerical simulation is presented and the conclusions are summarized in Sections 4 and 5, respectively.

2 Design of Networked Predictive Control System for Networked Bilinear Systems

Consider the following discrete-time nonlinear system

 ${x_{t + 1}} = A{x_t} + {u_t}N{x_t}$ (1)

where $x_t\in \mathbb{R}^n$ is the system state vector, $u_t\in \mathbb{R}$ is a scalar control input. The $x_t$ and $u_t$ multiplied together which leads the system (1) to be bilinear system. Meanwhile, A $\in$ $\mathbb{R}^{n\times n}$, $N\in \mathbb{R}^{n\times n}$ are constant system matrices.

To simplify the presentation of the networked bilinear predictive control scheme, some assumptions are made as follows.

Assumption 1: The communication delay $i_t$ in the forward channel is bounded by $M_d$ and the communication delay $k_t$ in the feedback channel is bounded by $N_d$, i.e., $i_t$ $\leq$ $M_d$ and $k_t \leq N_d$.

Assumption 2: The number of consecutive data dropouts in the forward channel and feedback channel are upper bounded by $M_l$ and $N_l$, respectively.

Assumption 3: The sensor, controller and actuator considered in the networked predictive control system (shown in Fig. 1) are time-synchronized.

Assumption 4: All the data packets sent from both the sensor and the controller are time-stamped to notify when they are sent.

Remark 1: From Assumptions 1 and 2, we can know that the sum of consecutive data dropouts and communication delay in the forward channel and feedback channel are all upper bounded, i.e., $M_d+M_l\leq \overline{M}$ and $N_d+N_l\leq \overline{N}$, where $\overline{M}$ and $\overline{N}$ are positive integers. Moreover, according to Assumption 3, $\overline{M}$ should not be lager than the length of the control prediction sequence, which is sent from the controller to the actuator to compensate for the network--induced issues actively.

Remark 2: In practical NCSs, the so called data dropout can be divided into three situations: 1) the packet fails to arrive at the destination in a certain transmission period, it will be treated as dropout based on the network protocols although it arrives in the future; 2) the "first sent later arrival", i.e., packet A is sent earlier than packet B, but it arrives later than packet B, then packet A would be treated as dropout; 3) the real data loss, i.e., the packet never reached. Additionally, the number of consecutive data dropouts should be finite in order to avoid the system becoming open loop. Thus, the Assumption 2 is reasonable.

Remark 3: From Assumptions 3 and 4, the network-induced delays in the feedback channel and forward channel can be determined, respectively, according to the time when the data packets arrive at the controller and actuator.

Remark 4: In the networked predictive control scheme, when the data dropouts happen in the forward channel within the current sampling period, the $(i+1)$th value of the predictive control sequence that received in the previous sampling period will be applied to the plant in the case of forward delay i. In the same way, when the data are lost in the feedback channel, the CPG will use the latest system state which received in the controller. Therefore, the data dropouts can be treated as "long delay" in the proposed compensation scheme.

Based on the aforementioned assumptions, the following networked predictive control method for networked bilinear systems is proposed.

When the communication delay in the feedback channel is $k_t$, $k_t \in \{0, 1, 2, \ldots, \overline{N}\}$, according to the system model (1), the one--step ahead state prediction can be designed by

 $\begin{array}{l} {{\hat x}_{t - {k_t} + 1|t - {k_t}}} = A{x_{t - {k_t}}} + {u_{t - {k_t}}}N{x_{t - {k_t}}}\\ \quad \quad \quad \quad = \left[ {A + {u_{t - {k_t}}}N} \right]{x_{t - {k_t}}}. \end{array}$ (2)

Following by (2), the state predictions from time $t-k_t+2$ to time t can be constructed by

 $\begin{array}{l} {{\hat x}_{t - {k_t} + 2|t - {k_t}}} = A{{\hat x}_{t - {k_t} + 1|t - {k_t}}} + {u_{t - {k_t} + 1}}N{{\hat x}_{t - {k_t} + 1|t - {k_t}}}\\ \quad \quad \quad \quad = \left[ {A + {u_{t - {k_t} + 1}}N} \right]{{\hat x}_{t - {k_t} + 1|t - {k_t}}}\\ \quad \quad \quad \quad \quad \vdots \\ \quad \quad {{\hat x}_{t|t - {k_t}}} = A{{\hat x}_{t - 1|t - {k_t}}} + {u_{t - 1}}N{{\hat x}_{t - 1|t - {k_t}}}\\ \quad \quad \quad \quad \quad = [A + {u_{t - 1}}N]{{\hat x}_{t - 1|t - {k_t}}} \end{array}$ (3)

which results in

 ${\hat x_{t|t - {k_t}}} = \prod\limits_{j = 1}^{{k_t}} {\left[ {A + {u_{t - j}}N} \right]} {x_{t - {k_t}}}.$ (4)

Moreover, when the communication delay $i_t$, $i_t \in\{0, 1, 2,$ $\ldots, \overline{M}\}$, also exists in the forward channel, the state predictions from time $t+1$ to time $t+i_t$ can be constructed as

 $\begin{array}{l} {{\hat x}_{t + 1|t - {k_t}}} = A{{\hat x}_{t|t - {k_t}}} + {u_{t|t - {k_t}}}N{{\hat x}_{t|t - {k_t}}}\\ \quad \quad \quad = \left[ {A + {u_{t|t - {k_t}}}N} \right]{{\hat x}_{t|t - {k_t}}}\\ {{\hat x}_{t + 2|t - {k_t}}} = A{{\hat x}_{t + 1|t - {k_t}}} + {u_{t + 1|t - {k_t}}}N{{\hat x}_{t + 1|t - {k_t}}}\\ \quad \quad \quad = \left[ {A + {u_{t + 1|t - {k_t}}}N} \right]{{\hat x}_{t + 1|t - {k_t}}}\\ \quad \quad \quad \quad \vdots \\ {{\hat x}_{t + {i_t}|t - {k_t}}} = A{{\hat x}_{t + {i_t} - 1|t - {k_t}}} + {u_{t + {i_t} - 1|t - {k_t}}}N{{\hat x}_{t + {i_t} - 1|t - {k_t}}}\\ \quad \quad \quad = \left[ {A + {u_{t + {i_t} - 1|t - {k_t}}}N} \right]{{\hat x}_{t + {i_t} - 1|t - {k_t}}}{\rm{ }} \end{array}$ (5)

which results in

 ${\hat x_{t + {i_t}|t - {k_t}}} = \prod\limits_{l = 1}^{{i_t}} {\left[ {A + {u_{t + {i_t} - l|t - {k_t}}}N} \right]} {\hat x_{t|t - {k_t}}}.$ (6)

Combining with (4), the predictive state of time $t+i_t$ can be obtained at time $t-k_t$

 ${\hat x_{t + {i_t}|t - {k_t}}} = \prod\limits_{l = 1}^{{i_t}} {\left[ {A + {u_{t + {i_t} - l|t - {k_t}}}N} \right]} \prod\limits_{j = 1}^{{k_t}} {\left[ {A + {u_{t - j}}N} \right]} {x_{t - {k_t}}}.$ (7)

Remark 5: From (7), it is clear that the predictive state of time $t+i_t$ consists of system state $x_{t-k_t}$, the past control input from time $t-k_t$ to time $t-1$, and the predictive control inputs from time t to time $t+i_t-1$. It is noteworthy that all these items are known for time t, so predictive state (7) is theoretically feasible.

Considering the optimization problems depend on state vector and control vector

 \begin{align}\label{eq10} & \min\ J(k)=\overline{X}^TQ_x\overline{X}+\overline{U}^TR_u\overline{U}\notag\\ & \begin{cases} x_{t+1}=Ax_t+u_tNx_t &\\ \| u_{t+i}\| \le \theta_i,&i\in \{0,1,2,\ldots,\overline{M}-1\}\\ \end{cases} \end{align} (8)

where $\overline{X}$ and $\overline{U}$ represent the state vector and control vector, respectively, $Q_x$ and $R_u$ are constant positive definite matrixes and $\|u_{t+i-1}\| \le \theta_i$ represents the input constraints.

Obviously, (8) is a non-convex optimization problem because of the nonlinear form of system (1) and predictive state (7). Therefore, setting the derivative of performance index with respect to the control vector be zero directly did not obtain the desired control values. It is the main difficulty in dealing with networked predictive control of nonlinear system, which prompts us to explore new algorithms.

3 Two Algorithms for Solving the Non-convex Optimization Problem

Observe the system predictive state (7), we can find that the order of each control input in the expression is up to one order. It means that the system state can be represented as the first order polynomial of one control input when the other inputs are known. This characteristic is decided by the special structure of bilinear system (1), which can be used as the breakthrough point to solve the non-convex optimization problem (8).

Based on the above observation, two algorithms are proposed as follows. Algorithm 1, which is called the stepwise iterative algorithm, begins with a given initial sequence of control inputs. Only one control input is calculated every time and the results obtained previously will be used to calculate the latter ones. Substituting the obtained control inputs into performance index, and iterating the whole process until the variations of index meets the given threshold value $\lambda$. Algorithm 2, which is called the approximate forward stagewise algorithm, deals with the one-dimensional (scalar) quadratic optimization, firstly. Then the dimensions of control input vector and state vector are added step by step until all the predictive control inputs are derived. Every time a new control input is calculated, the old inputs obtained before are updated, thus the approximate values will be approaching the optimal solution gradually.

3.1 Stepwise Iterative Algorithm

An initial sequence of control inputs is given firstly

 ${{\bar{U}}_{0}}={{\left[ {{u}_{t|t-{{k}_{t}}}},{{u}_{t+1|t-{{k}_{t}}}},\ldots ,{{u}_{t+\bar{M}|t-{{k}_{t}}}} \right]}^{T}}.$

Now, assume that the communication delay in the feedback channel is $k_t$, $k_t\in\{0, 1, 2, \ldots, \overline{N}\}$, and the communication delay in the forward is $i_t$, $i_t\in\{0, 1, 2, \ldots, \overline{M}\}$.

Let

 ${{\bar{X}}_{\bar{N},\bar{M}}}={{\left[ \hat{x}_{t+1|t-{{k}_{t}}}^{T},\hat{x}_{t+2|t-{{k}_{t}}}^{T},\ldots ,\hat{x}_{t+\bar{M}+1|t-{{k}_{t}}}^{T} \right]}^{T}}.$ (9)

According to (7), we have known that the state prediction $\hat{x}_{t+i_t|t-k_t}$, $i_t=1, 2, \ldots, \overline{M}+1$, is constant or the first order polynomial of the single control input. Thus, considering the future time $t+i$, the state $\hat{x}_{t+i_t|t-k_t}$ can be expressed by $u_{t+i|t-k_t}$ by using the given control inputs in $\overline{U}_0$, i.e.,

 \begin{align} &{{{\hat{x}}}_{t+{{i}_{t}}|t-{{k}_{t}}}}={{K}_{{{i}_{t}},i}}+{{V}_{{{i}_{t}},i}}{{u}_{t+i|t-{{k}_{t}}}} \\ &\qquad \qquad \qquad \qquad \qquad {{i}_{t}}=1,2,\ldots ,\bar{M}+1 \\ \end{align} (10)

where $K_{i_t, i}$ and $V_{i_t, i}$ are all consist of partial initial values in $\overline{U}_0$, optimal control values that have been calculated, the past control inputs and the system state $x_{t-k_t}$.

When it is time to calculate $u^{\star}_{t+i|t-k_t}$, the predictive control sequence $[u^{\star}_{t|t-k_t}, u^{\star}_{t+1|t-k_t}, \ldots, u^{\star}_{t+i-1|t-k_t}]^T$ has been obtained. Then the initial sequence $\overline{U}_0$ can be updated as the following

 \begin{align} &\bar{U}_{t+i-1}^{\star }=\ [u_{t|t-{{k}_{t}}}^{\star },\ldots ,u_{t+i-1|t-{{k}_{t}}}^{\star },{{u}_{t+i|t-{{k}_{t}}}},\ldots , \\ &\quad \quad \quad {{u}_{t+{{i}_{t}}|t-{{k}_{t}}}},\ldots ,{{u}_{t+\bar{M}|t-{{k}_{t}}}}{{]}^{T}}. \\ \end{align}

Further according to (7), the $\hat{x}_{t+i_t|t-k_t}$ can be rewritten as

 \begin{align} &{{{\hat{x}}}_{t+{{i}_{t}}|t-{{k}_{t}}}}=\prod\limits_{l=1}^{{{i}_{t}}-i-1}{\left[ A+{{u}_{t+{{i}_{t}}-l|t-{{k}_{t}}}}N \right]}\left[ A+{{u}_{t+i|t-{{k}_{t}}}}N \right] \\ &\ \quad \ \quad \ \quad \times \prod\limits_{m={{i}_{t}}-i+1}^{{{i}_{t}}}{\left[ A+u_{t+{{i}_{t}}-m|t-{{k}_{t}}}^{\star }N \right]} \\ &\ \quad \ \quad \ \quad \ \times \prod\limits_{n=1}^{{{k}_{t}}}{\left[ A+{{u}_{t-n}}N \right]}{{x}_{t-{{k}_{t}}}}. \\ \end{align} (11)

Setting $u_{t+i|t-k_t}=0$ in (11), then we can get

 \begin{align} &{{K}_{{{i}_{t}},i}}=\ {{{\hat{x}}}_{t+{{i}_{t}}|t-{{k}_{t}}}}[{{u}_{t+i|t-{{k}_{t}}}}=0] \\ &\quad \quad =\prod\limits_{l=1}^{{{i}_{t}}-i-1}{\left[ A+{{u}_{t+{{i}_{t}}-l|t-{{k}_{t}}}}N \right]}A \\ &\ \quad \quad \quad \times \prod\limits_{m={{i}_{t}}-i+1}^{{{i}_{t}}}{\left[ A+u_{t+{{i}_{t}}-m|t-{{k}_{t}}}^{\star }N \right]} \\ &\ \quad \quad \quad \times \prod\limits_{n=1}^{{{k}_{t}}}{\left[ A+{{u}_{t-n}}N \right]}{{x}_{t-{{k}_{t}}}}. \\ \end{align} (12)

And setting $u_{t+i|t-k_t}=1$ in (11), we can get

 \begin{align} &{{V}_{{{i}_{t}},i}}={{{\hat{x}}}_{t+{{i}_{t}}|t-{{k}_{t}}}}[{{u}_{t+i|t-{{k}_{t}}}}=1]-{{K}_{{{i}_{t}},i}} \\ &\quad \quad =\prod\limits_{l=1}^{{{i}_{t}}-i-1}{\left[ A+{{u}_{t+{{i}_{t}}-l|t-{{k}_{t}}}}N \right]}N \\ &\times \prod\limits_{m={{i}_{t}}-i+1}^{{{i}_{t}}}{\left[ A+u_{t+{{i}_{t}}-m|t-{{k}_{t}}}^{\star }N \right]}\prod\limits_{n=1}^{{{k}_{t}}}{\left[ A+{{u}_{t-n}}N \right]}{{x}_{t-{{k}_{t}}}}. \\ \end{align} (13)

Remark 6: It is noticed that the $K_{i_t, i}$ and $V_{i_t, i}$ in (10) are also related to the future control inputs. Therefore, when $i<i_t$, the $K_{i_t, i}$ and $V_{i_t, i}$ are calculated by using (12) and (13). When $i\geq i_t$, $\hat{x}_{t+i_t|t-k_t}$ is constant and can be expressed by the calculated control inputs, i.e.,

 ${{\hat{x}}_{t+{{i}_{t}}|t-{{k}_{t}}}}=\prod\limits_{l=1}^{{{i}_{t}}}{\left[ A+u_{t+{{i}_{t}}-l|t-{{k}_{t}}}^{\star }N \right]}\prod\limits_{j=1}^{{{k}_{t}}}{\left[ A+{{u}_{t-j}}N \right]}{{x}_{t-{{k}_{t}}}}$

and then we hav

 ${{K}_{{{i}_{t}},i}}={{\hat{x}}_{t+{{i}_{t}}|t-{{k}_{t}}}},\ {{V}_{{{i}_{t}},i}}=0,\quad \quad {{i}_{t}}=1,2,\ldots ,i.$

Based on the above discussion, the system state vector $\overline{X}_{\overline{N}, \overline{M}}$ in (9) can be rewritten as

 \begin{align} &{{{\bar{X}}}_{\bar{N},\bar{M},i}}=\left[ \begin{matrix} {{K}_{1,i}}+{{V}_{1,i}}{{u}_{t+i|t-{{k}_{t}}}} \\ \vdots \\ {{K}_{\bar{M}+1,i}}+{{V}_{\bar{M}+1,i}}{{u}_{t+i|t-{{k}_{t}}}} \\ \end{matrix} \right] \\ &\quad \quad \quad ={{K}_{i}}+{{V}_{i}}{{u}_{t+i|t-{{k}_{t}}}} \\ \end{align}

where

 ${{K}_{i}}=\left[ \begin{matrix} {{K}_{1,i}} \\ \vdots \\ {{K}_{\bar{M}+1,i}} \\ \end{matrix} \right],\ {{V}_{i}}=\left[ \begin{matrix} {{V}_{1,i}} \\ \vdots \\ {{V}_{\bar{M}+1,i}} \\ \end{matrix} \right].$

Consider the following performance index

 \begin{align} &{{J}_{\bar{N},\bar{M},i}}=\ \bar{X}_{\bar{N},\bar{M},i}^{T}{{Q}_{\bar{N},\bar{M}}}{{{\bar{X}}}_{\bar{N},\bar{M},i}}+\bar{U}_{\bar{N},\bar{M},i}^{T}{{R}_{\bar{N},\bar{M}}}{{{\bar{U}}}_{\bar{N},\bar{M},i}} \\ &\quad \quad \quad =\ {{[{{K}_{i}}+{{V}_{i}}{{u}_{t+i|t-{{k}_{t}}}}]}^{T}}{{Q}_{\bar{N},\bar{M}}}[{{K}_{i}}+{{V}_{i}}{{u}_{t+i|t-{{k}_{t}}}}] \\ &\quad \quad \quad +\bar{U}_{\bar{N},\bar{M},i}^{T}{{R}_{\bar{N},\bar{M}}}{{{\bar{U}}}_{\bar{N},\bar{M},i}} \\ \end{align} (14)

where $Q_{\overline{N}, \overline{M}}$ and $R_{\overline{N}, \overline{M}}$ are positive definite constant matrixes.

Now, let the partial derivative of $J_{\overline{N}, \overline{M}, i}$ with respect to $u_{t+i|t-k_t}$ be zero, i.e., $\frac{\partial J_{\overline{N}, \overline{M}, i}}{\partial u_{t+i|t-k_t}}=0$, we can obtain

 $u_{t+i|t-{{k}_{t}}}^{\star }=-K_{i}^{T}{{Q}_{\bar{N},\bar{M}}}{{V}_{i}}{{\left[ V_{i}^{T}{{Q}_{\bar{N},\bar{M}}}{{V}_{i}}+\sum\limits_{j=1}^{\bar{M}+1}{{{R}_{i,j}}} \right]}^{-1}}.$ (15)

Updating $\overline{U}^{\star}_{t+i-1}$ by replacing $u_{t+i|t-k_t}$ with $u^{\star}_{t+i|t-k_t}$, we will obtain $\overline{U}^{\star}_{t+i}$, i.e.,

 \begin{align} &\bar{U}_{t+i}^{\star }=\ [u_{t|t-{{k}_{t}}}^{\star },\ldots ,u_{t+i-1|t-{{k}_{t}}}^{\star }, \\ &\quad \ \quad \ u_{t+i|t-{{k}_{t}}}^{\star },{{u}_{t+i+1|t-{{k}_{t}}}},\ldots ,{{u}_{t+\bar{M}|t-{{k}_{t}}}}{{]}^{T}}. \\ \end{align}

Let $\overline{U}_0=\overline{U}^{\star}_{t+i}$ as the new predictive control initial values, and continuing the above process until $\overline{U}^{\star}_{t+\overline{M}}$ is obtained, where $\overline{U}^{\star}_{t+\overline{M}}$ is described as follows.

 $\begin{array}{*{35}{l}} \bar{U}_{t+\bar{M}}^{\star }=[u_{t|t-{{k}_{t}}}^{\star },&\ldots ,u_{t+{{i}_{t}}|t-{{k}_{t}}}^{\star },\ldots ,u_{t+\bar{M}|t-{{k}_{t}}}^{\star }{{]}^{T}}. \\ \end{array}$

After the predictive control sequence $\overline{U}^{\star}_{t+\overline{M}}$ is transmitted from the controller side to the actuator side via the forward channel, it becomes $\overline{U}^{\star}_{t+\overline{M}-i_t}$ because of the communication delay $i_t$ in the forward channel, where $\overline{U}^{\star}_{t+\overline{M}-i_t}$ can be expressed by

 \begin{align} &\bar{U}_{t+\bar{M}-{{i}_{t}}}^{\star } \\ &\quad ={{[u_{t-{{i}_{t}}|t-{{i}_{t}}-{{k}_{t}}}^{\star },\ldots ,u_{t|t-{{i}_{t}}-{{k}_{t}}}^{\star },\ldots ,u_{t+\bar{M}-{{i}_{t}}|t-{{i}_{t}}-{{k}_{t}}}^{\star }]}^{T}}. \\ \end{align}

To compensate for the communication delay $i_t$ actively, the predictive control value $u^{\star}_{t|t-i_t-k_t}$ will be applied to the system at the current time t.

The detailed steps of Algorithm 1 are given in the following.

 Algorithm 1 Stepwise iterative algorithm Step 1. Given an initial predictive control sequence $\overline{U}_0$ and a threshold $\lambda$. Step 2. Calculate $J_{\rm old}(k)$ by using $\overline{U}_0$. Step 3. Determine the communication delay $k_t$ in the feedback channel. For the forward communication delay $i=0, 1, \ldots, \overline{M}$: Step 3.1 Calculate $K_i$ according to (12); Step 3.2 Calculate $V_i$ according to (13); Step 3.3 Calculate $u^{\star}_{t+i|t-k_t}$ according to (15). Step 4. Calculate $J_{\rm new}(k)$ by using the obtained $\overline{U}^{\star}_{t+\overline{M}}$. Step 5. If $[J_{\rm old}(k)-J_{\rm new}(k)]/J_{\rm old}(k)>\lambda$, let $J_{\rm old}(k)=J_{\rm new}(k),$ $\overline{U}_0=\overline{U}^{\star}_{t+\overline{M}}$, and go to Step 3; Otherwise, stop iterating and let $\overline{U}^{\star}_{t+\overline{M}}$ as the final optimal predictive control sequence.

Each control input in the predictive control sequence calculated by Algorithm 1 will make the performance index to be minimal, which leads the final optimal predictive control sequence to approximate the optimal solution with any small error after several iterations. However, it should be noted that Algorithm 1 will cost lots of computation time if the initial sequence $\overline{U}_0$ is far away from the optimal values.

Remark 7: As the system state converge asymptotically to the equilibrium state over time, the control inputs should be reduced to avoid system oscillation. Therefore, before calculating a new control sequence for the next sampling period, some processing should be done to the initial sequence $\overline{U}_0$, such as be multiplied by a weight $\alpha$, where $0<$ $\alpha$ $<$ $1$.

3.2 Approximate Forward Stagewise Algorithm

Algorithm 1 usually requires many iterations, which result lots of computation and thus it is time consuming. For higher calculation efficiency, Algorithm 2 is proposed as follows.

For convenience, we consider that the forward communication delay is $i_t-1$, $i_t\in \{2, 3, \ldots, \overline{M}+1\}$ and the feedback communication delay is $k_t$, $k_t\in \{0, 1, 2, \ldots, \overline{N}\}$. Assume that partial control inputs have been obtained as follows

 $\bar{U}_{{{k}_{t}},{{i}_{t}}-2}^{\star }={{\left[ u_{t|t-{{k}_{t}}}^{\star },u_{t+1|t-{{k}_{t}}}^{\star },\ldots ,u_{t+{{i}_{t}}-2|t-{{k}_{t}}}^{\star } \right]}^{T}}.$

According to (6), we have

 \begin{align} &{{{\hat{x}}}_{t+{{i}_{t}}|t-{{k}_{t}}}} \\ &\quad =\left[ A+{{u}_{t+{{i}_{t}}-1|t-{{k}_{t}}}}N \right]\prod\limits_{j=2}^{{{i}_{t}}}{\left[ A+u_{t+{{i}_{t}}-j|t-{{k}_{t}}}^{\star }N \right]}{{{\hat{x}}}_{t|t-{{k}_{t}}}} \\ \end{align}

where only $u_{t+i_t-1|t-k_t}$ is unknown.

Further it is easy to find that $\hat{x}_{t+i_t|t-k_t}$ is the first order polynomial of $u_{t+i_t-1|t-k_t}$ under the condition $\overline{U}^{\star}_{k_t, i_t-2}$ is known. Therefore, $\hat{x}_{t+i_t|t-k_t}$ can be redescribed as

 ${{\hat{x}}_{t+{{i}_{t}}|t-{{k}_{t}}}}={{K}_{{{i}_{t}}}}+{{V}_{{{i}_{t}}}}{{u}_{t+{{i}_{t}}-1|t-{{k}_{t}}}}$ (16)

where $K_{i_t}$ and $V_{i_t}$ can be calculated as follows

 \begin{align} &{{K}_{{{i}_{t}}}}={{{\hat{x}}}_{t+{{i}_{t}}|t-{{k}_{t}}}}[{{u}_{t+{{i}_{t}}-1|t-{{k}_{t}}}}=0] \\ &\quad =A\prod\limits_{j=2}^{{{i}_{t}}}{\left[ A+u_{t+{{i}_{t}}-j|t-{{k}_{t}}}^{\star }N \right]}{{{\hat{x}}}_{t|t-{{k}_{t}}}} \\ \end{align} (17)
 \begin{align} &{{V}_{{{i}_{t}}}}={{{\hat{x}}}_{t+{{i}_{t}}|t-{{k}_{t}}}}[{{u}_{t+{{i}_{t}}-1|t-{{k}_{t}}}}=1]-{{K}_{{{i}_{t}}}} \\ &\quad =N\prod\limits_{j=2}^{{{i}_{t}}}{\left[ A+u_{t+{{i}_{t}}-j|t-{{k}_{t}}}^{\star }N \right]}{{{\hat{x}}}_{t|t-{{k}_{t}}}}. \\ \end{align} (18)

Obviously, $K_{i_t}$ and $V_{i_t}$ are comprised of the control inputs in $\overline{U}^{\star}_{k_t, i_t-2}$ and the predictive state $\hat{x}_{t|t-k_t}$, so they can be calculated.

Let

 \begin{align} &{{{\bar{U}}}_{{{k}_{t}},{{i}_{t}}-1}}={{\left[ {{u}_{t|t-{{k}_{t}}}},{{u}_{t+1|t-{{k}_{t}}}},\ldots ,{{u}_{t+{{i}_{t}}-1|t-{{k}_{t}}}} \right]}^{T}} \\ &{{{\bar{X}}}_{{{k}_{t}},{{i}_{t}}}}={{\left[ \hat{x}_{t+1|t-{{k}_{t}}}^{T},\hat{x}_{t+2|t-{{k}_{t}}}^{T},\ldots ,\hat{x}_{t+{{i}_{t}}|t-{{k}_{t}}}^{T} \right]}^{T}} \\ \end{align}

where $\overline{U}_{k_t, i_t-1}\in \mathbb{R}^{i_t}$ and $\overline{X}_{k_t, i_t}\in \mathbb{R}^{i_tn}$.

Based on the discussion above, we have

 ${{\bar{X}}_{{{k}_{t}},{{i}_{t}}}}={{\bar{K}}_{{{i}_{t}}}}+{{\bar{V}}_{{{i}_{t}}}}{{\bar{U}}_{{{k}_{t}},{{i}_{t}}-1}}$ (19)

where

 ${{\bar{K}}_{{{i}_{t}}}}=\left[ \begin{matrix} {{K}_{1}} \\ \vdots \\ {{K}_{{{i}_{t}}}} \\ \end{matrix} \right],\ \ {{\bar{V}}_{{{i}_{t}}}}=\left[ \begin{matrix} {{V}_{1}}&\ldots &0 \\ \vdots &\ddots &\vdots \\ 0&\ldots &{{V}_{{{i}_{t}}}} \\ \end{matrix} \right]$

Defining the following performance index

 ${{J}_{{{k}_{t}},{{i}_{t}}-1}}=\bar{X}_{{{k}_{t}},{{i}_{t}}}^{T}{{Q}_{{{k}_{t}},{{i}_{t}}}}{{\bar{X}}_{{{k}_{t}},{{i}_{t}}}}+\bar{U}_{{{k}_{t}},{{i}_{t}}-1}^{T}{{R}_{{{k}_{t}},{{i}_{t}}}}{{\bar{U}}_{{{k}_{t}},{{i}_{t}}-1}}$

where $Q_{k_t, i_t}\in \mathbb{R}^{i_tn\times i_tn}$ and $R_{k_t, i_t}\in \mathbb{R}^{i_t\times i_t}$ are positive definite matrices.

Let the partial derivative of $J_{k_t, i_t-1}$ with respect to $\overline{U}_{k_t, i_t-1}$ be zero, i.e., $\frac{\partial J_{k_t, i_t-1}}{\partial \overline{U}_{k_t, i_t-1}}=0$, we have

 $\bar{U}_{{{k}_{t}},{{i}_{t}}-1}^{\star T}=-\bar{K}_{{{i}_{t}}}^{T}{{Q}_{{{k}_{t}},{{i}_{t}}}}{{\bar{V}}_{{{i}_{t}}}}{{\left[ \bar{V}_{{{i}_{t}}}^{T}{{Q}_{{{k}_{t}},{{i}_{t}}}}{{{\bar{V}}}_{{{i}_{t}}}}+{{R}_{{{k}_{t}},{{i}_{t}}}} \right]}^{-1}}.$ (20)

Remark 8: For the case $i_t=1$, i.e., there is no delay in the forward channel, the predictive state (16) can be stated as

 ${{\hat{x}}_{t+1|t-{{k}_{t}}}}=[A+{{u}_{t|t-{{k}_{t}}}}N]{{\hat{x}}_{t|t-{{k}_{t}}}}$

and the performance index $J_{k_t, 0}$ becomes

 ${{J}_{{{k}_{t}},0}}=\hat{x}_{t+1|t-{{k}_{t}}}^{T}{{Q}_{{{k}_{t}},1}}{{\hat{x}}_{t+1|t-{{k}_{t}}}}+u_{t|t-{{k}_{t}}}^{T}{{R}_{{{k}_{t}},1}}{{u}_{t|t-{{k}_{t}}}}.$

Let the partial derivative of $J_{k_t, 0}$ with respect to $u_{t|t-k_t}$ be zero, we can obtain $u^{\star}_{t|t-k_t}$. Accordingly, the Algorithm 2 can be activated to calculate $\overline{U}^{\star}_{k_t, i_t-1}$.

Assume that the communication delay $k_t$ in the feedback channel is up to $\overline{N}$, then the procedure above should be continued until all the predictive control inputs are calculated, i.e., $\overline{U}^{\star}_{\overline{N}, \overline{M}}$ is obtained as following

 $\bar{U}_{\bar{N},\bar{M}}^{\star }={{\left[ u_{t|t-\bar{N}}^{\star },u_{t+1|t-\bar{N}}^{\star },\ldots ,u_{t+\bar{M}-1|t-\bar{N}}^{\star },u_{t+\bar{M}|t-\bar{N}}^{\star } \right]}^{T}}.$

Next, the predictive control inputs $\overline{U}^{\star}_{\overline{N}, \overline{M}}$ will be transmitted to the actuator via the forward channel.

Remark 9: After obtaining $\overline{U}^{\star}_{k_t, i_t-1}$ by the promised Algorithm 2, $\overline{U}^{\star}_{k_t, i_t}$ will be calculated next. The new predictive value $u^{\star}_{t+i_t|t-k_t}$ will be achieved and an update to $\overline{U}^{\star}_{k_t, i_t-1}$ is done. Continue as above, the final $\overline{U}^{\star}_{\overline{N}, \overline{M}}$ will be more close to the optimal solution. In addition, if we consider to reduce the complexity of Algorithm 2, only one optimal control value $u^{\star}_{t+i_t|t-k_t}$ can be calculated each time, then all the predictive values constitute the final predictive control sequence.

The detailed steps of Algorithm 2 are given in the following.

 Algorithm 2 Approximate forward stagewise algorithm Step 1. Calculate $u^{\star}_{t|t-k_t}$ according to Remark 8 for the case of no delay in the forward channel. Step 2. For the forward communication delay $i_t=1, 2, \ldots, \overline{M}$: Step 2.1 Add $\hat{x}^T_{t+i_t+1|t-k_t}$ into state vector $\overline{X}$; Step 2.2 Add $u_{t+i_t|t-k_t}$ into input vector $\overline{U}$; Step 2.3 Extend matrixes $Q_{k_t, 0}$ and $R_{k_t, 0}$ to the corresponding dimensions; Step 2.4 Calculate $\overline{K}_{i_t}$ and $\overline{V}_{i_t}$ according to (17) and (18); Step 2.5 Obtain $\overline{U}^{\star}_{k_t, i_t}$ according to (20). Step 3. Take $\overline{U}^{\star}_{\overline{N}, \overline{M}}$ as the final predictive control sequence and send it to the actuator via the forward channel.

Analyzing the process of Algorithm 2, we can know that $\overline{M}$ times calculation are needed, and $\overline{M}(\overline{M}+1) /2$ control inputs will be derived. Compared to Algorithm 2, $\overline{M}$ control inputs will be calculated in each iteration in Algorithm 1, and many iterations are always needed before the variations of performance index less than the given threshold value, which leads to more calculation. Therefore, Algorithm 2 always has the faster calculation speed.

4 Simulation

In this section, an illustrative example is provided to verify the design scheme and algorithms developed in this paper.

Consider the following discrete-time bilinear system

 ${{x}_{t+1}}=A{{x}_{t}}+B{{u}_{t}}+{{u}_{t}}N{{x}_{t}}$ (21)

with the system matrixes

 \begin{align} &A=\left[ \begin{array}{*{35}{r}} 0.65&0.24 \\ -0.59&-0.27 \\ \end{array} \right] \\ &B=\left[ \begin{array}{*{35}{r}} -0.43 \\ 0.35 \\ \end{array} \right] \\ &N=\left[ \begin{array}{*{35}{r}} 1.49&0.34 \\ 1.08&-0.41 \\ \end{array} \right]. \\ \end{align}

The initial state of the system is chosen as $x_0=$ $[0.45 \quad-0.25]^T$, and the upper bound of communication delay and consecutive data dropouts in the forward channel and feedback channel are all 3, i.e., $\overline{M}=3$ and $\overline{N}=3$.

To illustrate the effectiveness of the networked predictive control scheme proposed in this paper, we consider the following two cases.

Case Ⅰ: Local control (offline) method presented in [25]. Here, the linear state-feedback controller is designed as $u_t$ = $Kx_t$, where the feedback gain $K=[-0.15 \quad0.34]$. Clearly, it can be obtained that $\| A+BK\|=0.8083 <1$, which means that the designed linear state-feedback law can guarantee the uniformly asymptotic stability of zero state.

Case Ⅱ: Networked predictive control scheme proposed in this paper. Applying the networked predictive control scheme proposed in Section Ⅱ and solving the predictive control sequence by using Algorithm 1 and Algorithm 2 presented in Section Ⅲ. The constant matrixes R and Q in performance index are unit matrixes with appropriate dimensions. The initial control inputs in $\overline{U}_0$ are set to be random numbers in $(-1, 1)$ and the threshold value $\lambda=0.5$. In addition, let $\overline{U}_0=0.5*\overline{U}_{\rm new}$ according to Remark 7, where $\overline{U}_{\rm new}$ represents the new initial control sequence in each sampling period.

The system (21) is controlled by the two kinds of controllers stated above, respectively. The state curves of $x_1$ and $x_2$, which are marked by Algorithm 1, Algorithm 2 and local control, are shown in Fig. 2 and Fig. 3, respectively.

 Figure 2 State x1 under Algorithm 1, Algorithm 2 and local control.
 Figure 3 State x2 under Algorithm 1, Algorithm 2 and local control.

The simulation results show that the predictive control sequences, calculated by Algorithm 1 and Algorithm 2, can guarantee the stability of the networked bilinear system while random time delays exist both in the forward channel and feedback channel. Meanwhile, the state curves under Algorithm 1 and Algorithm 2 have the similar trend with the local linear state-feedback controller [25], which shows the effectiveness of the proposed algorithms. Further according to Figs. 2 and 3, we can find that the system state under Algorithm 2 has the faster convergence speed and better control performance compared with Algorithm 1. In addition, the time spent on Algorithm 1 and Algorithm 2 are 0.159 s and 0.067 s, respectively, which shows that Algorithm 2 has lower computational complexity.

5 Conclusion

In this paper, the networked predictive control of bilinear systems is investigated. Two algorithms, called the stepwise iterative (SI) algorithm and approximate forward stagewise (AFS) algorithm, are proposed to solve a non-convex optimization problem to obtain the predictive control sequence. The SI algorithm begins with a given initial sequence of control inputs, and only one control input is calculated every time. Keeping on the iterative process until the variations of performance index meets the given threshold value $\lambda$. For faster calculation speed and better control performance, the AFS algorithm is proposed. A group of optimal values are calculated every time and the dimensions of control vector and state vector are added gradually until all the predictive control values are obtained. Then the control sequence will be transmitted to the actuator via the forward channel, and the appropriate control value will be selected according to the actual forward communication delay. Finally, a simulation demonstrates the effectiveness of the proposed predictive control scheme and optimization algorithms.

References
 1 R. A. Gupta and M. Y. Chow, "Networked control system:overview and research trends, " IEEE Trans. Ind. Electron., vol.57, no.7, pp.2527-2535, Nov.2010. https://www.researchgate.net/publication/224076713_Networked_Control_System_Overview_and_Research_Trends 2 W. P. M. H. Heemels, A. R. Teel, N. van de Wouw, and D. Nesic, "Networked control systems with communication constraints:tradeoffs between transmission intervals, delays and performance, " IEEE Trans. Autom. Control, vol.55, no. 8, pp.1781-1796, Feb.2010. https://www.researchgate.net/publication/224112751_Networked_Control_Systems_With_Communication_Constraints_Tradeoffs_Between_Transmission_Intervals_Delays_and_Performance 3 R. N. Yang, P. Shi, G. P. Liu, and H. J. Gao, "Network-based feedback control for systems with mixed delays based on quantization and dropout compensation, " Automatica, vol.47, no.12, pp.2805-2809, Dec.2011. http://www.sciencedirect.com/science/article/pii/S000510981100433X 4 M. F. Ren, J. H. Zhang, M. Jiang, and J. L. Xu, "Minimum (h, φ)-entropy control for non-Gaussian stochastic networked control systems and its application to a networked DC motor control system, " IEEE Trans. Control Syst. Technol., vol.23, no.1, pp.406-411, Jan.2015. https://www.researchgate.net/publication/280235329_Minimum_Entropy_Control_for_Non-Gaussian_Stochastic_Networked_Control_Systems_and_Its_Application_to_a_Networked_DC_Motor_Control_System 5 G. Nikolakopoulos, L. Dritsas, and S. S. Delshad, "Combined networked switching output feedback control with D-region stability for performance improvement, " Int. J. Control, vol.87, no.6, pp.1172-1180, Jun.2014. https://www.researchgate.net/publication/263495014_Combined_networked_switching_output_feedback_control_with_Inline_graphic-region_stability_for_performance_improvement 6 Luck R., Ray A.. Experimental verification of a delay compensation algorithm for integrated communication and control systems. Int. J. Control., 1994, 59(6): 1357-1372. DOI:10.1080/00207179408923135 7 H. J. Gao and T. W. Chen, "Network-based output tracking control, " IEEE Trans. Autom. Control, vol.53, no.3, pp.655-667, Apr.2008. https://www.researchgate.net/publication/3033073_Network-Based_Output_Tracking_Control 8 G. P. Liu, J. X. Mu, and D. Rees, "Networked predictive control of systems with random communication delay, " in Proc. UKACC Control, Bath, U. K. , 2004. 9 G. P. Liu, "Predictive controller design of networked systems with communication delays and data loss, " IEEE Trans. Circuits Syst. Ⅱ Exp. Briefs, vol.57, no.6, pp.481-485, Jun. 2010. https://www.researchgate.net/publication/224142479_Predictive_Controller_Design_of_Networked_Systems_With_Communication_Delays_and_Data_Loss 10 Y. Q. Xia, W. Xie, B. Liu, and X. Y. Wang, "Data-driven predictive control for networked control systems, " Inf. Sci., vol. 235, pp.45-54, Jun.2013. http://www.sciencedirect.com/science/article/pii/S0020025512000758 11 Z. H. Pang, G. P. Liu, and D. H. Zhou, "Design and performance analysis of incremental networked predictive control systems, " IEEE Trans. Cyber., vol.46, no.6, pp.1400-1410, Jun. 2016. https://www.researchgate.net/publication/280116504_Design_and_Performance_Analysis_of_Incremental_Networked_Predictive_Control_Systems 12 J. H. Zhang, J. Lam, and Y. Q. Xia, "Output feedback delay compensation control for networked control systems with random delays, " Inf. Sci., vol.265, pp.154-166, May2014. http://www.sciencedirect.com/science/article/pii/S0020025513008712 13 Liu G. P., Mu J. X., Rees D., Chai S. C.. Design and stability analysis of networked control systems with random communication time delay using the modified MPC. Int. J. Control, 2006, 79(4): 288-297. DOI:10.1080/00207170500533288 14 Mu J. X, Liu G. P., Rees D.. "Design of robust networked predictive control systems, " in Proc. 2005 American Control Conf. , Portland, OR, USA, 2005(1): 638-643. 15 Mohler R.R.. Bilinear Control Processes:with Applications to Engineering, Ecology and Medicine. Orlando: Academic Press, 1973. 16 Marchesini G., Mitter S.K.. Mathematical Systems Theory. Berlin Heidelberg: Springer, 1976. 17 Y. B. Zhao, G. P. Liu, and D. Rees, "A predictive control-based approach to networked Hammerstein systems:design and stability analysis, " IEEE Trans. Syst. Man Cybern. B, Cybern., vol.38, no.3, pp.700-708, Jun.2008. https://www.ncbi.nlm.nih.gov/pubmed/18558535/# 18 Y. B. Zhao, G. P. Liu, and D. Rees, "A predictive control based approach to networked wiener systems, " Int. J. Innov. Comput. Inf. Control., vol.4, no.11, pp.2793-2802, Nov. 2008. https://www.researchgate.net/publication/264840709_A_predictive_control_based_approach_to_networked_wiener_systems 19 H. Ouyang, G. P. Liu, D. Rees, and W. Hu, "Predictive control of networked non-linear control systems, " Proc. Inst. Mech. Eng. J. Syst. Control Eng., vol.221, no.3, pp.453-463, May2007. https://www.researchgate.net/publication/245389365_Predictive_control_of_networked_non-linear_control_systems 20 G. Pin and T. Parisini, "Networked predictive control of uncertain constrained nonlinear systems:recursive feasibility and input-to-state stability analysis, " IEEE Trans. Autom. Control, vol.56, no.1, pp.72-87, Jan.2011. https://www.researchgate.net/publication/220387079_Networked_Predictive_Control_of_Uncertain_Constrained_Nonlinear_Systems_Recursive_Feasibility_and_Input-to-State_Stability_Analysis 21 G. P. Liu, "Design and analysis of networked non-linear predictive control systems, " IET Control Theory Appl., vol.9, no.11, pp.1740-1745, Jul.2015. 22 P. E. Gill and E. Wong, "Sequential quadratic programming methods", in Mixed Integer Nonlinear Programming, J. Lee and S. Leyffer, Eds. , New York, USA, Springer, 2012. 23 M. Hotz and C. Vogel, "Linearization of time-varying nonlinear systems using a modified linear iterative method, " IEEE Trans. Signal Process., vol.62, no.10, pp.2566-2579, May 2014. http://www.oalib.com/paper/4081839 24 A. Bidram, F. L. Lewis, and A. Davoudi, "Synchronization of nonlinear heterogeneous cooperative systems using input-output feedback linearization, " Automatica., vol.50, no.10, pp.2578-2585, Oct.2014. http://dl.acm.org/citation.cfm?id=2941006.2941225 25 Z. L. Mu, F. Liu, and Z. L. Qiu, "Analysis of stabilizing control of discrete bilinear system, " Control Theory Appl., vol.17, no.2, pp.267-269, Apr.2000. http://en.cnki.com.cn/Article_en/CJFDTOTAL-KZLY200002025.htm