以网络作为通信传输介质的闭环反馈控制系统称为网络控制系统，又称网络化控制系统(networked control systems, NCS).网络控制系统的优点主要在于节约成本、便于实现资源共享、系统可靠性和模块性较高等^[1-3].正因如此，对于网络控制系统的研究在近年逐渐成为学术界的一大热点.

将网络引入控制系统在获得众多好处的同时，也使得对于其结构及数学模型的分析变得相对困难.这是因为，含有网络模块的控制系统不得不面对网络本身就存在的众多问题.而其中最主要的两个问题是网络诱导时延和数据包的丢失，它们极有可能成为控制系统性能下降的诱因，甚至造成整个控制系统的不稳定.因此，时延和丢包的问题在网络控制系统的数学模型建立中是不得不考虑的两大因素^[4].

网络控制系统的一般结构如图 1所示.系统的3个主要节点是传感器、控制器和执行器.而根据传感器节点、控制器节点和执行器节点工作方式的不同，所建立的模型可以分为离散模型、连续模型和混合模型^[5].文献[6]在离散模型的基础上，对时延的网络控制系统进行了分析，采用了针对二次型性能指标的最优控制策略.文献[7]对网络控制系统的有界丢包进行了研究，并使用了两种控制方法，其一保证了存在随机丢包时网络控制系统的大范围渐进稳定，另一种则保证了存在丢包的马尔可夫模型下的均方稳定性.然而，同时考虑时延和丢包因素的研究却不多，文献[8]在数据包无时许错乱的前提下，考虑时延和丢包的问题建立了网络控制系统的数学模型并进行了控制器设计.

本文结合了影响网络控制系统的两大因素，建立了时延和丢包的网络控制系统的离散数学模型，采用随机最优控制算法计算最优的二次型性能指标和相应的控制增益.

1 系统模型分析

假设被控对象为线性时不变系统，则系统模型可表示为

$ \mathit{\boldsymbol{\dot x}}\left( t \right) = \mathit{\boldsymbol{Ax}}\left( t \right) + \mathit{\boldsymbol{Bu}}\left( t \right), $

(1)

其中x (t)∈R_n, u (t)∈R_m分别表示系统状态向量和控制信号向量，而A ∈R_n×n，B∈R_n×m表示系统矩阵.

1.1 假设

本文建立的带有时延和丢包的网络控制系统的数学模型，需要作出如下假设：

1) 传感器为时间驱动，控制器和执行器为事件驱动；

2) 传感器与控制器之间的时延和控制器与执行器的之间的时延相互独立，且可以将两者之和作为系统总时延进行分析，网络诱导总时延有界；

3) 如果最近的控制信号被执行，则之前的控制信号被丢弃^[9].

1.2 数学模型建立

考虑网络控制系统的丢包因素，文献[10]提出了引入二元随机变量φ来表示在采样周期内控制命令是否传递给了被控对象，本文同样采用这种方法.如果控制信号被接受到，则其为单位矩阵；反之控制信号数据丢失则为零矩阵，φ∈R_m×m.同时考虑两种网络诱导时延之和τ (t)= τ_sc(t)+ τ_ca(t)，则式(1)在引入丢包和时延因子之后被重写为

$ \mathit{\boldsymbol{\dot x}}\left( t \right) = \mathit{\boldsymbol{Ax}}\left( t \right) + \mathit{\boldsymbol{B\varphi }}\left( t \right)\mathit{\boldsymbol{u}}\left( {t - \tau \left( t \right)} \right). $

(2)

由于控制器的驱动方式为事件驱动，设采样周期时长为T.根据文献[6]引入的时延限制变量h(h为大于等于1的整数)，使得两种网络诱导时延之和有界：

$ 0 \le \tau < hT. $

考虑网络控制系统的驱动方式，将式(2)离散化之后，带有时延和丢包的网络控制系统在采样周期[KT, (K+1)T]内数学模型表示为

