-
非线性动力学的稀疏辨识(Spark Identification of Nonlinear Dynamics, SINDY)[12, 14]是由Brunton等开发的利用稀疏识别的方式,从含有噪声的观测数据中识别动力学系统的控制方程。
考虑一个自治系统:
$$ \dot{\boldsymbol{x}}\left(t\right)=\boldsymbol{f}\left(\boldsymbol{x}\left(t\right)\right) $$ (1) $$ \boldsymbol{y}\left(t\right)=\boldsymbol{g}\left(\boldsymbol{x}\right(t\left)\right) $$ (2) 式中:
$\boldsymbol{x}\left(t\right)={\left[\begin{array}{c}{x}_{1}\left(t\right)\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}{x}_{2}\left(t\right)\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\cdots \hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}{x}_{n}\left(t\right)\end{array}\right]}^{{\rm{T}}}\in {\mathbb{{{R}}}}^{n}$ 代表系统在时刻t的状态;函数$ \boldsymbol{f}\left(\boldsymbol{x}\right(t\left)\right) $ 表示系统的运动方程;函数$ \boldsymbol{g}\left(\boldsymbol{x}\right(t\left)\right) $ 表示系统的输出方程。为了从带有噪声的测量数据中获得$ \boldsymbol{f}\left(\boldsymbol{x}\right(t\left)\right) $ 函数,通过测量状态及状态的导数形成两个状态及状态导数矩阵。$$ \boldsymbol{X}=\left[\begin{array}{c}{\boldsymbol{x}}^{{\rm{T}}}\left({t}_{1}\right)\\ {\boldsymbol{x}}^{{\rm{T}}}\left({t}_{2}\right)\\ \vdots \\ {\boldsymbol{x}}^{{\rm{T}}}\left({t}_{m}\right)\end{array}\right]=\left[\begin{array}{cccc}{x}_{1}\left({t}_{1}\right)& {x}_{2}\left({t}_{1}\right)& \cdots & {x}_{n}\left({t}_{1}\right)\\ {x}_{1}\left({t}_{2}\right)& {x}_{2}\left({t}_{2}\right)& \cdots & {x}_{n}\left({t}_{2}\right)\\ \vdots & \vdots & \ddots & \vdots \\ {x}_{1}\left({t}_{m}\right)& {x}_{2}\left({t}_{m}\right)& \cdots & {x}_{n}\left({t}_{m}\right)\end{array}\right] $$ (3) $$ \dot{\boldsymbol{X}}=\left[\begin{array}{c}{\dot{\boldsymbol{x}}}^{{\rm{T}}}\left({t}_{1}\right)\\ {\dot{\boldsymbol{x}}}^{{\rm{T}}}\left({t}_{2}\right)\\ \vdots \\ {\dot{\boldsymbol{x}}}^{{\rm{T}}}\left({t}_{m}\right)\end{array}\right]=\left[\begin{array}{cccc}{\dot{x}}_{1}\left({t}_{1}\right)& {\dot{x}}_{2}\left({t}_{1}\right)& \cdots & {\dot{x}}_{n}\left({t}_{1}\right)\\ {\dot{x}}_{1}\left({t}_{2}\right)& {\dot{x}}_{2}\left({t}_{2}\right)& \cdots & {\dot{x}}_{n}\left({t}_{2}\right)\\ \vdots & \vdots & \ddots & \vdots \\ {\dot{x}}_{1}\left({t}_{m}\right)& {\dot{x}}_{2}\left({t}_{m}\right)& \cdots & {\dot{x}}_{n}\left({t}_{m}\right)\end{array}\right] $$ (4) 然后构建一个函数库,该函数库包含了一些由线性及非线性函数组成的基函数:
$$ {\varTheta }\left(\boldsymbol{X}\right)=\left[\begin{array}{cccccccccc}& & & & & & & & & \\ 1& \boldsymbol{X}& {\boldsymbol{X}}^{{P}_{2}}& {\boldsymbol{X}}^{{P}_{3}}& \cdots & \mathrm{s}\mathrm{i}\mathrm{n}\left(\boldsymbol{X}\right)& \mathrm{c}\mathrm{o}\mathrm{s}\left(\boldsymbol{X}\right)& \mathrm{s}\mathrm{i}\mathrm{n}\left(2\boldsymbol{X}\right)& \mathrm{c}\mathrm{o}\mathrm{s}\left(2\boldsymbol{X}\right)& \cdots \\ & & & & & & & & & \end{array}\right] $$ (5) 式中:
$ {\boldsymbol{X}}^{{P}_{2}} $ 、$ {\boldsymbol{X}}^{{P}_{3}} $ 表示高阶的多项式项;$ {P}_{2} $ 表示二次多项式项。$ {\varTheta }\left(\boldsymbol{X}\right) $ 的每一列代表一个候选基函数,基函数可以自由选择,取决于对系统的具体分析。对于多数动力学系统只包括基函数中的几项[15]。因此,可以使用稀疏回归的方法来求得模型的系数,该方法平衡了系数的稀疏性和模型的准确性。设模型系数为$ {\varXi }=\left[\begin{array}{c}{\xi }_{1}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}{\xi }_{2}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\cdots \hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}{\boldsymbol{\xi }}_{n}\end{array}\right] $ ,$ {\xi }_{k} $ 表示基函数库$ {\varTheta }\left(\boldsymbol{X}\right) $ 中第k个状态函数的系数。非线性动力学系统可以表示为:$$ \widehat {\boldsymbol{X}} = \varTheta \left( {\boldsymbol{X}} \right)\varXi $$ (6) 模型精度的成本函数定义如下:
$$ {\xi }_{k}={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{{\xi }_{k}}{\left\| {\dot{{\boldsymbol{X}}}}_{k}-{\xi }_{k}{{\varTheta }}^{{\rm{T}}}\left({\boldsymbol{X}}\right) \right\|}_{2}+\lambda {\left\| {\xi }_{k} \right\|}_{1} $$ (7) 式中:
${\dot{\boldsymbol{X}}}_{k}$ 表示$ \dot{\boldsymbol{X}} $ 的第k行;$ \lambda $ 为惩罚因子,使用该系数来平衡模型复杂度与精度之间的关系。为了获得准确的模型系数,文中采用序列阈值最小二乘法[14]来辨识模型的参数,该算法对噪声具有非常强的鲁棒性,甚至当测量数据与噪声数据十分接近时,序列阈值最小二乘法也能表现出非常好的精度,并且可以快速收敛。序列阈值最小二乘法的计算过程是从
$ {\varXi } $ 的最小二乘解开始,然后对所有小于$ \lambda $ 值的系数进行阈值化。阈值化后,非零系数再进行最小二乘求解,继续以上执行过程,直到非零系数收敛为止,流程图如图1所示。为了确定系统输入对系统的影响,建立精确的动力系统输入输出模型,Brunton Steven L等 [12]将SINDY推广到SINDYc。SINDYc算法与SINDY算法思路相同,区别在于SINDYc算法扩大了基函数库,包括了输入U与状态X的交叉项。系统方程变为:
$$ \dot{\boldsymbol{X}}={\varXi }{{\varTheta }}^{{{\rm{T}}}}({{\boldsymbol{X}}},{{\boldsymbol{U}}}) $$ (8) 与SINDY使用相同的方法可以获得系统参数
$ \mathrm{\Xi } $ ,从而获得精确的系统动力学方程。 -
由于实际的望远镜伺服系统模型十分复杂,包括各种非线性因素,为了简化分析,文中建立了简化的望远镜伺服系统的二质量动力学模型。由于望远镜在跟踪过程中运动速度很慢,方位轴系运动与高度轴系运动的耦合作用很小,因此可以将方位轴系运动与高度轴系运动独立分析。以方位轴系为例,图2为望远镜方位轴系的简化二质量动力学模型。
方位轴系的运动可以由以下两个方程描述:
$$ \begin{array}{c}{T}_{m}={J}_{m}{s}^{2}{\theta }_{m}+{B}_{m}s{\theta }_{m}+{B}_{L}\left(s{\theta }_{m}-s{\theta }_{L}\right)+{K}_{L}\left({\theta }_{m}-{\theta }_{L}\right)\end{array} $$ (9) $$ {T}_{L}={J}_{L}{s}^{2}{\theta }_{L}+{B}_{L}\left(s{\theta }_{L}-s{\theta }_{m}\right)+{K}_{L}\left({\theta }_{L}-{\theta }_{m}\right) $$ (10) 式中:
$ {T}_{m} $ 为电机输出力矩;$ {J}_{m} $ 为方位轴的转动惯量;$ {B}_{m} $ 为方位轴的阻尼;$ {B}_{L} $ 为方位轴与负载之间的阻尼;$ {T}_{L} $ 为负载受到的外部扰动力矩;$ {J}_{L} $ 为负载转动惯量;$ {K}_{L} $ 为方位轴与负载之间的扭转刚度;$ {\theta }_{m} $ 为方位轴转动角度;$ {\theta }_{L} $ 为负载转动角度。函数库
$\mathrm{\varTheta }(X,U)$ 中的函数形式有很多选择,文中择线性函数作为模型的候选函数。选取四个状态变量$ \left[\begin{array}{c}{X}_{0}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}{X}_{1}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}{X}_{2}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}{X}_{3}\end{array}\right] $ ,分别为:方位轴电机转动角度$ {\theta }_{m} $ (″);方位轴电机转速$ {\dot{\theta }}_{m} $ ((″)/s);方位轴负载转动角度$ {\theta }_{L} $ (″);方位轴负载转速$ {\dot{\theta }}_{L} $ ((″)/s)。写成如下形式:$$ {\boldsymbol{X}}=\left[\begin{array}{c}{\theta }_{m}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}{\dot{\theta }}_{m}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}{\theta }_{L}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}{\dot{\theta }}_{L}\end{array}\right] $$ (11) $$ \dot{\boldsymbol{X}}=\left[\begin{array}{c}{\dot{\mathrm{\theta }}}_{\mathrm{B}}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}{\ddot{\mathrm{\theta }}}_{\mathrm{B}}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}{\dot{\mathrm{\theta }}}_{\mathrm{L}}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}{\ddot{\mathrm{\theta }}}_{\mathrm{L}}\end{array}\right] $$ (12) 首先给定位置、速度的初始值X0=[−2 0 2 1],并通过计算获得模拟数据。为了验证该算法对具有随机误差数据的模型预测的有效性,在模拟数据中增加了高斯伪随机噪声扰动信号,随机噪声扰动信号的幅值分别设定为信号幅值的0%、5%、10%、20%、30%及50%。图3为状态变量变化曲线,其中红色虚线为模型预测曲线,蓝色实线为数值模拟曲线。从图中可以看出,红色虚线的预测曲线在不同的噪声水平下都有很好的跟随效果,在10%噪声水平以上时,红色虚线和蓝色实线曲线出现了一定的差异,在10%噪声水平以下时,红、蓝曲线的最大偏差值全部在信号幅值的5%以内。表1为各个状态跟随偏差。
表 1 状态随偏差的变化
Table 1. State change deviation
Noise level 0 5% 10% 20% 30% 50% X0 0.130% 4.908% 2.635% 4.600% 9.794% 14.263% X1 0.060% 1.112% 2.642% 7.372% 14.341% 31.482% X2 4.120% 3.524% 4.192% 3.343% 2.169% 0.937% X3 0.004% 0.485% 0.109% 4.122% 10.839% 28.409% 辨识方法获得的模型是真实系统的近似描述,如果辨识出的模型与实际系统在相同输入的前提下,系统状态(输出)曲线基本一致,则可以认为辨识模型可以准确地描述真实的系统。为了定量地衡量辨识模型与真实系统的相似程度,需要建立一个评判标准,文中采用决定系数
$ {R}^{2} $ 衡量辨识模型的准确度。$ {R}^{2} $ 定义如下:$$ {R}^{2}\equiv 1-\frac{S{S}_{\mathrm{r}\mathrm{e}\mathrm{s}}}{S{S}_{\mathrm{t}\mathrm{o}\mathrm{t}}} $$ (13) $$S{S_{{\rm{res}}}} = \sum\limits_i {{{({y_i} - {f_i})}^2} = \sum\limits_i {e{i^2}} } $$ (14) $$ S{S_{{\rm{tot}}}} = \sum\limits_i {{{({y_i} - \overline y )}^2}} $$ (15) 式中:
${SS}_{{\rm{res}}}$ 为数据的残差平方和;${{SS}}_{{\rm{tot}}}$ 为数据的总平方和;$ {f}_{i} $ 为回归值;$ {y}_{i} $ 为测量值;$ \overline{y} $ 为平均值。通过上述模拟结果拟合
$ {R}^{2} $ 数据曲线如图4所示,横坐标为噪声水平,纵坐标为决定系数$ {R}^{2} $ 。可以看出噪声水平在20%以内时,决定系数在0.99以上,SINDYc算法拟合的模型可以准确地模拟实际系统的输出。 -
为了验证SINDYc算法,在南极望远镜实验平台上进行了辨识实验,实验设备包括Power UMAC主控制器、ETEL力矩电机、海德汉光栅鼓等,望远镜方位轴系的控制方式采用直接PWM控制,控制结构采用典型的位置、速度和电流环三环结构,三环的控制算法在主控制器中实现。系统辨识实验原理如图5所示,速度环与位置环设置为开路。在实验过程中,是否能对系统施加合适的激励是能否成功辨识的一个关键因素。实验设计两种类型的输入信号以对比两种方案。方案一输入信号采用方波信号,如图6所示,幅值为最大输出量的10%;方案二输入信号采用5 Hz正弦信号,如图7所示,同时以2257 Hz的采样率对编码器输出信号进行采样。
使用SINDYc算法的辨识结果如图8和图9所示,辨识精度指标分别为R2=0.9857和R2=0.9952。从图中可以看出,在两类数据信号的激励下,SINDYc算法都能很好地识别出系统的动态跟随特性,而采用正弦输入信号的激励下,辨识出的模型精度比方波信号激励下的模型精度高0.95%。
Model prediction of optical infrared telescope servo system in extreme environment
-
摘要: 在极端环境下光学红外望远镜伺服系统建模的实际工作过程中,实测的望远镜状态数据经常包含有各种噪声。为了减小噪声对模型辨识精度的影响,提出了一种基于带有控制的非线性动力学稀疏辨识(Spark Identification of Nonlinear Dynamics with Control, SINDYc)算法的稀疏辨识方法。针对望远镜伺服系统,对SINDYc算法进行了理论分析和数值模拟,对比了在不同的噪声水平下,望远镜伺服系统预测模型的状态变量曲线,并拟合了不同噪声水平下辨识模型的决定系数曲线。基于南极望远镜实验平台,设计了正弦和方波信号作为激励信号进行模型辨识实验,对SINDYc算法的建模准确性进行了实验验证。数值模拟模型的预测输出结果显示:SINDYc算法在20%噪声水平以下时,模型辨识精度在0.99以上;在10%噪声水平以下时,状态变化跟随最大偏差值在信号幅值的5%以内。辨识实验数据表明,在两种不同信号激励下望远镜伺服系统模型预测的辨识精度分别为0.9857与0.9952,证实了基于SINDYc算法的稀疏辨识方法的有效性和准确性。该方法辨识出的系统模型可以为未来的南极大口径光学红外望远镜控制系统的分析及控制器设计提供很好的分析模型。Abstract: In practical terms of model identification for optical infrared telescope servo system, the measured states often contain noise. In order to improve accuracy of modeling, a sparse identification method based on Spark Identification of Nonlinear Dynamics with Control (SINDYc) algorithm was proposed. For the telescope servo system, theoretical analysis and numerical simulation of the SINDYc algorithm were carried out. The state variable curves of the telescope servo system model under different noise levels were compared, and the determination coefficient curve of the identification model under different noise levels was fitted. Based on the antarctic telescope experimental platform, sine and square wave signals were designed as excitation signals for model identification experiments, and the modeling accuracy of the SINDYc algorithm was experimentally verified. The prediction results of the numerical simulation model show that the identification accuracy of SINDYc algorithm is above 0.99 when the noise level is below 20%. When the noise level is below 10%, the maximum deviation is within 5% of the signal amplitude. The identification experimental data show that the accuracy of model prediction is 0.9857 and 0.9952 under the excitation of two different signals. The effectiveness and accuracy of the sparse identification method based on the SINDYc algorithm are confirmed by numerical simulation and experimental validation. The obtained model by this identification method can provide a good analytical model for the analysis and controller design of future large-aperture optical infrared telescope control systems.
-
表 1 状态随偏差的变化
Table 1. State change deviation
Noise level 0 5% 10% 20% 30% 50% X0 0.130% 4.908% 2.635% 4.600% 9.794% 14.263% X1 0.060% 1.112% 2.642% 7.372% 14.341% 31.482% X2 4.120% 3.524% 4.192% 3.343% 2.169% 0.937% X3 0.004% 0.485% 0.109% 4.122% 10.839% 28.409% -
[1] Sun Haotian, Du Fujia, Zhang Zhiyong. Large aperture fast steering mirror servo system based on model predictive control [J]. Infrared and Laser Engineering, 2020, 49(2): 0214001. (in Chinese) [2] Fang Lianwei, Shi Shouxia, Jiang Zhiyong. Servo mechanism parameter identification of fast steering mirror based on flexible supports [J]. Infrared and Laser Engineering, 2020, 50(5): 20200303. (in Chinese) [3] Wang Hao, Liu Jinghong, Deng Yongting, et al. Control model identification of opto-electronic tracking turntable [J]. Infrared and Laser Engineering, 2016, 45(6): 0617007. (in Chinese) [4] Xia Peipei, Deng Yongting, Wang Zhiqian, et al. Model identification for K mirror turntable of 2 m telescope [J]. Infrared and Laser Engineering, 2018, 47(3): 0318001. (in Chinese) [5] Zhou Wangping, Wang Ya, Gong Liyuan. Modeling and decoupling of rack motion of slant mount telescopes [J]. Computer Simulation, 2017, 34(11): 248-252. (in Chinese) doi: 10.3969/j.issn.1006-9348.2017.11.054 [6] Deng Yongting, Li Hongwen, Chen Tao. Dynamic analysis of two meters telescope mount control system [J]. Optics and Precision Engineering, 2018, 26(3): 654-661. (in Chinese) [7] Deng Yongting, Liu Jun, Li Hongwen, et al. Control system of 4 meters telescope based on segmented permanent magnet arc synchronous motor [J]. Optics and Precision Engineering, 2020, 28(3): 591-600. (in Chinese) [8] Dong Quanrui, Chen Tao, Gao Shijie, et al. Identification of opto-electronic fine tracking systems based on an improved differential evolution algorithm [J]. Chinese Optics, 2020, 13(6): 1314-1323. (in Chinese) doi: 10.37188/CO.2020-0021 [9] Jiménez-Garcı́a S, Mario M E, Benı́tez-Read J S, et al. Modelling, simulation, and gain scheduling control of large radio telescopes [J]. Simulation Practice & Theory, 2000, 8(3-4): 141-160. [10] Enrico Cascone, Dario Mancini, Pietro Schipani. Galileo telescope model identification [J]. Telescope Control Systems II, 1997, 3112: 343-350. doi: 10.1117/12.284232 [11] Gawronski W, Bartos R. Modeling and analysis of the DSS-14 antenna control system [J]. Telecommunications & Data Acquisition Progress Report, 1996, 124: 113-142. [12] Brunton S L, Proctor J L, Kutz J N. Sparse identification of nonlinear dynamics with control (SINDYc) [J]. IFAC-Papers OnLine, 2016, 49(18): 710-715. doi: 10.1016/j.ifacol.2016.10.249 [13] Zhang L, Schaeffer H. On the convergence of the SINDy algorithm [J]. Multiscale Modeling & Simulation, 2019, 17(3): 948-972. [14] Brunton S L, Proctor J L, Kutz J N. Discovering governing equations from data by sparse identification of nonlinear dynamical systems [J]. Proceedings of the National Academy of Sciences, 2016, 113(15): 3932-3937. doi: 10.1073/pnas.1517384113 [15] Sun Chong, Tian Tian, Zhu Xiaocheng, et al. Sparse identification of nonlinear unsteady aerodynamics of the oscillating airfoil [J]. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 2021, 235(7): 809-824.