-
摘要: 高精度时间同步是任务关键型工业网络控制系统的核心支撑技术, 针对工业环境中普遍存在周期性振动等扰动信号导致晶振频率漂移, 影响时间同步精度的问题, 基于扩展比例积分(Proportional integral, PI)观测器, 提出一种新型的抗扰补偿器结构, 用于消除周期性扰动的影响, 并构建了相应的精细抗干扰反馈控制方法, 用于实现高精度时间同步. 与传统的扰动观测器相比, 所提出的扩展$ \mathrm{P}\mathrm{I} $抗扰补偿器克服了传统扰动观测器零点不变局限性, 提出了零点配置方法, 以充分利用闭环系统的传递函数矩阵(Transfer function matrix, TFM)在系统零点处降秩的特性, 实现了对于特定频率扰动信号的补偿作用. 并给出了所提出的控制器和抗扰补偿器的稳定性证明和控制器参数的稳定域. 通过基于实测参数的无线网络仿真实验, 验证了在$5\;\mathrm{g}$周期性振动干扰下, 提出的方法明显优于传统滤波器和补偿器, 达到了同步误差在4 μs以内, 实现了高精度时间同步.Abstract: Precision time synchronization is one of the enabling technologies for mission-critical industrial networked control systems. To address the precision degradation due to crystal oscillator's frequency drift caused by external disturbances and environment changes, such as the periodic acceleration and vibration which are and commonly found in industrial environments, this paper proposes a novel disturbance compensation structure based on an extended proportional integral (PI) observer for effectively eliminating the effects of periodic disturbances. A disturbance rejection feedback control method is proposed to achieve the precision time synchronization. The proposed extended PI disturbance compensator overcomes the zero-point invariance limitation of the conventional disturbance observer and a zero-point configuration is adopted to make full use of the rank deficiency of the transfer function matrix (TFM) at zeros of the closed-loop system to achieve the compensation effect for frequency-specific disturbance signals. The stability of the proposed controller and disturbance compensator is proved and the convergence domain of the controller parameters is also given. Through the simulation experiments of the wireless network based on measured parameters, it is verified that the proposed method significantly outperforms the conventional observers and compensators and synchronization errors within 4 μs are achieved under $ 5\;\mathrm{g} $ periodic vibration.
-
工业系统通常是任务关键型时间敏感系统, 随着工业物联网的发展, 如智能电网、车联网、机器人/无人机群协同、分布式工业控制系统等都需要时间敏感网络的支撑. 时间同步作为时间敏感网络关键技术之一, 对目标跟踪定位、协同控制、数据融合、介质访问控制等起着关键作用. 通常的嵌入式物联网节点是以石英晶体振荡器作为时钟源, 典型石英晶体振荡器的精度通常在几十$ \mathrm{p}\mathrm{p}\mathrm{m} $(Part per million)至$ 100\;\mathrm{p}\mathrm{p}\mathrm{m} $左右[1], 且由于制造误差、温度、压力、加速度等环境变化而普遍存在频率漂移, 造成时钟精度降低[2]. 因此通常采用时间同步算法, 对漂移时钟进行调控, 使得各个漂移时钟的时间值达到较高精度的一致性, 以满足应用需求.
现有时间同步协议主要采用带时间戳的包交换方式来获取节点间的时间差值, 并基于该时间差对漂移时钟的相位(即时间值)或者频率进行调整, 以减少时间差值, 达到同步. 常见的时间同步算法有, 基于单向包交换的“接收者−接收者”同步, 如参考广播同步协议 (Reference broadcast synchronization, RBS)[3]等, 虽然可避免双向包交换中普遍存在的传输延迟不对称性问题, 但是面临着只能在同级节点之间同步, 不适用于大规模网络. 为克服扩展性差的不足, 目前大多数时间同步协议采用的是双向包交换, 如IEEE1588协议[4]等.
从控制论的角度来看, 时间同步的本质是一个动态系统的状态估计和反馈控制问题, 即利用时间戳包交换获取到包含了测量噪声的时间差测量值, 估计出时间差的准确值, 然后对漂移时钟的相位或者频率进行调整. 在状态估计方面, 采用较多的是线性拟合对时间差进行估计如泛洪时间同步 (Flooding time synchronization protocol, FTSP)[5]、多参考节点时间同步方法[6], 以及Wang等[7]提出的基于卡尔曼滤波的时钟频率漂移的估计方法, Yang等[8]提出的贝叶斯估计的时间同步算法等. 在时钟调控方面, 传统的时间同步算法如FTSP, RBS等是一种控制增益为1的比例控制, 直接将测量或者估计得到的时间误差调整本地时钟. 新近的研究如Carli等[9]提出的通过比例积分控制器调节时钟相位和时钟频率, 其中主要通过积分作用抑制时钟噪声和延时抖动的影响. 然而其面对复杂网络收敛速度较为缓慢, 为了更快收敛, Yildirim[10]提出了一种自适应时钟同步方法, 根据时钟偏差的大小动态调整控制器参数, 能更快达到同步, 但其还是局限于比例积分控制器的调整策略, 泛化性有待提高. 而后提出的PISync时间同步算法[11], 在自适应比例积分控制器的基础上, 结合泛洪的分布式协议改善了系统的稳态误差和可扩展性, 使得时间同步算法适用于更加复杂的网络化拓扑结构.
目前时间同步的研究主要是考虑时钟相位噪声以及通信延迟、计算机处理延迟等造成的随机误差, 取得了较好的效果. 但对于工业环境中由于往复运动、振动等造成晶振信号中出现特定频率周期性扰动[12], 导致时间同步精度变差的研究还没有. 如旋转机械、周期性振动、飞行器在飞行过程中机身出现低频振动[13], 都是一些有限带宽特定频率的扰动信号, 限制了时间同步精度. 因此需要对于特定频率的扰动进行有针对性的处理, 实现高精度的时间同步.
控制理论的研究中, 对于特定扰动的处理目前有内模抗扰控制(逆模型)[14]、自适应评判方法[15]、基于扰动观测器(Disturbance observer, DOB)[16−17]的抗扰控制等, 通常需要预先知道扰动的动力学模型, 构建扰动的状态空间模型, 并将扰动作为扩展的状态, 构建增广的龙伯格观测器, 实现对扰动的估计, 进而予以相应的反馈补偿, 消除扰动的影响. 针对模型未知的扰动, 主要有自抗扰控制 (Active disturbance rejection control, ADRC)[18−19], 采用扩张状态观测器(Extended state observer, ESO)实现对未知扰动的估计. Gao[20]采用极点配置的思想将ADRC参数与频率联系起来, 实现了更为简化的参数整定方法. 但是值得指出的是这些方法主要通过极点优化来实现抗扰性能的提高. Dai等[21]针对故障检测中的特定频率扰动, 提出了一种扩展PI观测器和相应的零极点联合优化方法, 利用传递函数矩阵(Transfer function matrix, TFM)的传输零点, 使得传递函数矩阵降秩, 对于特定频率的扰动取得了更好的抑制效果.
本文针对时间同步中特定频率扰动的抗扰控制问题, 基于零点配置的思想, 提出了一种新的扩展PI抗扰补偿器结构, 采用零极点联合优化方法, 且不依赖于扰动的精确模型, 实现对本地漂移时钟的反馈补偿控制, 提高了周期性扰动下时间同步系统的精度和抗扰性, 满足时间敏感工业应用的需要. 本文的主要创新点如下:
1) 提出了一种新的基于扩展PI抗扰补偿器的高精度时间同步系统框架, 将扩展PI观测器和二自由度内膜控制的原理相结合. 区别于内膜控制[14]、DOB[16], 本文提出的方法不依赖于扰动的具体模型, 且较Dai等[21]提出的扩展PI观测器而言, 原理性区别在于本文的动态反馈回路是作用于实际的被控对象, 而不是作用于被控对象的模型, 这也是和传统观测器的根本区别.
2)控制器和补偿器的设计采用了二自由度控制的思想, 更好地兼顾了抗扰性和控制精度. 内环为抗扰补偿作用, 主要通过所设计的抗扰补偿器, 有效地减少了外部环境变化引发周期性扰动的影响, 使得带补偿器的有扰时钟能逼近于理想时钟, 从而提高了时间同步系统的抗扰性和鲁棒性. 外环为相位跟踪作用, 主要通过所设计的控制器消除时间同步系统稳态误差, 提高同步精度. 在内外环共同作用下, 能同时兼顾并维持较好的动态性能和稳态误差.
3)本文扩展PI抗扰补偿器的反馈回路设计采用零极点优化的思想, 使得含有该补偿回路的带扰动被控系统逐步逼近于理想无扰动系统, 从而实现了对特定频率扰动的补偿和抑制. 并采用零点配置方法, 利用多变量系统的零点, 更好地抑制了特定频率的扰动, 并更好地实现了周期性扰动下的抗扰控制.
本文结构如下: 第1节构建了工业网络控制系统中的节点时钟模型; 第2节提出了带抗扰补偿器的时间同步系统框架; 第3节和第4节分别介绍了扩展PI抗扰补偿器和控制器的稳定性证明及参数优化方法; 最后通过与传统滤波器和扰动观测器对比, 验证了扩展PI抗扰补偿器的良好的抗扰补偿作用.
1. 问题描述
现有大多数通信网络是基于包交换的数字通信网络. 基于包交换的时钟同步问题, 是多个节点之间通过带时间信息的包交换, 获取各个时钟之间的时间偏差, 进而对本地时钟的相位和频率进行相应调整, 以实现所有节点的时钟具有相同的时间值.
1.1 时钟模型
工业网络系统中, 本地时间信息通常是由本地的以某个特定频率运行的晶体振荡器提供. 精准时钟是指晶体振荡器频率偏差始终为零的理想时钟; 晶振频率受晶体自身因素以及外界环境影响而存在偏差和相位噪声的时钟, 则为实际有偏时钟.
设时钟的时钟值和时钟速率分别用变量$ c $和$ \gamma $来表示, 采样间隔为$ h $, 则在第$ k $个采样时刻$ {t}_{k}=kh $, 时钟的时间值可以表达为
$$ c\left(k+1\right)=c\left(k\right)+\gamma \left(k\right)h+{\omega }_{c}\left(k\right) $$ (1) 其中, $ {\omega }_{c}\left(k\right) $是时钟的相位噪声, $ \gamma \left(k\right)h $为时钟经过采样间隔$ h $后的时间值的增量. 值得指出的是, 对于理想时钟的速率$ \gamma \left(k\right)\equiv 1 $, 而实际工业网络嵌入式系统中广泛采用的晶体振荡器时钟, 由于制造容差、环境温度、机械振动以及器件老化等原因, 其频率会发生漂移, 因此对于实际时钟有$ \gamma \left(k\right)\ne 1 $, 时钟速率在1左右漂移. 定义时钟频率偏斜$ \chi \left(k\right)= \gamma \left(k\right)-1 $, 用于表征实际时钟速率相对于理想标称时钟速率的偏移, 单位为$ \mathrm{p}\mathrm{p}\mathrm{m} $.
值得指出的是, 时钟速率$ \gamma $变化存在一定的随机性, 如经典的Brownian随机运动模型[22], 同时, 工作环境, 如温度、供电、机械振动等变化时, 晶振频率也会随之发生波动. 如晶振受到外力(外界加速度)作用时, 会发生形变, 造成频率偏斜, 相当于产生了频率调制, 频率的变化和加速度成比例关系[23]. 周期性的加速度会导致频率周期性波动, 所形成的频率扰动可以用一个具有特定频率特征的周期性扰动信号来描述. 因此, 时钟速率$ \gamma \left(k\right) $可建模为
$$ \gamma \left(k\right)={\gamma }_{0}+f\left(k\right)+\xi \left(k\right)\tag{2a} $$ 其中, $ {\gamma }_{0} $为一个恒定的或者缓慢变化频率偏差, 如常见的无线传感网节点平台Atmel SMART SAM R21, $ {\gamma }_{0}-1 $通常在$ \pm 20\;\mathrm{p}\mathrm{p}\mathrm{m} $~ $ \pm 50\;\mathrm{p}\mathrm{p}\mathrm{m} $[1], $ f\left(k\right) $用于表征加速度等环境变化导致的周期性扰动, 可以用一组正弦信号来描述$ f\left(k\right)=\sum {a}_{i}\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{i}k\right) $, 振幅$ {a}_{i} $和频率$ {\omega }_{i} $由晶振本身特性和外界加速度扰动等决定[23]. $ \xi \left(k\right) $描述热噪声的随机扰动, 通常可用一个Brownian随机运动过程描述[23].
不失一般性, 时钟速率的动态特性可以用差分方程描述, 写作$ \gamma \left(k+1\right)=\gamma \left(k\right)+{\omega }_{\gamma }\left(k\right) $, 其中$ {\omega }_{\gamma }\left(k\right) $为时钟速率噪声, 由扰动$ f\left(k\right) $和噪声$ \xi \left(k\right) $决定, 可定义为
$$ {\omega }_{\gamma }\left(k\right)=\gamma \left(k+1\right)-\gamma \left(k\right)={f'} \left(k\right)+{\xi '} \left(k\right)\tag{2b} $$ 其中, ${f'} \left(k\right)=f\left(k+1\right)-f\left(k\right)$为周期性信号$ f\left(k\right) $的差分形式, 对应连续系统的微分形式, 因此, ${f'} \left(k\right)$同样可以由一个周期性信号来描述. 同理, ${\xi '} \left(k\right)$可以由一个高斯过程来描述.
引入时钟状态变量$x\left(k\right)={\left[c\left(k\right)\;\gamma \left(k\right)\right]}^{\mathrm{T}}$, 定义时钟扰动信号$d\left( k \right) = {\left[ {{\omega _c}\left( k \right)}\;\;{{\omega _\gamma }\left( k \right)} \right]^{\rm{T}}}$, 则一个无调控的自由时钟, 其状态空间模型可以写为
$$ x\left(k+1\right)={\boldsymbol{A}}x\left(k\right)+{\boldsymbol{B}}u(k)+{\boldsymbol{F}}d\left(k\right) $$ (3) 其中, 状态转移矩阵${\boldsymbol{A}}=\left[ {\begin{aligned} 1\;\;h\\ 0\;\;1 \end{aligned}} \right]$, 状态扰动输入矩阵$ {\boldsymbol{F}}={{\boldsymbol{I}}}_{2\times 2} $为$ 2\times 2 $单位阵, 时钟扰动信号$ d\left(k\right) $包含相位噪声$ {\omega }_{c}\left(k\right) $和时钟速率噪声$ {\omega }_{\gamma }\left(k\right) $. 相位噪声$ {\omega }_{c}\left(k\right) $是一个随机噪声, 可以用一个零均值的高斯过程来描述[1]; 时钟速率噪声$ {\omega }_{\gamma }\left(k\right) $由式(2b)描述.
1.2 基于时间戳的观测方程
在实际工业物联网硬件系统中, 本地时间是通过读取本地时钟的寄存器值, 获取时钟的状态值$ c\left(k\right) $, 称作时间戳. 在读取寄存器获取时间戳的过程中, 不可避免地存在某些操作延迟, 如总线占用、函数调用、操作系统任务调度等, 造成所获得的时间戳与实际值有偏差. 因此, 第$ k $次本地时钟值进行观测时, 记有误差的时间戳为$ y\left(k\right) $, 时钟观测过程可以描述为
$$ y\left(k\right)=c\left(k\right)+\eta \left(k\right) $$ (4) 其中, $ \eta \left(k\right) $是时间戳获取过程中随机延迟等导致的观测偏差, 通常情况下可用一个高斯过程来描述[1].
综合式(3)和式(4), 并取${\boldsymbol{C}}=\left[ 1\;\;0 \right]$, 则时间同步的本地时钟模型可以用状态方程描述为
$$ \left\{ {\begin{aligned} &{x\left( {k + 1} \right) = {\boldsymbol{A}}x\left( k \right) + {\boldsymbol{F}}d\left( k \right)}\\ &{y\left( k \right) = {\boldsymbol{C}}x\left( k \right) + \eta \left( k \right)} \end{aligned}} \right. $$ (5) 值得指出的是, 式(5)描述的是自由时钟, 即时钟按照自身有偏差的频率进行更新, 与标准时间的偏差通常会不断变大, 因此要实现时间同步, 则需要设计负反馈控制器对本地的有偏差时钟进行调控, 以使其维持在较小的时间同步误差范围内. 时间同步控制器的作用是通过控制输入量进而调整时钟的频率或相位, 设$ u\left(k\right) $为第$ k $时刻的控制输入量, 则受控时钟模型可以描述为
$$ \left\{ {\begin{aligned} &{x\left( {k + 1} \right) = {\boldsymbol{A}}x\left( k \right) + {\boldsymbol{B}}u\left( k \right) + {\boldsymbol{F}}d\left( k \right)}\\ &{y\left( k \right) = {\boldsymbol{C}}x\left( k \right) + \eta \left( k \right)} \end{aligned}} \right. $$ (6) 实际工业物联网嵌入式系统中难以直接调整晶体振荡频率, 常见的控制策略是相位调整, 取$ {\boldsymbol{B}}= {\left[1\;0\right]}^{\mathrm{T}} $. 控制输入$ u\left(k\right) $的取值主要是根据本地时间测量值$ y\left(k\right) $与参考时钟的测量值的差值来确定, 本文提出的PI控制器的设计将在第2节详细介绍.
1.3 基于包交换协议的同步误差获取
时间同步的核心技术之一是如何获取本地时钟与参考时钟的时间差. 在工业网络分布式系统中, 参考时钟与本地时钟分布于不同节点, 需要通过特定的网络通信协议来获取参考时钟的信息. 由于现有大多数通信网络是基于包交换的数字通信网络, 本节主要讨论工业物联网中如何通过包交换协议获取时钟偏差. 在第$ k $时刻, 漂移时钟$ i $和参考时钟之间的同步误差的真值, 即时间偏差$ {\theta }_{i}^{{*}}\left(k\right) $, 定义为
$$ {\theta }_{i}^{{*}}\left(k\right)={c}_{i}\left(k\right)-{c}_{0}\left(k\right) $$ (7) 其中, $ {c}_{i}\left(k\right) $和$ {c}_{0}\left(k\right) $分别表示第$ k $时刻漂移时钟$ i $的值和主参考时钟的时间值. 结合式(6)和式(7), 第$ i $个受控时钟的偏差模型可以表达为
$$ {\theta }_{i}^{{*}}\left(k+1\right)={\theta }_{i}^{{*}}\left(k\right)+{u}_{i}\left(k\right)+{\chi }_{i}\left(k\right)h+{\omega }_{ci}\left(k\right) $$ (8) 其中, $ {\theta }_{i}^{{*}}\left(k\right) $表示第$ i $个时钟节点在第$ k $时刻与主节点的时钟相位偏差, $ {u}_{i}\left(k\right) $为第$ i $个时钟节点在$ k $时刻相位补偿量, ${\chi}_i \left(k\right)h$表示第$ i $个时钟节点在$ k $时刻由于自身时钟频率偏斜带来的扰动, $ {\omega }_{ci}\left(k\right) $为第$ i $个时钟的相位噪声.
工业物联网中的时间同步协议通常以双向包交换为主, 如$ \mathrm{I}\mathrm{E}\mathrm{E}\mathrm{E}\;1588\;\mathrm{P}\mathrm{T}\mathrm{P} $(Precision time protocol)协议. 不失一般性, 本文以$ \mathrm{P}\mathrm{T}\mathrm{P} $机制为例, 通过交换带有时间戳的同步包完成点对点时钟偏差$ \theta \left(k\right) $的获取.
如图1所示, 通过包交换过程[24], 最终得到其时钟偏差值$ {\theta }_{i}\left(k\right) $为
$$ {\theta }_{i}\left(k\right)=\frac{\left(t_2-t_1\right)-\left(t_4-t_3\right)}{2} $$ (9) 由于存在时间戳的不确定性, 时间戳$t_1,{t}_2,{t}_3$和${t}_4$中包含有观测噪声$ \eta \left(k\right) $, 更重要的是包交换通信中普遍存在的传输延迟抖动导致双向包交换不对称, 由$ \mathrm{P}\mathrm{T}\mathrm{P} $等时间同步协议获得的$ \theta \left(k\right) $是一个真值$ {\theta }^{{*}}\left(k\right) $带偏差的观测值, 由此可得第$ i $个时钟的偏差观测值为
$$ {\theta }_{i}\left(k\right)={\theta }_{i}^{{*}}\left(k\right)+{\tau }_{i}\left(k\right) $$ (10) 其中, $ {\tau }_{i}\left(k\right) $为第$ i $个时钟节点在$ k $时刻包交换通信中双向包交换不对称、时间戳不确定性等带来的扰动误差.
由于这些误差的影响, 从而会导致同步精度的恶化. 为进一步提高系统精度, 本文提出了带抗扰补偿的时间同步系统设计方案.
2. 基于抗扰补偿器的时间同步系统
对于一个点对点网络的时间同步, 最终目标就是两个时钟的时间值$ {c}_{i}\left(k\right) $和$ {c}_{0}\left(k\right) $相等, 亦即时钟同步误差${{\mathrm{l}\mathrm{i}\mathrm{m}}_{k\to {\infty }}}{\theta }_{i}^{{*}}\left(k\right)=0$. 但由于漂移时钟存在频率偏斜、相位噪声等, 影响了时间同步的精度. 为了进一步提高时间同步系统的精度和抗扰性, 满足时间敏感工业应用的需要, 本文提出了扩展PI抗扰补偿器, 主要用于减少扰动信号$ d\left(k\right) $对时间同步的影响, 尤其是针对$ \gamma \left(k\right) $中的周期性扰动$ f\left(k\right) $设计了一种新型的扩展PI补偿器, 提出了零点配置的补偿器优化设计方法. 针对$ \gamma \left(k\right) $中缓慢变化的扰动$ {\gamma }_{0} $和延时均值导致的相位测量扰动, 采用了PI控制器消除稳态误差, 并证明了带补偿器的PI控制闭环系统的稳定性条件.
本文提出的带抗扰补偿的时间同步的系统如图2所示, 通过设计虚线框中的抗扰补偿器, 并选取合适的参数, 将残差信号通过传递函数矩阵$\boldsymbol{Q}\left(z\right)$叠加到被控时钟模型的输入上, 以此实现对于被控时钟中周期性扰动信号$ f\left(k\right) $的抑制作用; 补偿后的漂移时钟可以近似看作一个无扰动的理想时钟, 然后基于理想时钟模型设计控制器实现点对点的时间同步, 并通过积分控制消除缓慢频率漂移$ {\gamma }_{0} $、时间戳不确定性、包传输延时不对称等形成的偏差.
所提出的抗扰时间同步系统可以看作是一个内外环控制的双闭环系统. 外环考虑主从两个时钟, 通过第$ k $次包交换过程得到从时钟$ i $和主时钟的偏差$ {\theta }_{i}\left(k\right) $, 偏差信号作为外环控制器的输入, 通过控制器的作用输出$ {u}_{i}^{0}\left(k\right) $调节时钟的相位, 以此达到主从时钟的同步; 内环抗扰补偿器则主要是对于时钟模型中的扰动和噪声进行一定的补偿作用, 由两部分构成, 包含一个无相位噪声、无频率偏斜和测量扰动的理想时钟, 用于描述在控制输入$ {u}_{i}^{0}\left(k\right) $的作用下, 一个无噪声无扰动的理想时钟的动态响应, 可表达为
$$ \left\{ {\begin{aligned} &{{{\hat x}_i}\left( {k + 1} \right) = {\boldsymbol{A}}{{\hat x}_i}\left( k \right) + {\boldsymbol{B}}u_i^0\left( k \right)}\\ &{{{\hat y}_i}\left( k \right) = {\boldsymbol{C}}{{\hat x}_i}\left( k \right)} \end{aligned}} \right. $$ (11) 其中, $ {u}_{i}^{0}\left(k\right) $为第$ i $个时钟节点外环控制器的输出值, 且${\hat{x}}_{i}\left(k\right)={\left[{\hat{c}}_{{i}}\left(k\right)\;{\hat{\gamma }}_{i}\left(k\right)\right]}^{\mathrm{T}}$, 其中${ \hat{c}}_{{i}}(k)$为理想时钟值, ${\hat{\gamma }}_{i}\left(k\right)$为理想时钟速率值, 即${\hat{\gamma }}_{i}\left(0\right)=1$, 等价于${\chi }_{i}\left(0\right)= 0$. 状态转移矩阵${\boldsymbol{A}}=\left[ {\begin{aligned} 1\;\;h\\ 0\;\;1 \end{aligned}} \right]$.
由于实际时钟$ {\chi }_{i}\left(0\right)\ne 0 $, 与理想时钟具有不同频率偏斜, 实际时钟还存在噪声和扰动的影响, 在$ {u}_{i}^{0}\left(k\right) $的作用下, 实际时钟的输出$ {y}_{i}\left(k\right) $与理想时钟模型的输出${\hat{y}}_{i}\left(k\right)$并不一致. 定义残差信号$ {r}_{i}\left(k\right) $用来描述实际时钟与理想时钟的偏差.
内环抗扰补偿器的第二个部分是一个动态系统, 其动态响应特性可以用一个离散域的传递函数矩阵$\left(\mathrm{T}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{f}\mathrm{e}\mathrm{r}\; \mathrm{f}\mathrm{u}\mathrm{n}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\;\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{x}, \mathrm{T}\mathrm{F}\mathrm{M}\right) \;\boldsymbol{Q}(z)$来描述. 动态系统$\boldsymbol{Q}(z)$以残差信号${r}_{i}(k)$作为输入, 生成一个补偿信号${\overline {v}}_{i}(k)$, 与原有的控制信号${u}_{i}^{0}(k)$相加后, 共同作用于有扰动的实际时钟系统. 补偿信号${\overline {v}}_{i} \left(k\right)$的主要目的是抵消扰动$ d\left(k\right) $导致时钟偏离理想时钟的作用, 使得在扰动不为零时, 实际时钟状态${x}_{i}(k)$与理想时钟模型$ {\hat{x}}_{i}\left(k\right) $仍然保持一致, 从而反向“补偿”了扰动的作用. 本文将$\boldsymbol{Q} (z)$称作补偿器反馈回路, 与传统的扰动观测器不同, 本文针对时钟频率偏斜中存在的扰动, 提出了一种扩展PI补偿回路和零极点联合配置方法, 可以在扰动的结构不清楚的情况下, 根据扰动频率即可完成设计, 实现对特定频率扰动信号的补偿消减作用, 达到更高精度的时间同步.
值得指出的是, 本文提出的如图2所示的抗扰时间同步系统, 属于一种抗扰性能和跟踪性能二自由度控制的思想[13], 即内环为抗扰补偿作用, 用于提高时间同步系统的鲁棒性; 外环为相位跟踪作用, 用于提供系统的时钟相位同步性能. 如后文所证明, 在$ k\to {\infty } $时, 所提出的带抗扰补偿器的时钟系统是稳态收敛的, 在稳态下能很好地消除扰动的影响, 其动态特性与无扰动理想时钟模型一致. 因此, 在稳态下, 控制器和补偿器可以分开独立设计, 设计控制器时, 将抗扰补偿器补偿后的时钟系统作为一个扰动为零的理想时钟系统即可. 但是考虑系统的动态特性, 抗扰补偿器和控制器的参数会存在一定的约束关系, 因此在设计控制器时, 通过综合考虑抗扰补偿器的动态特性, 对控制器参数进行优化设计, 可以在保证稳态下抗扰性能和追踪性能的基础上, 提高整个时间同步系统的动态响应特性.
3. 基于扩展PI的抗扰补偿器设计
基于前述抗扰补偿器的时间同步方案, 本文提出的采用扩展PI观测器的新型抗扰补偿器结构如图3所示.
对于某个带扰动不精确的漂移时钟$ i $, 所提出的扰动补偿器的无扰动的理想时钟模型如式(11)所示, 反馈回路$\boldsymbol{Q}\left(z\right)$是具有积分环节的动态系统. 整个抗扰补偿器可以视作一个$ {u}_{i}^{0}\left(k\right) $ 和$ {y}_{i}\left(k\right) $ 为输入、$ {v}_{i}\left(k\right) $为输出的动态系统. 设$ {r}_{i}\left(k\right) $为理想时钟模型的输出值${\hat{y}}_{i}\left(k\right)$与实际时间戳测量值$ {y}_{i}\left(k\right) $的残差信号, 表示为
$$ {r}_{i}\left(k\right)={\hat{y}}_{i}\left(k\right)-{y}_{i}\left(k\right) $$ (12) 反馈回路$\boldsymbol{Q}\left(z\right)$是由${{\boldsymbol{K}}}_{{{1}}},{{\boldsymbol{K}}}_{{{2}}},{{\boldsymbol{K}}}_{{{3}}},{{\boldsymbol{K}}}_{{{4}}}$参数矩阵所构成的一个动态系统, 可表示为
$$ \left\{ {\begin{aligned} &{{z_i}\left( {k + 1} \right) = {{\boldsymbol{K}}_{{1}}}{z_i}\left( k \right) + {{\boldsymbol{K}}_{{2}}}{r_i}\left( k \right)}\\ &{{v_i}\left( k \right) = {{\boldsymbol{K}}_{{3}}}{z_i}\left( k \right) + {{\boldsymbol{K}}_{{4}}}{r_i}\left( k \right)} \end{aligned}} \right. $$ (13) 其中, $ {z}_{i}\left(k\right) $为抗扰补偿器反馈回路的扩展状态变量. 该动态反馈回路是以残差信号$ {r}_{i}\left(k\right) $为输入、以$ {v}_{i}\left(k\right) $为输出的动态系统, $ {v}_{i}\left(k\right) $通过一个${{\boldsymbol{B}}}_{{c}}$矩阵作为输入信号, 调节系统状态值$ x\left(k\right) $, 使得系统输出值$ {y}_{i}\left(k\right) $与理想输出${\hat{y}}_{i}\left(k\right)$之间的残差$ {r}_{i}\left(k\right) $最小, 系统输出接近理想输出. 因为$ {\boldsymbol{B}} $和$ v $都为$ m\times n $的矩阵, 且$ {\boldsymbol{B}} $为列满秩, 因此存在$ n\times m $的伪逆矩阵${{\boldsymbol{B}}}_{{{c}}}$使得${\boldsymbol{B}}{{\boldsymbol{B}}}_{{{c}}}={{\boldsymbol{I}}}_{m{\times}m}$.
在补偿器输出的补偿量${{\boldsymbol{B}}}_{{{c}}}{v}_{i}\left(k\right)$的作用下, 实际作用于有扰漂移时钟上的控制输入为${u}_{i}^{0}\left(k\right)+ {{\boldsymbol{B}}}_{{{c}}}{v}_{i}\left(k\right)$, 因为${\boldsymbol{B}}{{\boldsymbol{B}}}_{{{c}}}$为单位阵, 代入时钟模型(6), 可得增加扩展PI抗扰补偿器后的被控对象的状态空间方程为
$$ \left\{ {\begin{split} &{{x_i}\left( {k + 1} \right) = {\boldsymbol{A}}{x_i}\left( k \right) + {\boldsymbol{B}}u_i^0\left( k \right) + {\boldsymbol{F}}{d_i}\left( k \right) + {v_i}\left( k \right)}\\ &{{y_i}\left( k \right) = {\boldsymbol{C}}{x_i}\left( k \right) + {\eta _i}\left( k \right)} \end{split}} \right. $$ (14) 抗扰补偿器反馈作用补偿值$ {v}_{i}\left(k\right) $主要是受待设计的扩展PI抗扰补偿器的4个反馈参数矩阵${{\boldsymbol{K}}}_{{{1}}},{{\boldsymbol{K}}}_{{{2}}},{{\boldsymbol{K}}}_{{{3}}},{{\boldsymbol{K}}}_{{{4}}}$决定, 是影响扰动补偿和同步精度的关键矩阵. 下面重点讨论所提出的扩展PI补偿器的收敛性和参数矩阵的优化设计方法.
3.1 抗扰补偿器稳定性分析
对时间同步系统, 本文提出的如式(13)和式(14)所示的扩展PI抗扰补偿器具有如下的稳定性定理.
定理 1. 对于受到特定周期性扰动$ d\left(k\right) $的时间同步系统$ \left(6\right)和\left(10\right) $, 以及对应的如式(13)和式(14)的扩展$ \mathrm{P}\mathrm{I} $抗扰补偿器, 当且仅当$\left[ {\begin{align} {\boldsymbol{A}}-{{\boldsymbol{K}}}_{{{4}}}{\boldsymbol{C}}\;\; -{{\boldsymbol{K}}}_{{{3}}}\\ {{\boldsymbol{K}}}_{{{2}}}{\boldsymbol{C}}\;\;\;\;\;\; {{\boldsymbol{K}}}_{{{1}}}\;\; \end{align}} \right]$的特征值在$ {\bf{Z}} $域的单位圆内, 该扩展PI抗扰补偿器是稳定的, 即该抗扰补偿器的输出误差${e}_{i}\left(k\right)= {\hat{x}}_{i}\left(k\right)- {x}_{i}\left(k\right)$可以指数收敛${\rm{li}}{{\rm{m}}_{{{k}} \to \infty }}{e}_{i}\left(k\right)=0$. 带补偿的时钟状态$ {x}_{i}\left(k\right) $可以指数收敛于理想时钟模型状态${\hat{x}}_{i}\left(k\right)$, ${\rm{li}}{{\rm{m}}_{{{k}} \to \infty }}{{\rm{E}}(x}_{i}\left(k\right))={{\rm{E}}(\hat{x}}_{i}\left(k\right)$), ${\rm{E}}$表示时间平均算子的期望, 即带补偿器的有扰时钟与理想时钟模型可以在时间平均指数渐进意义下保持一致.
证明. 令带补偿的实际时钟与理想时钟模型状态误差${e}_{i}\left(k\right)={\hat{x}}_{i}\left(k\right)-{x}_{i}\left(k\right)$, 将系统理想模型$ \left(11\right) $和实际系统(14)代入定义式$ \left(12\right) $中, 可得
$$ {r}_{i}\left(k\right)={\boldsymbol{C}}{e}_{i}\left(k\right)-{\eta }_{i}\left(k\right) $$ (15) 将式$\left(13\right)$和式$\left(15\right)$代入式$ \left(14\right) $, 可得
$$ \begin{split} {x}_{i}\left(k+1\right)=\;&{\boldsymbol{A}}{x}_{i}\left(k\right)+{\boldsymbol{B}}{u}_{i}^{0}\left(k\right)+{{\boldsymbol{K}}}_{{{3}}}{z}_{i}\left(k\right)+\\ &{{\boldsymbol{K}}}_{{{4}}}\left({\boldsymbol{C}}{e}_{i}\left(k\right)-{\eta }_{i}\left(k\right)\right)+{\boldsymbol{F}}{d}_{i}\left(k\right) \end{split} $$ (16) 将式$ \left(15\right) $代入式$ \left(13\right) $, 可得
$$ {z}_{i}\left(k+1\right)={{\boldsymbol{K}}}_{{{1}}}{z}_{i}\left(k\right)+{{\boldsymbol{K}}}_{{{2}}}\left({\boldsymbol{C}}{e}_{i}\left(k\right)-{\eta }_{i}\left(k\right)\right) $$ (17) 由于${e}_{i}\left(k\right)={\hat{x}}_{i}\left(k\right)-{x}_{i}\left(k\right)$, 结合式$ \left(16\right) $和理想模型$ \left(11\right) $, 可得
$$ \begin{split} {e}_{i}\left(k+1\right)=\;&\left({\boldsymbol{A}}-{{\boldsymbol{K}}}_{{{4}}}{\boldsymbol{C}}\right){e}_{i}\left(k\right)-{{\boldsymbol{K}}}_{{{3}}}{z}_{i}\left(k\right) +\\ &{{\boldsymbol{K}}}_{{{4}}}{\eta }_{i}\left(k\right)+{\boldsymbol{F}}{d}_{i}\left(k\right) \end{split} $$ (18) 联立式$\left(15\right),\;\left(17\right)和\left(18\right)$, 由实际时钟和扩展PI补偿器构成的整体系统可以用增广状态空间方程描述
$$ \left\{ {\begin{aligned} &\left(\begin{array}{c}{e}_{i}\left(k+1\right)\\ {z}_{i}\left(k+1\right)\end{array}\right)={\widetilde {\boldsymbol{A}}}\left(\begin{array}{c}{e}_{i}\left(k\right)\\ {z}_{i}\left(k\right)\end{array}\right)+{\widetilde {\boldsymbol{D}}}_{{{d}}}{\eta }_{i}\left(k\right)+\\ &\qquad\qquad\qquad\qquad\left(\begin{array}{c}{\boldsymbol{F}}{d}_{i}\left(k\right)\\ {{\bf{0}}}_{2\times 1}\end{array}\right)\\ &{r}_{i}\left(k\right)={\widetilde {\boldsymbol{C}}}\left(\begin{array}{c}{e}_{i}\left(k\right)\\ {z}_{i}\left(k\right)\end{array}\right)-{\eta }_{i}\left(k\right) \end{aligned}} \right. $$ (19) 其中, ${\widetilde {\boldsymbol{A}}}= \left[ {\begin{align} {\boldsymbol{A}} - {{\boldsymbol{K}}}_{{{4}}}{\boldsymbol{C}}\; -{{\boldsymbol{K}}}_{{{3}}}\\ {{\boldsymbol{K}}}_{{{2}}}{\boldsymbol{C}}\;\;\;\;\;\;\; {{\boldsymbol{K}}}_{{{1}}}\;\; \end{align}} \right] ,$ ${\widetilde {\boldsymbol{C}}}=\left[ {\boldsymbol{C}}\;\; {{\bf{0}}}_{1\times 2} \right]$, ${\widetilde {\boldsymbol{D}}}_{{{d}}}= {\left[ {{\boldsymbol{K}}}_{{{4}}};\;\; -{{\boldsymbol{K}}}_{{{2}}} \right]}^{\mathrm{T}}$.
不难看出, 对于补偿器状态$ {z}_{i}\left(k\right) $和估计误差$ {e}_{i}\left(k\right) $所构成的闭环动态系统$ \left(19\right) $, 当状态转移矩阵$\widetilde {{\boldsymbol{A}}}$的特征根在单位圆内, 系统是稳定的, 即无扰动$ {d}_{i}\left(k\right)=0 $时${\rm{li}}{{\rm{m}}_{{{k}} \to \infty }}{e}_{i}\left(k\right)=0$. 当存在随机噪声时, 即$ {d}_{i}\left(k\right)\ne 0 $, 因为随机噪声的数学期望为0且相互独立, 故系统误差$ {e}_{i}\left(k\right) $在统计平均的意义下是收敛的, 即
$$ \mathop {{\rm{lim}}}\limits_{{{k}} \to \infty } {{{\rm{E}}}}\left({x}_{i}\left(k\right)\right)={{{\rm{E}}}}\left({\hat{x}}_{i}\left(k\right)\right) $$ □ 上述结论在不考虑噪声和扰动时是正确的, 当仅考虑随机噪声在统计平均的意义下是一致的情况时, 如果存在非零均值扰动, 则理论上还不完善.
对于一个存在特定频率扰动的工业网络时间同步系统, 可以构造一个扩展PI抗扰补偿器, 该抗扰补偿器的零点由开环系统的零点和矩阵${{\boldsymbol{K}}}_{{{1}}}$的特征值组成. 即通过配置参数矩阵${{\boldsymbol{K}}}_{{{1}}}$便可以配置相应的零点到扰动信号的特征频率附近, 利用零极点对消的原理使传递函数矩阵降秩[21], 从而衰减扰动信号的影响.
3.2 抗扰补偿器零极点优化设计
抗扰补偿器的目的是为了使被控系统中即使存在扰动$ d $时, 其实际输出$ y $与期望的无扰动系统的输出仍保持一致. 即在$ d\ne 0 $时, 仍然具有$y\left(t\right)= \hat{y}\left(t\right)$, 也就是使得受扰系统的输入−输出的动力学特性符合理想的无扰动的动力学模型, 从而使得针对无扰动模型所设计的控制器保持良好的跟踪性能. 因此, 抗扰补偿器优化设计的目的就是通过选择合适的参数矩阵${{\boldsymbol{K}}}_{{{1}}},{{\boldsymbol{K}}}_{{{2}}},{{\boldsymbol{K}}}_{{{3}}},{{\boldsymbol{K}}}_{{{4}}}$, 使得在$ d\ne 0 $时, 残差信号$ {r}_{i}\left(k\right) $最小.
针对上述优化问题, 结合时间同步中存在特定频率的频率扰动, 本文提出了一种零点和极点联合配置的方法, 优化系统性能, 达到消除扰动的干扰, 使得添加了补偿器的动态系统(11), (13)和(14)尽可能接近无扰动系统的理想输出.
对于如式$ \left(19\right) $的系统方程, $ {G}_{d1}\left(z\right) $定义为扰动信号$ {\eta }_{i}\left(k\right) $到残差输出$ {r}_{i}\left(k\right) $的传递函数矩阵
$$ {G}_{d1}\left(z\right)={\widetilde {\boldsymbol{C}}}{\left(z{\boldsymbol{I}}-{\widetilde {\boldsymbol{A}}}\right)}^{-1}{\widetilde {\boldsymbol{D}}}_{{{d}}}-1 $$ (20) 定义$ {G}_{d2}\left(z\right) $为扰动信号$ d\left(k\right) $到残差输出$ {r}_{i}\left(k\right) $的传递函数矩阵
$$ {G}_{d2}\left(z\right)={\widetilde {\boldsymbol{C}}}{\left(z{\boldsymbol{I}}-{\widetilde {\boldsymbol{A}}}\right)}^{-1}\widetilde {{\boldsymbol{F}}} $$ (21) 其中, $\widetilde {{\boldsymbol{F}}}={\left[ {\bf{1}}\;\;{\bf{1}} \;\;{{{\bf{0}}_{2 \times 1}}} \right]^{\rm{T}}}$为式(19)中相应维度的状态扰动输入矩阵.
考虑被控系统所有的扰动信号对残差的影响, 取扰动传递函数矩阵$ {G}_{d}\left(z\right)={G}_{d1}\left(z\right)+{G}_{d2}\left(z\right) $.
1) 零点配置: 考虑测量扰动主要集中在有限频率的周期信号, 如果将零点配置在扰动频率处, 那么就可以有效地使得扰动信号进行衰减, 使得系统输出更接近无扰动系统的理想输出. 对于采样周期为$ T $的离散系统, 如果存在一个干扰频率为$ {\omega }_{i} $的时钟频率干扰信号, 则矩阵${{\boldsymbol{K}}}_{1}^{i,i}$可以构造为
$$ {{\boldsymbol{K}}}_{{{1}}}^{{{i}},{{i}}}=\left[ {\begin{array}{*{20}{c}}\mathrm{cos}\left({\omega }_{i}T\right)& \mathrm{sin}\left({\omega }_{i}T\right)\\ -\mathrm{sin}\left({\omega }_{i}T\right)& \mathrm{cos}\left({\omega }_{i}T\right) \end{array}} \right] $$ (22) 不难求出, 该矩阵的特征根为${{\rho }_{i}=\mathrm{e}}^{\pm {\rm{j}}{\omega }_{i}T}$. 如果时钟系统中存在$ n $个不同频率的扰动信号, 则${{\boldsymbol{K}}}_{{{1}}}$可以构造成一个$ 2n\times 2n $的分块矩阵
$$ {{\boldsymbol{K}}}_{{{1}}}={\rm{diag}}\left\{{\boldsymbol{K}}_{{{1}}}^{{{1}},{{1}}},{{\boldsymbol{K}}}_{{{1}}}^{{{2}},{{2}}}, \cdots ,{{\boldsymbol{K}}}_{{{1}}}^{{{i}},{{i}}}\right\} $$ (23) 因此, 矩阵${{\boldsymbol{K}}}_{{{1}}}$的特征根是各分块矩阵${{\boldsymbol{K}}}_{{{1}}}^{{{i}},{{i}}}$特征根的集合$ \left\{{\rho }_{i}\right\} $. 根据系统零点使得闭环系统的传递函数矩阵(TFM)降秩的原理[21], 通过上述设计, 将补偿器闭环系统的零点配置在扰动频率$ \left\{{\omega }_{i}\right\} $处, 利用在$ \left\{{\omega }_{i}\right\} $处由扰动信号到残差信号的TFM发生降秩, 从而可以有效削弱扰动信号对残差信号的影响, 有助于更好地达到实际系统的输出与无扰动系统的输出一致的目标.
2) 极点配置: 在利用零点配置确定反馈矩阵${{\boldsymbol{K}}}_{{{1}}}$后, 下一步通过极点配置和函数优化方法, 选择合适的${{\boldsymbol{K}}}_{{{2}}},{{\boldsymbol{K}}}_{{{3}}},{{\boldsymbol{K}}}_{{{4}}}$取值, 达到保证补偿器系统的稳定性, 并同时优化系统性能的目的, 即状态转移矩阵的特征根在单位圆内, 并同时优化参数选择${{\boldsymbol{K}}}_{{{2}}}, {{\boldsymbol{K}}}_{{{3}}}, {{\boldsymbol{K}}}_{{{4}}}$, 使得扰动信号对残差信号$ {r}_{i}\left(k\right) $的影响最小化.
综合考虑在系统稳定和相应参数矩阵关系约束的条件下, 可以得到测量扰动频率${\omega}_{i}$信号对$ {r}_{i}\left(k\right) $的影响, 并使得其对残差$ {r}_{i}\left(k\right) $的影响最小化, 可以得到$ {r}_{i}\left(k\right) $对于扰动信号的鲁棒性目标函数为
$$ \begin{split} & \mathop {{\rm{min}} }_{{{\boldsymbol{K}}_{{2}}},{{\boldsymbol{K}}_{{3}}},{{\boldsymbol{K}}_{{4}}}}{J}_{d}=\sum _{i=1}^{n}{m}_{i}{||{G}_{d}\left(z\right)||}_{2,z={{\rm{e}}}^{{\rm{j}}{\omega }_{i}T}} \\ &\;\mathrm{s}.\mathrm{t}. \;\;\;\lambda \left({\widetilde {\boldsymbol{A}}}\right)\le 1\\& \qquad {a}_{{K}_{3}}{a}_{{K}_{2}}+{b}_{{K}_{3}}{b}_{{K}_{2}}=0,{b}_{{k}_{1}}{b}_{{k}_{2}}=0,{b}_{{k}_{1}}{a}_{{k}_{2}}=0 \end{split} $$ (24) 其中, $ {m}_{i} $为权重, 其参数取值是通过频谱分析得到的相应频率扰动的幅值大小来选择. 幅值越大, 对应的$ {m}_{i} $权重分量也越大, 从而使得相应频率的扰动信号对输出的影响越小. 同时考虑约束条件使得系统稳定, 即状态转移矩阵$ {\widetilde {\boldsymbol{A}}} $的特征根在单位圆内, ${a}_{{K}_{3}}{a}_{{K}_{2}}+{b}_{{K}_{3}}{b}_{{K}_{2}}=0 ,{b}_{{k}_{1}}{b}_{{k}_{2}}=0,{b}_{{k}_{1}}{a}_{{k}_{2}}=0$为${{\boldsymbol{K}}}_{{{2}}},{{\boldsymbol{K}}}_{{{3}}}$矩阵的参数约束, 参数矩阵的具体参数约束设定详见定理2的时间同步控制器参数设计.
4. 时间同步控制器设计
由于点对点网络时间同步最终的控制目标是时钟偏差真值${{\rm{lim}}}_{k \to \infty } {\theta }_{i}^{{*}}\left(k\right)=0$, 且实际包交换获取时钟偏差的过程中会受到包交换通讯不对称影响, 从而带来时钟扰动偏差, 因此要实现点对点时钟同步必须要设计合理的控制器, 起到相位调节作用, 实现时钟相位跟踪效果, 达到时钟同步.
为更好地消除稳态下由于包交换延时不对称导致的稳态误差, 本文选取PI控制器作为同步系统的外环控制器, 以提高系统的跟踪性能. PI控制器以包交换获得的时钟偏差$ {\theta }_{i}\left(k\right) $作为输入信号, 其控制方程可以描述为
$$ \left\{ {\begin{aligned} &{{\omega_i}\left( {k + 1} \right) = {\omega_i}\left( k \right) + \beta {\theta _i}\left( k \right)}\\ &{u_i^0\left( k \right) = {\omega_i}\left( k \right) + \alpha {\theta _i}\left( k \right)} \end{aligned}} \right. $$ (25) 其中, ${\omega}_{i}\left(k\right)$为PI控制器的内部状态变量, 体现积分控制器的作用.
考虑包交换过程的噪声, 则PI控制器的方程可以改写为
$$ \left\{ {\begin{aligned} &{{\omega_i}\left( {k + 1} \right) = {\omega_i}\left( k \right) + \beta \left( {{\theta _i}\left( k \right) + {\tau _i}\left( k \right)} \right)}\\ &{u_i^0\left( k \right) = {\omega_i}\left( k \right) + \alpha \left( {{\theta _i}\left( k \right) + {\tau _i}\left( k \right)} \right)} \end{aligned}} \right. $$ (26) 且由于$ {r}_{i}\left(k\right)=-{\theta }_{i}\left(k\right) $, 为便于分析, 将抗扰补偿器反馈回路$ \left(13\right) $中的$ {r}_{i}\left(k\right) $替换为$ -{\theta }_{i}\left(k\right) $, 即
$$ \left\{ {\begin{aligned} &{z}_{i}\left(k+1\right)={{\boldsymbol{K}}}_{{{1}}}{z}_{i}\left(k\right)-{{\boldsymbol{K}}}_{{{2}}}{\theta }_{i}\left(k\right)\\ &{v}_{i}\left(k\right)={{\boldsymbol{K}}}_{{{3}}}{z}_{i}\left(k\right)-{{\boldsymbol{K}}}_{{{4}}}{\theta }_{i}\left(k\right) \end{aligned}} \right. $$ (27) 受控时钟偏差模型(8)代入PI控制器方程, 得
$$ \begin{split} {\theta }_{i}\left(k+1\right)=\;&\,{\theta }_{i}\left(k\right)+{\chi }_{i}\left(k\right)h+{\omega}_{i}\left(k\right) +\\ &\,\alpha ({\theta }_{i}\left(k\right)+{\tau }_{i}\left(k\right))+{\boldsymbol{C}}\left({{\boldsymbol{K}}}_{{{3}}}{z}_{i}\left(k\right)\right.-\\ &\left. {{\boldsymbol{K}}}_{{{4}}}{\theta }_{i}\left(k\right)\right)+{\omega }_{ci}(k) \end{split} $$ (28) 定义状态变量, 得到状态空间模型为
$$ \begin{split} \left[ {\begin{array}{*{20}{c}} {{\theta _i}(k + 1)}\\ {{\omega_i}(k + 1)}\\ {{z_i}(k + 1)} \end{array}} \right] = \;&\left[ {\begin{array}{*{20}{c}} {1 + \alpha - {\boldsymbol{C}}{{\boldsymbol{K}}_{{4}}}}&1&{{\boldsymbol{C}}{{\boldsymbol{K}}_{{3}}}}\\ \beta &1&{\boldsymbol{0}}\\ { - {{\boldsymbol{K}}_2}}&{\boldsymbol{0}}&{{{\boldsymbol{K}}_{{1}}}} \end{array}} \right]\times\\ &\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{\theta _i}\left( k \right)}\\ {{\omega_i}\left( k \right)} \end{array}}\\ {{z_i}\left( k \right)} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} 1&{ - \alpha }&{\boldsymbol{0}}\\ {{\boldsymbol{0}}}&{ - \beta }&{\boldsymbol{0}}\\ {\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}} \end{array}} \right]\times\\ &\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{\chi _i}\left( k \right)h + {\omega _{ci}}\left( k \right)}\\ {{\tau _i}\left( k \right)} \end{array}}\\ {\boldsymbol{0}} \end{array}} \right] \\[-15pt]\end{split} $$ (29) 其中, $ {z}_{i}\left(k\right) $为2 × 1状态变量; ${{\boldsymbol{K}}}_{{{2}}}$、${{\boldsymbol{K}}}_{{{4}}}$为2 × 1矩阵; ${{\boldsymbol{K}}}_{{{1}}}$、${{\boldsymbol{K}}}_{{{3}}}$为2 × 2矩阵.
取状态变量${p}_{i}\left(k\right)={\left[{\theta }_{i}\left(k\right)\;{\omega}_{i}\left(k\right)\;{z}_{i}\left(k\right)\right]}^{\mathrm{T}}$, 那么式(29)可以化简为
$$ {p}_{i}\left(k+1\right)={\boldsymbol{G}}{p}_{i}\left(k\right)+{\boldsymbol{H}}{o}_{i}\left(k\right) $$ (30) 其中, ${\boldsymbol{G}} = \left[ {\begin{aligned} {1 + \alpha - {\boldsymbol{C}}{{\boldsymbol{K}}_{{4}}}}\;\;1\;\;{{\boldsymbol{C}}{{\boldsymbol{K}}_{{3}}}}\\ \beta \;\;\;\;\;\;\;\;\;\;\;\;{1}\;\;\;\;\,{\boldsymbol{0}}\;\ \;\;\\ { - {{\boldsymbol{K}}_2}}\;\;\;\;\;\;\;\;\;{\boldsymbol{0}}\;\;\;{{{\boldsymbol{K}}_{{1}}}}\;\; \end{aligned}} \right]$为一个$ 4\times 4 $的矩阵; ${o_i}\left( k \right) = {\left[ {{\chi _i}\left( k \right)h + {\omega _{ci}}\left( k \right)}\;\;{{\tau _i}\left( k \right)} \;\;{\boldsymbol{0}} \right]^{\rm{T}}}$.
为分析控制系统的稳定性, 先给出引理1.
引理 1[25]. 若$ {\boldsymbol{A}},{\boldsymbol{B}},{\boldsymbol{C}},{\boldsymbol{D}} $是$ n\times n $矩阵, 且${\boldsymbol{C}}{\boldsymbol{D}}={\boldsymbol{DC}} $, 则有分块矩阵的行列式, 可计算为
$$ \left| {\begin{aligned} {\boldsymbol{A}}\;\; {\boldsymbol{B}}\\ {\boldsymbol{C}}\;\; {\boldsymbol{D}} \end{aligned}} \right| = |{\boldsymbol{A}}{\boldsymbol{D}}- {\boldsymbol{B}}{\boldsymbol{C}}|$$ 为了更好地说明抗扰补偿器参数矩阵${{\boldsymbol{K}}}_{{{1}}},{{\boldsymbol{K}}}_{{{2}}}, {{\boldsymbol{K}}}_{{{3}}},{{\boldsymbol{K}}}_{{{4}}}$的参数设定, 特别地, 对于单频扰动即${{\boldsymbol{K}}}_{{{1}}}$为2$\times $2的矩阵, 定义
$$ \begin{split} &{{\boldsymbol{K}}_{{i}}} = \left[ {\begin{array}{*{20}{c}} {{a_{{K_i}}}}&{{b_{{K_i}}}}\\ {{c_{{K_i}}}}&{{d_{{K_i}}}} \end{array}} \right],{{\boldsymbol{K}}_{{j}}} = \left[ {\begin{array}{*{20}{c}} {{a_{{K_j}}}}\\ {{c_{{K_j}}}} \end{array}} \right],\\ &\qquad \qquad \qquad \quad \quad i = 1,3;\;\;j = 2,4\end{split}$$ 对于本文提出的带扰动补偿的PI控制器, 有充分条件如下.
定理 2. 对于如图3所示带PI抗扰补偿器的时间同步系统, 采用如式(22)和式(24)的零极点优化的方法设计补偿器增益, 当满足以下两个条件时
1)${{\boldsymbol{K}}}_{{{2}}},{{\boldsymbol{K}}}_{{{3}}}$矩阵满足 ${a}_{{K}_{3}}{a}_{{K}_{2}}+{b}_{{K}_{3}}{b}_{{K}_{2}}=0,$${b}_{{k}_{1}}{b}_{{k}_{2}}= 0,{b}_{{k}_{1}}{a}_{{k}_{2}}=0 $
2)控制器增益$ \alpha $和$ \beta $的取值范围为
$$ \left\{ {\begin{aligned} &{0 < \alpha < {a_{{K_4}}} - 4}\\ &{0 < \beta < \frac{{{{\left( {{a_{{K_4}}} - \alpha } \right)}^2}}}{4}} \end{aligned}} \right.\qquad\quad\;\; $$ 或
$$ \left\{ {\begin{aligned} &{0 < \alpha < {a_{{K_4}}} - 4}\\ &{\frac{{{{\left( {{a_{{K_4}}} - \alpha } \right)}^2}}}{4} < \beta < {a_{{K_4}}} - \alpha } \end{aligned}} \right. $$ 则所设计的时间同步系统是渐进稳定的.
证明. 对于如式(30)所示的状态空间方程, 其状态转移矩阵$ {\boldsymbol{G}} $可以设定为分块矩阵, 即
$$ {\boldsymbol{G}}=\left[ {\begin{array}{*{20}{c}} {\boldsymbol{A}}& {\boldsymbol{B}}\\ {\boldsymbol{C}}& {\boldsymbol{D}} \end{array}} \right] $$ (31) 其中, ${\boldsymbol{A}}=\left[ {\begin{aligned} 1+\alpha -{a}_{{K}_{4}}\;\;1\\ \beta \;\;\;\;\;\qquad 1 \end{aligned}} \right]$, ${\boldsymbol{B}}=\left[ {\begin{aligned}{a}_{{K}_{3}}\;\; {b}_{{K}_{3}}\\ 0\;\;\;\;\;\; 0\;\; \end{aligned}} \right]$, ${\boldsymbol{C}}= \left[ {\begin{aligned} -{a}_{{K}_{2}}\;\; 0\\ -{b}_{{K}_{2}}\;\; 0 \end{aligned}} \right] ,$ 特别地, 由于${{\boldsymbol{K}}}_{{{1}}}$矩阵有${a}_{{K}_{1}}={d}_{{K}_{1}}, {c}_{{K}_{1}}= -{b}_{{K}_{1}}$, 因此, ${\boldsymbol{D}}=\left[ {\begin{aligned} {a}_{{K}_{1}}\;\;\;\; {b}_{{K}_{1}}\\ -{b}_{{K}_{1}}\;\; {a}_{{K}_{1}} \end{aligned}} \right] .$
根据引理1可得, 当${b}_{{k}_{1}}{b}_{{k}_{2}}= 0,{b}_{{k}_{1}}{a}_{{k}_{2}}=0$时, 闭环系统(30)的特征多项式为
$$ \begin{split} \left|z{\boldsymbol{I}}-{\boldsymbol{G}}\right|=\;&\left|\left(z{\boldsymbol{I}}-{\boldsymbol{A}}\right)\left(z{\boldsymbol{I}}-{\boldsymbol{D}}\right)-{\boldsymbol{B}}{\boldsymbol{C}}\right|= \\ \;& \left| \left[ {\begin{array}{*{20}{c}} {z - 1 - \alpha + {a_{{K_4}}}}&{ - 1}\\ { - \beta }&{z - 1} \end{array}} \right] \times \right.\\ &\left[ {\begin{array}{*{20}{c}} {z - {a_{{K_1}}}}&{ - {b_{{K_1}}}}\\ {{b_{{K_1}}}}&{z - {a_{{K_1}}}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{a_{{K_3}}}}&{{b_{{K_3}}}}\\ 0&0 \end{array}} \right] \times \\ &\left.\left[ {\begin{array}{*{20}{c}} { - {a_{{K_2}}}}&0\\ { - {b_{{K_2}}}}&0 \end{array}} \right] \right| = \left| {\begin{array}{*{20}{c}} {g_1}&{g_2}\\ {g_3}&{g_4} \end{array}} \right|\\[-15pt] \end{split} $$ (32) 其中,
$$ \begin{split} g_1=&\left(z-1-\alpha +{a}_{{K}_{4}}\right)-\\&{b}_{{K}_{1}}+ {a}_{{K}_{3}}{a}_{{K}_{2}}+ {b}_{{K}_{3}}{b}_{{K}_{2}} \\ g_2=&-{b}_{{K}_{1}}\left(z-1-\alpha +{a}_{{K}_{4}}\right)-(z-{a}_{{K}_{1}}) \\ g_3=&-\beta \left(z-{a}_{{K}_{1}}\right)+{b}_{{K}_{1}}(z-1)\\ g_4= &\;{b}_{{K}_{1}}\beta +\left(z- 1\right)\left(z-{a}_{{K}_{1}}\right) \end{split} $$ 展开式(32)并整理, 可得特征多项式方程为
$$ \begin{array}{l}\left({a}_{{K}_{3}}{a}_{{K}_{2}}+{b}_{{K}_{3}}{b}_{{K}_{2}}\right)\left(\left(z-1\right)\left(z-{a}_{{K}_{1}}\right)+{b}_{{K}_{1}}\beta \right)+ \\ \quad({z}^{2}-2{a}_{{K}_{1}}z+1)\left(\left(z-1-\alpha {a}_{{K}_{4}}\right)\left(z-1\right)+\beta \right)=0 \end{array} $$ 由上所示, 可知当条件$ {a}_{{K}_{3}}{a}_{{K}_{2}}+{b}_{{K}_{3}}{b}_{{K}_{2}}=0 $满足时, 该特征方程的解由两部分组成, 即
$$ \left({z}^{2}-2{a}_{{K}_{1}}z+1\right)=0\qquad\qquad\quad\;\; $$ (33) $$ \left(z-1-\alpha +{a}_{{K}_{4}}\right)\left(z-1\right)+\beta =0 $$ (34) 求解式(33)可知, $ {a}_{{K}_{1}}^{2} < 1 $恒成立, 方程不存在实数解, 因此控制器参数的稳点域由式(34)确定, 求解式(34)可得控制器参数的取值范围为
$$\tag{35a} \left\{ {\begin{aligned} &{0 < \alpha < {a_{{K_4}}} - 4}\\ &{0 < \beta < \frac{{{{\left( {{a_{{K_4}}} - \alpha } \right)}^2}}}{4}} \end{aligned}} \right.\qquad\quad\;\; $$ 或
$$ \tag{35b} \left\{ {\begin{aligned} &{0 < \alpha < {a_{{K_4}}} - 4}\\ &{\frac{{{{\left( {{a_{{K_4}}} - \alpha } \right)}^2}}}{4} < \beta < {a_{{K_4}}} - \alpha } \end{aligned}} \right. $$ 对于式(30)所示的系统, 其脉冲传递函数可以表达为
$$ {p}_{i}\left(z\right)={\left(z{\boldsymbol{I}}-{\boldsymbol{G}}\right)}^{-1}\left(0-{\boldsymbol{H}}{o}_{i}\left(z\right)\right) $$ (36) 由终值定理可以求得
$$ \begin{split} &\mathop {{\rm{lim}}}\limits_{k \to \infty } {p_i}\left( k \right) = \mathop {{\rm{lim}}}\limits_{z \to 1} \left( {z - 1} \right){p_i}\left( z \right)=\\ &\quad\;\mathop {{\rm{lim}}}\limits_{z \to 1} \left(z-1\right){\left(z{\boldsymbol{I}}-{\boldsymbol{G}}\right)}^{-1}\left(0-{\boldsymbol{H}}{o}_{i}\left(k\right)\right)= {\bf{0}} \end{split} $$ (37) 由此可得, 其时钟偏差最终可以收敛为0, 且抗扰补偿器和外部控制器的扩展状态最终都可以稳定到0.
□ 定理2中, ${{\boldsymbol{K}}}_{{{2}}},{{\boldsymbol{K}}}_{{{3}}}$满足${a}_{{K}_{3}}{a}_{{K}_{2}}+{b}_{{K}_{3}}{b}_{{K}_{2}}=0, {b}_{{k}_{1}}{b}_{{k}_{2}}=0,{b}_{{k}_{1}}{a}_{{k}_{2}}=0$是稳定性的一个充分条件, 可以在设计抗扰补偿器时予以考虑, 详细可参考第3.2节抗扰补偿器的极点配置(式(24)), 以此实现基于抗扰补偿器的时间同步系统设计和稳定性分析.
5. 仿真验证
为了验证本文提出的带抗扰补偿的时间同步系统的性能, 本节的实验基于硬件实验平台获取的晶振时钟参数, 搭建仿真系统进行验证. 所选硬件平台为Atmel SMART SAM R21[1], 根据硬件实测数据, 包交换过程存在过程延时约为500 μs, 过程噪声$ \tau \left(k\right) $标准差为$ {\sigma }_{\tau }=0.29\times {10}^{-6} $, 设定初始相位偏差${\theta }_{0}=500$μs, 初始时钟偏斜$ {\gamma }_{0}=20\;\mathrm{p}\mathrm{p}\mathrm{m} $, 并考虑相位噪声$ {\omega }_{c}\left(k\right) $和测量噪声$ \eta \left(k\right) $标准差为${\sigma }_{c}= {\sigma }_{\eta }=1\times {10}^{-6}$. 根据Schriegel等[26]基于晶振振动实验, 在5 g重力加速度下时钟晶振的频率漂移约为3$ \mathrm{p}\mathrm{p}\mathrm{m} $, 且考虑工业现场振动多为低频振动, 式$\left({2\rm{b}}\right)$中频率漂移噪声${f'}(k)$设定为$f=0.1\;\mathrm{H}{\rm{z}}$、幅值为$ 3\;\mathrm{p}\mathrm{p}\mathrm{m} $的正弦周期信号, 并考虑热噪声${\xi '}(k)$为$ {\sigma }_{\gamma }= 0.29\times {10}^{-6} $ 的随机噪声, 因此, 式$\left(2{\rm{b}}\right)$可以表达为
$$ {\omega }_{\gamma }\left(kT\right)=\left(3\mathrm{sin}\left(0.2\pi kT\right)+\omega \right){\rm{ppm}} $$ (38) 其中, $\omega$为$ {\sigma }_{\gamma }=0.29\times {10}^{-6} $的白噪声, $ T $为采样周期, 本文中仿真设定同步周期$T=1\;{\rm{s}}$.
主时钟考虑采用由卫星信号为基准的标准时钟信号, 即为理想时钟, 通过设定各个噪声和扰动的值, 对实际模型与主时钟进行对比, 得到开环时钟偏差信号, 如图4所示为自由时钟偏差$ \theta $的变化曲线. 在不加控制作用的情况下, 时钟偏差逐渐变大, 且自由时钟频率偏斜$ \chi \left(k\right) $如图5所示为一个正弦趋势变化的曲线.
利用参考时钟获取自由时钟偏差信号, 并进行快速傅里叶变换 (Fast Fourier transform, FFT), 得到时钟偏差信号的频谱图, 如图6所示. 从图6中可以看出, 其主频是频率$ f=0.1\;\mathrm{H}\mathrm{z} $的扰动信号且幅值约为$ 3.5\times {10}^{-6} $, 其他频点的信号较小, 主要是白噪声的影响, 与仿真设定噪声信号相符.
5.1 抗扰补偿器参数设计
基于扩展PI的抗扰补偿器的参数设计主要是对于4个参数矩阵${{\boldsymbol{K}}}_{{{1}}},{{\boldsymbol{K}}}_{{{2}}},{{\boldsymbol{K}}}_{{{3}}},{{\boldsymbol{K}}}_{{{4}}}$的参数选择, 通过扰动频率频谱分析得到其主频$ f=0.1\;\mathrm{H}\mathrm{z} $, 通过式$ \left(22\right) $即可求出${{\boldsymbol{K}}}_{{{1}}}$的值. 求解出${{\boldsymbol{K}}}_{{{1}}}$之后, 通过式(24)在系统稳定性和矩阵参数的约束条件下, 求得最优解${{\boldsymbol{K}}}_{{{2}}},{{\boldsymbol{K}}}_{{{3}}},{{\boldsymbol{K}}}_{{{4}}}$.
本文中, 对于特定频率$ f=0.1\;\mathrm{H}\mathrm{z} $的扰动信号, $ \omega =2\pi f, $通过式$ \left(22\right) $可知, 存在${{\boldsymbol{K}}}_{{{1}}}$为$2\times2 $的矩阵, 使得${{\boldsymbol{K}}}_{{{1}}}$矩阵的特征根为${\mathrm{e}}^{\pm 0.2\pi {\rm{j}}}$.
$$ {{\boldsymbol{K}}}_{{{1}}}=\left[ {\begin{array}{*{20}{c}} \mathrm{cos}\left(0.2\pi \right)& -\mathrm{sin}\left(0.2\pi \right)\\ \mathrm{sin}\left(0.2\pi \right)& \mathrm{cos}\left(0.2\pi \right) \end{array}} \right] $$ (39) 通过极点配置, 配置${{\boldsymbol{K}}}_{{{2}}},{{\boldsymbol{K}}}_{{{3}}},{{\boldsymbol{K}}}_{{{4}}}$矩阵的参数使得系统稳定, 最终求得扩展PI抗扰补偿器的极点, 分别是$0.8458\pm 0.5155{\rm{j}}$、$0.6891\pm 0.5874{\rm{j}}$, 都在单位圆内, 满足系统稳定条件.
为了对比性能, 搭建两个参考系统. 第1个参考系统采用卡尔曼滤波替换本文提出的扩展PI抗扰补偿器, PI控制器不变; 第2个参考系统仅包含PI控制器, 扰动补偿器部分采用传统扰动观测器−比例积分扰动观测器, 值得指出的是, 其为扩展PI抗扰补偿器的一种特例, 即当式(13)中${{\boldsymbol{K}}}_{{{1}}}= {{\boldsymbol{K}}}_{{{3}}}={\boldsymbol{I}}$为单位阵时, 扰动补偿器部分即为比例积分扰动观测器.
5.2 同步精度对比
同步误差的收敛曲线如图7所示. 从图7中可以看出, 传统扰动观测器的超调量和调节时间明显大于扩展PI抗扰补偿器和卡尔曼滤波的结果. 而添加卡尔曼滤波和扩展PI抗扰补偿器的超调量和调节时间等动态参数都比较接近, 因此从时间同步系统的动态跟踪特性来看, 扩展PI抗扰补偿器和卡尔曼滤波可以使得系统以较小的超调量更快地达到收敛.
稳态下的时间同步误差如图8所示, 工业物联网要求的高精度时间同步通常是指系统的稳态误差在要求的精度范围内, 即考虑的是系统的稳态指标. 图8所示的是3种时间同步方案下, 系统达到稳定后, 稳态误差波动范围.
从图8中可以看出, 传统扰动观测器和添加卡尔曼滤波后的时钟偏差都是正弦波形式振荡, 其主要是由于特定频率的扰动信号并没有被滤除, 其中传统扰动观测器相比于卡尔曼滤波周期性扰动信号表现不很明显的原因是卡尔曼滤波器对于白噪声的抑制作用较为明显, 凸显出较为规律的正弦扰动信号. 而扩展PI抗扰补偿器的时钟偏差波动曲线明显可以看出主要是随机扰动造成, 周期性扰动的影响已经基本被消除, 且其时间同步精度可以达到4 μs以内, 相比于传统方法精度提高了两倍.
5.3 同步误差频谱分析
为了更好地说明本文所提时间同步方法对于特定频谱扰动的抑制作用, 我们对时钟同步误差$ \theta $进行FFT, 得到时钟同步误差$ \theta $幅频特性对比图, 如图9所示.
从图9所示的同步误差的频谱分析可以看出, 传统扰动观测器和卡尔曼滤波后$ f=0.1\;\mathrm{H}\mathrm{z} $处的扰动仍然是主频, 说明其对特定频率的信号并没有明显的衰减作用, 只是对所有噪声有一定的衰减作用. 且从图9中可以看出, 卡尔曼滤波对于白噪声的抑制效果更好, 因此图8中卡尔曼滤波的稳态同步误差表现出较为规律的正弦扰动信号. 而针对本文考虑的低频噪声, 传统扰动观测器和卡尔曼滤波的效果不是很理想. 通过扩展PI抗扰补偿器输出的时钟偏差幅频图可以看出, 其对于特定频率的扰动具有很好的抑制作用, 这也进一步说明了图8中扩展PI抗扰补偿器稳态误差较小的原因.
综上分析, 扩展PI抗扰补偿器相比于卡尔曼滤波器和传统扰动观测器具有更好的抗干扰特性, 特别是对于特定频率的周期信号, 可以达到很好的抑制效果. 由此也进一步验证了本文提出的扩展PI抗扰补偿器具有良好的扰动抑制特性.
6. 结束语
考虑点对点时钟同步系统存在周期性的频率漂移、测量扰动和包交换随机延迟等情况, 本文构建了一种带扩展PI抗扰补偿器的时间同步系统, 通过零极点优化配置实现对特定频率的扰动信号的抑制作用, 克服了传统扰动观测器的零点固定的局限性, 并给出了外部控制器的参数取值范围. 实验仿真采用点对点网络进行分析, 通过与传统扰动观测器和卡尔曼滤波等传统方法的补偿效果对比, 凸显了扩展PI抗扰补偿器的优越性, 且本文提出的方案验证了最终达到时钟同步精度为4 μs, 实现了高精度时间同步.
-
-
[1] Zong Y, Dai X, Gao Z. Proportional-integral synchronisation for non-identical wireless packet-coupled oscillators with delays. IEEE Transactions on Industrial Electronics, 2021, 68(11): 11598-11608 doi: 10.1109/TIE.2020.3036228 [2] Elsts A, Fafoutis X, Duquennoy S, et al. Temperature-resilient time synchronization for the internet of things. IEEE Transactions on Industrial Informatics, 2018, 14(5): 2241-2250 doi: 10.1109/TII.2017.2778746 [3] Elson J. Fine-grained network time synchronization using reference broadcasts. ACM SIGOPS Operating Systems Review, 2003, 36(S1): 147-163 [4] Giorgi G, Narduzzi C. Performance analysis of kalman-filter-based clock synchronization in IEEE 1588 networks. IEEE Transactions on Instrumentation and Measurement, 2011, 60(8): 2902-2909 doi: 10.1109/TIM.2011.2113120 [5] Maróti M, Kusy B, Simon G, Lédeczi Á. The flooding time synchronization protocol. In: Proceedings of International Conference on Embedded Network Sensor System. Baltimore, USA: DBLP, 2004. [6] 孙子文, 吴梦芸, 白勇. 基于多参考节点的WSN时间同步方法. 控制工程, 2016, 23(01): 12-17Sun Zi-Wen, Wu Meng-Yun, Bai Yong. A time synchronization method for WSN based on multiple reference nodes, Control Engineering of China, 2016, 23(01): 12-17 [7] Wang H, Yu F, Li M, Zhong Y. Clock skew estimation for timestamp-free synchronization in industrial wireless sensor networks, IEEE Transactions on Industrial Informatics, 2021, 60(8): 2902–2909 [8] Yang T, Y Niu, Yu J. Clock synchronization in wireless sensor networks based on bayesian estimation. IEEE Access, 2020, 8: 69683-69694 doi: 10.1109/ACCESS.2020.2984785 [9] Carli R, Chiuso A, Schenato L, Zampieri S. Optimal synchronization for networks of noisy double integrators. IEEE Transactions on Automatic Control, 2011, 56(5): 1146-1152 doi: 10.1109/TAC.2011.2107051 [10] Yildirim K S. Gradient descent algorithm inspired adaptive time synchronization in wireless sensor networks. IEEE Sensors Journal, 2016, 16(13): 5463-5470 doi: 10.1109/JSEN.2016.2555996 [11] Yildirim K S, Carli R, Schenato L. Adaptive proportional-integral clock synchronization in wireless sensor networks. IEEE Transactions on Control Systems Technology, 2018, 26(2): 610-623 doi: 10.1109/TCST.2017.2692720 [12] Filler R L. The effect of vibration on quartz crystal resonators. Nasa Sti/Recon Technical Report N, 1980. [13] 杨东岳, 梅杰. 有向图中基于扰动观测器的线性多智能体系统一致性. 自动化学报, 2018, 44(6): 1037−1044Yang Dong-Yue, Mei Jie. Disturbance observer based consensus of linear multi-agent systems under a directed graph. Acta Automatica Sinica, 2018, 44(6): 1037−1044 [14] Zhou K, Zhang R. A new controller architecture for high performance, robust, and fault-tolerant control. IEEE Transactions on Automatic Control, 2001, 46(10): 1613-1618 doi: 10.1109/9.956059 [15] Wang D, Qiao J, Cheng L. An approximate neuro-optimal solution of discounted guaranteed cost control design. IEEE Transactions on Cybernetics, 2022, 52(1): 77−86 [16] F Ye, Sun B, Ou L, Zhang W. Disturbance observer-based control for consensus tracking of multi-agent systems with input delays from a frequency domain perspective. Systems & Control Letters, 2018, 114: 66-75 [17] Chen W H, Yang J, Guo L, Li S H. Disturbance-observer-based control and related methods—an overview. IEEE Transactions on Industrial Electronics, 2016, 63(2): 1083-1095 doi: 10.1109/TIE.2015.2478397 [18] 韩京清. 自抗扰控制器及其应用. 控制与决策, 1998, 13(1): 19–23 doi: 10.3321/j.issn:1001-0920.1998.01.005Han Jing-Qing. Auto-disturbances-rejection controller and It’s applications. Control and Decision, 1998, 13(1): 19–23 doi: 10.3321/j.issn:1001-0920.1998.01.005 [19] 王怡怡, 赵志良. 二自由度无人直升机的非线性自抗扰姿态控制. 自动化学报, 2021, 47(8): 1951-1962Wang Yi-Yi, Zhao Zhi-Liang. Nonlinear active disturbance rejection attitude control of two-DOF unmanned helicopter. Acta Automatica Sinica, 2021, 47(8): 1951-1962 [20] Gao Z. Scaling and bandwidth-parameterization based controller tuning. In: Proceedings of the American Control Conference. Denver, USA: IEEE, 2003. 4989–4996 [21] Dai X, Gao Z, Breikin T, Wang H. Zero assignment for robust H2/H∞ fault detection filter design. IEEE Transactions on Signal Processing, 2009, 57(4): 1363-1372 doi: 10.1109/TSP.2008.2010598 [22] Galleani, Lorenzo. A tutorial on the two-state model of the atomic clock noise. Metrologia, 2008, 45(6): S175-S182. doi: 10.1088/0026-1394/45/6/S23 [23] Kulpinski R J. Minimum variance methods for synchronization of airborne clocks. In: Proceedings of the 30th Annual Symposium on Frequency Control. IEEE, 1976. 401−413 [24] Lévesque M, Tipper D. A survey of clock synchronization over packet-switched networks. IEEE Communications Surveys & Tutorials, 2016, 18(4): 2926-2947 [25] Silvester J R. Determinants of Block Matrices. Mathematical Gazette, 2000, 84(501): 460-467 doi: 10.2307/3620776 [26] Schriegel S, Jasperneite J. Investigation of industrial environmental influences on clock sources and their effect on the synchronization accuracy of IEEE 1588. In: Proceedings of the IEEE International Symposium on Precision Clock Synchronization for Measurement. IEEE, 2007. 50−55 -