$ \begin{array}{l} {\mathit{\boldsymbol{x}}_{k + 1}} = {\mathit{\boldsymbol{A}}_s}{\mathit{\boldsymbol{x}}_k} + \mathit{\boldsymbol{B}}_0^k{\mathit{\boldsymbol{\varphi }}_0}{\mathit{\boldsymbol{u}}_k} + \mathit{\boldsymbol{B}}_1^k{\mathit{\boldsymbol{\varphi }}_1}{\mathit{\boldsymbol{u}}_{k - 1}} + \cdots + \\ \mathit{\boldsymbol{B}}_h^k{\mathit{\boldsymbol{\varphi }}_k}{\mathit{\boldsymbol{u}}_{k - h}}, \end{array} $

(3)

其中, x_k+1= x(kT)，u_k=u(kT), A_s=e^AT.

系统控制信号时序图如图 2所示，则系统在采样周期内可能接收到控制信号的时刻为

$ kT + {t_i},{t_i} = {\tau _i} - iT,i = 0,1, \cdots ,h. $

同时，由于最近的控制信号被执行时之前的控制信号被丢弃，引入变量θ(x)判断在采样周期[KT, (K+1)T]内控制信号是否可能被执行，那么

$ \mathit{\boldsymbol{B}}_i^k = \int_{\tau _i^k - iT}^{\tau _{i - 1}^k - \left( {i - 1} \right)T} {{{\rm{e}}^{A\left( {T - s} \right)}}{\rm{d}}s\mathit{\boldsymbol{B}} \cdot \mathit{\boldsymbol{\theta }}\left( {T + \tau _{i - 1}^k - \tau _i^k} \right) \cdot \mathit{\boldsymbol{\theta }}\left( {\tau _i^k - iT} \right)} , $

$ \forall i = 0,1, \cdots ,h. $

而对于本周期发出的控制信号，如果其时延大于一个周期，则无法在本周期内执行.所以，其系数应表示为

$ \mathit{\boldsymbol{B}}_0^i = \int_{\tau _0^k}^T {{{\rm{e}}^{A\left( {T - s} \right)}}{\rm{d}}s\mathit{\boldsymbol{B}} \cdot \mathit{\boldsymbol{\theta }}\left( {T - \tau _0^k} \right)} , $

$ \mathit{\boldsymbol{\theta }}\left( x \right) = \left\{ \begin{array}{l} 0,x \le 0,\\ {\mathit{\boldsymbol{I}}_m},x > 0. \end{array} \right. $

$ {\mathit{\boldsymbol{\varphi }}_i} = \left\{ \begin{array}{l} 0,{\mathit{\boldsymbol{u}}_{k - i}}在\left[ {kT,\left( {k + 1} \right)T} \right]内接收,\\ {\mathit{\boldsymbol{I}}_m},{\mathit{\boldsymbol{u}}_{k - i}}在\left[ {kT,\left( {k + 1} \right)T} \right]内丢失. \end{array} \right. $

引入新的状态向量为原状态向量与周期内控制信号的增广矩阵得到

$ {\mathit{\boldsymbol{z}}_k} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{x}}_k}}\\ {{\mathit{\boldsymbol{u}}_{k - 1}}}\\ \vdots \\ {{\mathit{\boldsymbol{u}}_{k - h}}} \end{array}} \right], $

参照式(3)的系统矩阵，由新状态向量z_k表示的系统模型为

$ {\mathit{\boldsymbol{z}}_{k + 1}} = {\mathit{\boldsymbol{A}}_k}{\mathit{\boldsymbol{z}}_k} + {\mathit{\boldsymbol{B}}_k}{\mathit{\boldsymbol{u}}_k}, $

(4)

其中

$ {\mathit{\boldsymbol{A}}_k} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{A}}_s}}&{\mathit{\boldsymbol{B}}_1^k{\mathit{\boldsymbol{\varphi }}_1}}& \cdots &{\mathit{\boldsymbol{B}}_i^k{\mathit{\boldsymbol{\varphi }}_i}}& \cdots &{\mathit{\boldsymbol{B}}_h^k{\mathit{\boldsymbol{\varphi }}_h}}\\ 0&0& \cdots&\cdots&\cdots &0\\ 0&{{\mathit{\boldsymbol{I}}_m}}& \cdots&\cdots&\cdots &0\\ \vdots &0&{{\mathit{\boldsymbol{I}}_m}}& \cdots&\cdots &0\\ \vdots&\vdots &0&{}&0& \vdots \\ 0&0& \cdots&\cdots &{{\mathit{\boldsymbol{I}}_m}}&0 \end{array}} \right], $

$ {\mathit{\boldsymbol{B}}_k} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{B}}_0^k{\mathit{\boldsymbol{\varphi }}_0}}\\ {{\mathit{\boldsymbol{I}}_m}}\\ 0\\ 0\\ \vdots \\ 0 \end{array}} \right]. $

2 随机最优控制 2.1 假设

本文采用的随机最优控制算法的前提条件包括：系统状态可以完全观测；时延的丢包之间相互独立，且两者分布函数已知.

2.2 随机最优控制算法

针对带有时延和丢包的网络控制系统的离散数学模型，选取与文献[9]相同的二次型代价函数作为最优控制的性能指标.

采样周期[KT, (K+1)T]内的性能指标可以表示为

$ {\mathit{\boldsymbol{J}}_k} = E\left[ {\sum\limits_{i = k}^\infty {\left( {\mathit{\boldsymbol{x}}_i^{\rm{T}}Q{\mathit{\boldsymbol{x}}_i} + \mathit{\boldsymbol{u}}_i^{\rm{T}}R{\mathit{\boldsymbol{u}}_i}} \right)} } \right],k = 0,1, \cdots , $

其中Q为半正定矩阵，而R为正定矩阵.

引入增广矩阵z_k，重写代价函数后性能指标表示为

$ {\mathit{\boldsymbol{J}}_k} = E\left[ {\sum\limits_{i = k}^\infty {\left( {\mathit{\boldsymbol{z}}_i^{\rm{T}}{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{z}}}{\mathit{\boldsymbol{z}}_i} + \mathit{\boldsymbol{u}}_i^{\rm{T}}{\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{z}}}{\mathit{\boldsymbol{u}}_i}} \right)} } \right],k = 0,1, \cdots , $

其中

$ {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{z}}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{Q}}&0& \cdots &0\\ 0&{\frac{\mathit{\boldsymbol{R}}}{h}}& \cdots &0\\ \vdots&\vdots &{}& \vdots \\ 0&0& \cdots &{\frac{\mathit{\boldsymbol{R}}}{h}} \end{array}} \right],而\;{\mathit{\boldsymbol{R}}_z} = \frac{\mathit{\boldsymbol{R}}}{h}. $

其中Q_z仍为半正定矩阵，而R_z仍为正定矩阵.

若不考虑网络控制系统的噪声，则其最优控制律可以写成

$ {\mathit{\boldsymbol{u}}_k} = - \mathit{\boldsymbol{L}}{\mathit{\boldsymbol{z}}_k}, $

(5)

其中L称为最优增益，而由最优控制理论，随机代价函数可以表示为

$ {\mathit{\boldsymbol{J}}_k} = E\left[ {\mathit{\boldsymbol{z}}_k^{\rm{T}}{\mathit{\boldsymbol{P}}_k}{\mathit{\boldsymbol{z}}_k}} \right],{\mathit{\boldsymbol{P}}_k} \ge 0, $

其中P_k成为随机黎卡提方程的解.

那么，由离散系统动态规划方法，原性能指标可以被重写为

$ \begin{array}{l} {\mathit{\boldsymbol{J}}_k} = E\left[ {\mathit{\boldsymbol{z}}_k^{\rm{T}}{\mathit{\boldsymbol{P}}_k}{\mathit{\boldsymbol{z}}_k} + \mathit{\boldsymbol{u}}_k^{\rm{T}}{\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{z}}}{\mathit{\boldsymbol{u}}_k} + \sum\limits_{i = k + 1}^\infty {\left( {\mathit{\boldsymbol{z}}_i^{\rm{T}}{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{z}}}{\mathit{\boldsymbol{z}}_i} + } \right.} } \right.\\ \left. {\left. {\mathit{\boldsymbol{u}}_i^{\rm{T}}{\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{z}}}{\mathit{\boldsymbol{u}}_i}} \right)} \right] = E\left[ {\mathit{\boldsymbol{z}}_k^{\rm{T}}{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{z}}}{\mathit{\boldsymbol{z}}_k} + \mathit{\boldsymbol{u}}_k^{\rm{T}}{\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{z}}}{\mathit{\boldsymbol{u}}_k} + {\mathit{\boldsymbol{J}}_{k + 1}}} \right] = \mathit{\boldsymbol{z}}_k^{\rm{T}}{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{z}}}{\mathit{\boldsymbol{z}}_k} + \\ \mathit{\boldsymbol{u}}_k^{\rm{T}}{\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{z}}}{\mathit{\boldsymbol{u}}_k} + E\left[ {\mathit{\boldsymbol{z}}_{k + 1}^{\rm{T}}{\mathit{\boldsymbol{P}}_{k + 1}}{\mathit{\boldsymbol{z}}_{k + 1}}} \right], \end{array} $

将式(4)代入后得到

$ \begin{array}{l} {\mathit{\boldsymbol{J}}_k} = {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{z}}_k}}\\ {{\mathit{\boldsymbol{u}}_k}} \end{array}} \right]^{\rm{T}}}\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_z}}&0\\ 0&{{\mathit{\boldsymbol{R}}_z}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{z}}_k}}\\ {{\mathit{\boldsymbol{u}}_k}} \end{array}} \right] + \\ E\left( {{{\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{z}}_k}}\\ {{\mathit{\boldsymbol{u}}_k}} \end{array}} \right]}^{\rm{T}}}\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{A}}_z^{\rm{T}}}\\ {\mathit{\boldsymbol{B}}_z^{\rm{T}}} \end{array}} \right]{\mathit{\boldsymbol{P}}_{k + 1}}{{\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{A}}_z^{\rm{T}}}\\ {\mathit{\boldsymbol{B}}_z^{\rm{T}}} \end{array}} \right]}^{\rm{T}}}\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{z}}_k}}\\ {{\mathit{\boldsymbol{u}}_k}} \end{array}} \right]} \right) = {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{z}}_k}}\\ {{\mathit{\boldsymbol{u}}_k}} \end{array}} \right]^{\rm{T}}} \times \\ \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_z} + E\left( {\mathit{\boldsymbol{A}}_z^{\rm{T}}{\mathit{\boldsymbol{P}}_{k + 1}}{\mathit{\boldsymbol{A}}_z}} \right)}&{E\left( {\mathit{\boldsymbol{A}}_z^{\rm{T}}{\mathit{\boldsymbol{P}}_{k + 1}}{\mathit{\boldsymbol{B}}_z}} \right)}\\ {E\left( {\mathit{\boldsymbol{B}}_z^{\rm{T}}{\mathit{\boldsymbol{P}}_{k + 1}}{\mathit{\boldsymbol{A}}_z}} \right)}&{{\mathit{\boldsymbol{R}}_z} + E\left( {\mathit{\boldsymbol{B}}_z^{\rm{T}}{\mathit{\boldsymbol{P}}_{k + 1}}{\mathit{\boldsymbol{B}}_z}} \right)} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{z}}_k}}\\ {{\mathit{\boldsymbol{u}}_k}} \end{array}} \right], \end{array} $

根据动态规划方法，将J_k对u_k求偏微分，易求得最优增益为

$ \mathit{\boldsymbol{L}} = {\left[ {{\mathit{\boldsymbol{R}}_z} + E\left( {\mathit{\boldsymbol{B}}_z^{\rm{T}}{\mathit{\boldsymbol{P}}_{k + 1}}{\mathit{\boldsymbol{B}}_z}} \right)} \right]^{ - 1}}E\left( {\mathit{\boldsymbol{B}}_z^{\rm{T}}{\mathit{\boldsymbol{P}}_{k + 1}}{\mathit{\boldsymbol{A}}_z}} \right), $

从而可以由式(5)计算出最优控制律.

3 仿真研究 3.1 数值模型建立

为了便于对比加入丢包算子之后的响应曲线，本文选择与文献[6]相同的数值模型，即线性化后的倒立摆简化的状态方程^[11]：

$ \mathit{\boldsymbol{\dot x}} = \mathit{\boldsymbol{Ax}} + \mathit{\boldsymbol{Bu}} = \left[ {\begin{array}{*{20}{c}} 1&0\\ 0&1 \end{array}} \right]\mathit{\boldsymbol{x + }}\left[ {\begin{array}{*{20}{c}} 0\\ 1 \end{array}} \right]\mathit{\boldsymbol{u}}. $

以此不稳定的被控对象建立网络控制系统，设定时延限制h=2，即0≤ τ < 2T，采样周期T=0.1 s，将网络控制系统离散化后得到模型：

$ {\mathit{\boldsymbol{x}}_{k + 1}} = {\mathit{\boldsymbol{A}}_s}{\mathit{\boldsymbol{x}}_k} + \mathit{\boldsymbol{B}}_0^k{\mathit{\boldsymbol{\varphi }}_0}{\mathit{\boldsymbol{u}}_k} + \mathit{\boldsymbol{B}}_1^k{\mathit{\boldsymbol{\varphi }}_1}{\mathit{\boldsymbol{u}}_{k - 1}} + \cdots + \mathit{\boldsymbol{B}}_2^k{\mathit{\boldsymbol{\varphi }}_2}{\mathit{\boldsymbol{u}}_{k - 2}}, $

其中，

$ \mathit{\boldsymbol{B}}_0^i = \int_{\tau _0^k}^T {{{\rm{e}}^{A\left( {T - s} \right)}}{\rm{d}}s\mathit{\boldsymbol{B}} \cdot \mathit{\boldsymbol{\theta }}\left( {T - \tau _0^k} \right)} , $

$ \mathit{\boldsymbol{B}}_1^k = \int_{\tau _1^k - T}^{\tau _0^k} {{{\rm{e}}^{A\left( {T - s} \right)}}{\rm{d}}s\mathit{\boldsymbol{B}} \cdot \mathit{\boldsymbol{\theta }}\left( {T + \tau _0^k - \tau _1^k} \right) \cdot \mathit{\boldsymbol{\theta }}\left( {\tau _1^k - T} \right)} , $

$ \mathit{\boldsymbol{B}}_2^k = \int_{\tau _2^k - 2T}^{\tau _1^k - T} {{{\rm{e}}^{A\left( {T - s} \right)}}{\rm{d}}s\mathit{\boldsymbol{B}} \cdot \mathit{\boldsymbol{\theta }}\left( {T + \tau _1^k - \tau _2^k} \right) \cdot \mathit{\boldsymbol{\theta }}\left( {\tau _2^k - 2T} \right)} . $

选择仿真参数为$\mathit{\boldsymbol{Q}} = \left[ \begin{array}{l} 1\;\;0\\ 0\;\;1 \end{array} \right], $R=0.5，系统初值x₀=[－3 5]^T，时延为[0, 2T)内λ=0.8T的泊松分布序列，丢包服从伯努利分布.

3.2 仿真结果分析

在本文的数值模型基础上，经过Matlab仿真得到的图 3(a)和图 3(b)分别给出了一组丢包率约为0.1和约为0.2时的系统状态响应曲线.仿真结果表明，丢包率约为0.1的系统模型在随机最优控制算法的控制下效果明显好于丢包率约为0.2的系统模型.事实上，当该模型中丢包率小于等于0.2时，系统状态均能逐渐达到稳定值并最终收敛.而多次仿真实验都显示出一旦丢包率大于等于0.3，系统将无法达到稳定，状态发散.这是因为当系统丢包率超过临界值后，随机黎卡提方程无解，文献[12]详细地论证了随机黎卡提方程有解的条件，并给出了丢包率的临界值的分析.通过仿真结果和以上分析不难得出结论：该倒立摆数学模型丢包率的临界值应该在[0.2, 0.3]的区间内.

4 结语

网络控制系统的时延和丢包问题在本文所描述的数学模型中都得到了体现.在系统状态可完全观测的前提条件下，本文针对二次型性能指标，采用随机最优控制算法得到了系统的最优增益的表达形式.仿真结果表明了随机最优控制对于时延和丢包的网络控制系统的有效性和倒立摆模型的丢包率临界值的范围.考虑系统白噪声时的数学模型和控制算法还有待于进一步研究.