基于梯度方法的线性多延时系统参数辩识

Parameter Identification of Linear Multi-Delay System Based on Gradient Method

摘要: 针对一类线性多延时系统，提出一种基于梯度和最小二乘方法的参数辨识算法。该算法首先将系数矩阵和延时量的联合估计问题转化为仅关于延时量的优化问题，然后基于动量梯度法求出延时量的估计值，最后利用矩阵方程的最小二乘解得出系数矩阵的估计值，仿真结果表明，该算法有较好的辨识效果。

Abstract: For a class of linear systems with multiple delays, a new parameter identification algorithm based on gradient and least square method is proposed. Firstly, the joint estimation problem of coefficient matrix and time delay is transformed into the optimization problem which only estimates time delay. Then, the estimated value of time delay is obtained based on momentum gradient method. Finally, by using the least square solution of matrix equation, the estimated value of coefficient matrix is obtained. The algorithm has good identification effect.

1. 引言

随着社会的发展和科技的进步，延时系统在生物、经济和自动化控制等诸多领域日益受到重视 [1] [2]。参数辨识是设计、分析与控制实际系统的重要前提 [3]。

时域方法 [4] [5] [6] 是延时系统参数辨识领域的重要方法。基于梯度的辨识算法近年来受到许多学者的关注。文献 [7] [8] [9] 应用变分方法计算梯度并提出一种参数辨识算法；文献 [10] 求得了准则函数对参数的梯度所满足的系统，并由此给出了一种辨识算法。

本文受文献 [11] 启发，研究了一类线性多延时系统的系数和延时量的辨识问题。本文的主要工作是设计了关于系数矩阵和延时量的准则函数，利用最小二乘原理给出了固定延时量时系数矩阵的最小二乘解，从而将系数矩阵和延时量的同时估计问题看作仅关于延时量的估计问题，利用梯度方法求解优化问题得出延时量，最终用估计出的延时量求解矩阵方程求得系数矩阵。

2. 问题描述

考虑以下线性延时系统：

$\{\begin{array}{ll}\frac{\text{d}}{\text{d}t}x\left(t\right)={\displaystyle \underset{i=0}{\overset{N}{\sum }}}{A}_{i}x\left(t-{\tau }_i\right)+{\displaystyle \underset{i=0}{\overset{N}{\sum }}}{B}_{i}u\left(t-{\tau }_i\right),\hfill & t\in \left[0,T\right];\hfill \\ x\left(t\right)=\phi \left(t\right),\hfill & -\tau \le t\le 0.\hfill \end{array}$ (1)

其中 $x\left(t\right)\in {ℝ}^n$ 表示状态向量； $u\left(t\right)\in {ℝ}^p$ 表示控制向量； $A_i\in {ℝ}^{n\times n}$， $B_i\in {ℝ}^{n\times p},\left(i=0,1,2,\cdots ,N\right)$ 是系统的参数矩阵；初始函数 $\phi \left(t\right)\in {ℝ}^n$ ；这里假定系统的延时量满足 $0={\tau }_0<{\tau }_1<\cdots <{\tau }_N=\tau$。

记 $\Theta =\left[\begin{array}{cccc}{A}_0& A_1& \cdots & B_N\end{array}\right]$， $\bar{\tau }={\left[\begin{array}{ccc}{\tau }_1& \cdots & {\tau }_N\end{array}\right]}^{\text{T}}$。

定义2.1 设矩阵 $Q\in {ℝ}^{n\times m}$，称 ${\Vert Q\Vert }_F=\sqrt{trace\left(Q^{\text{T}}Q\right)}$ 为Q的Frobenius范数。

命题2.1 设 $Q=\left[\begin{array}{ccc}{q}_1& \cdots & q_m\end{array}\right]\in {ℝ}^{n\times m}$，则

$\begin{array}{c}{\Vert Q\Vert }_F=\sqrt{trace\left({\left[\begin{array}{ccc}{q}_1& \cdots & q_n\end{array}\right]}^{\text{T}}\left[\begin{array}{ccc}{q}_1& \cdots & q_n\end{array}\right]\right)}\\ =\sqrt{{\displaystyle \underset{i=1}{\overset{m}{\sum }}}{q}_i^{\text{T}}{q}_i}=\sqrt{{\displaystyle \underset{i=1}{\overset{m}{\sum }}}{\Vert q_i\Vert }_2^2}={\Vert {\left[\begin{array}{ccc}{q}_1^{\text{T}}& \cdots & q_m^{\text{T}}\end{array}\right]}^{\text{T}}\Vert }_2\end{array}$

引理2.1 [12] 若A是列满秩矩阵，B是行满秩矩阵，矩阵 $X\mathrm{,}F$ 满足

$AXB=F,$

那么在最小二乘意义下上式有唯一解 $X={\left(A^{\text{T}}A\right)}^{-1}{A}^{\text{T}}F{B}^{\text{T}}{\left(B{B}^{\text{T}}\right)}^{-1}$。

命题2.2 设矩阵函数 $Q\left(t\right)$ 可逆且可导，则 $\frac{\text{d}}{\text{d}t}{Q}^{-1}\left(t\right)=-Q^{-1}\left(t\right)\frac{\text{d}}{\text{d}t}Q\left(t\right)Q^{-1}\left(t\right)$。

命题2.3 设矩阵函数 $Q\left(t\right)\in {ℝ}^{n\times m}$ 可导，则

$\frac{\text{d}}{\text{d}t}{\Vert Q\left(t\right)\Vert }_F^2=\frac{\text{d}}{\text{d}t}tr\left(Q^{\text{T}}\left(t\right)Q\left(t\right)\right)=tr\left(\frac{\text{d}}{\text{d}t}\left(Q^{\text{T}}\left(t\right)Q\left(t\right)\right)\right)=2tr\left(Q^{\text{T}}\left(t\right)\frac{\text{d}}{\text{d}t}Q\left(t\right)\right).$

本文的目的是通过求解如下的优化问题辨识参数 $\Theta \mathrm{,}\bar{\tau }$，

$\begin{array}{c}\underset{\bar{\tau }\mathrm{,}\Theta }{min}J\left(\bar{\tau }\mathrm{,}\Theta \right)={\displaystyle \underset{l=1}{\overset{m}{\sum }}}{\Vert x\left(t_l\right)-x\left(0\right)-\left[{\displaystyle \underset{i=0}{\overset{N}{\sum }}}{A}_i{\displaystyle {\int }_0^{t_l}}x\left(t-{\tau }_i\right)\text{d}t+{\displaystyle \underset{i=0}{\overset{N}{\sum }}}{B}_i{\displaystyle {\int }_0^{t_l}}u\left(t-{\tau }_i\right)\text{d}t\right]\Vert }_2^2\\ ={\displaystyle \underset{l=1}{\overset{m}{\sum }}}{\Vert x\left(t_l\right)-x\left(0\right)-\Theta \varphi \left(\bar{\tau },l\right)\Vert }_2^2.\end{array}$ (2)

其中

$\varphi \left(\bar{\tau },l\right)={\left[\begin{array}{cccccc}{\displaystyle {\int }_0^{t_l}}{x}^{\text{T}}\left(t\right)\text{d}t& \cdots & {\displaystyle {\int }_0^{t_l}}{x}^{\text{T}}\left(t-{\tau }_N\right)\text{d}t& {\displaystyle {\int }_0^{t_l}}{u}^{\text{T}}\left(t\right)\text{d}t& \cdots & {\displaystyle {\int }_0^{t_l}}{u}^{\text{T}}\left(t-{\tau }_N\right)\text{d}t\end{array}\right]}^{\text{T}}.$ (3)

3. 算法的推导

固定延时量 $\bar{\tau }$，式(2)中 $\Theta$ 有唯一最小二乘解。从而对优化问题(2)的求解可转化为只与 $\bar{\tau }$ 有关的优化问题。接下来，利用新目标函数关于 $\bar{\tau }$ 的梯度采用动量梯度法估计延时量。最后，求系数矩阵关于式(2)的最小二乘解。

3.1. 梯度的计算

定理3.1 对于式(2)所述的优化问题，当 $\bar{\tau }$ 固定时有唯一最小二乘解

$\hat{\Theta }\left(\bar{\tau }\right)=\tilde{X}{\Phi }^{\text{T}}\left(\bar{\tau }\right){\left(\Phi \left(\bar{\tau }\right){\Phi }^{\text{T}}\left(\bar{\tau }\right)\right)}^{-1}.$ (4)

其中

$\tilde{X}=\left[\begin{array}{ccc}x\left(t_1\right)-x\left(0\right)& \cdots & x\left(t_m\right)-x\left(0\right)\end{array}\right],$ (5)

$\Phi \left(\bar{\tau }\right)=\left[\begin{array}{ccc}\varphi \left(\bar{\tau }\mathrm{,1}\right)& \cdots & \varphi \left(\bar{\tau }\mathrm{,}m\right)\end{array}\right]\mathrm{.}$ (6)

证明。

$J\left(\bar{\tau },\Theta \right)={\displaystyle \underset{l=1}{\overset{m}{\sum }}}{\Vert x\left(t_l\right)-x\left(0\right)-\Theta \varphi \left(\bar{\tau },l\right)\Vert }_2^2.$

由命题2.1，上式可写为 $J\left(\bar{\tau }\mathrm{,}\Theta \right)={\Vert \tilde{X}-\Theta \Phi \left(\bar{\tau }\right)\Vert }_F^2$。

由引理2.1上式中 $\Theta$ 有唯一最小二乘解。

${\Theta }^{*}=\hat{\Theta }\left(\bar{\tau }\right)=\tilde{X}{\Phi }^{\text{T}}\left(\bar{\tau }\right){\left(\Phi \left(\bar{\tau }\right){\Phi }^{\text{T}}\left(\bar{\tau }\right)\right)}^{-1}.$ (7)

$\square$

根据定理3.1可知，优化问题(2)可以转化为

$\underset{\bar{\tau }}{\min }{J}_{*}\left(\bar{\tau }\right)=\underset{\bar{\tau }}{\min }{\Vert \tilde{X}-\hat{\Theta }\left(\bar{\tau }\right)\Phi \left(\bar{\tau }\right)\Vert }_F^2.$ (8)

接下来分三个定理求 $\frac{\partial }{\partial \bar{\tau }}{J}_{\mathrm{<mo>\; *\; </mo>}}\left(\bar{\tau }\right)$。

定理3.2 设 $\varphi \left(\bar{\tau }\mathrm{,}l\right)\in {ℝ}^{\left(N+1\right)\left(n+p\right)}$，则

$\frac{\partial }{\partial {\tau }_j}\varphi \left(\bar{\tau },l\right)=\left[\begin{array}{c}O\\ -x\left(t_l-{\tau }_j\right)+x\left(-{\tau }_j\right)\\ O\\ -u\left(t_l-{\tau }_j\right)+u\left(-{\tau }_j\right)\\ O\end{array}\right].$

证明。

$\frac{\partial }{\partial {\tau }_j}\varphi \left(\bar{\tau },l\right)=\frac{\partial }{\partial {\tau }_j}\left[\begin{array}{c}{\displaystyle {\int }_0^{t_l}}x\left(t\right)\text{d}t\\ ⋮\\ {\displaystyle {\int }_0^{t_l}}x\left(t-{\tau }_N\right)\text{d}t\\ {\displaystyle {\int }_0^{t_l}}u\left(t\right)\text{d}t\\ ⋮\\ {\displaystyle {\int }_0^{t_l}}u\left(t-{\tau }_N\right)\text{d}t\end{array}\right]=\left[\begin{array}{c}\frac{\partial }{\partial {\tau }_j}{\displaystyle {\int }_0^{t_l}}x\left(t\right)\text{d}t\\ ⋮\\ \frac{\partial }{\partial {\tau }_j}{\displaystyle {\int }_0^{t_l}}x\left(t-{\tau }_N\right)\text{d}t\\ \frac{\partial }{\partial {\tau }_j}{\displaystyle {\int }_0^{t_l}}u\left(t\right)\text{d}t\\ ⋮\\ \frac{\partial }{\partial {\tau }_j}{\displaystyle {\int }_0^{t_l}}u\left(t-{\tau }_N\right)\text{d}t\end{array}\right]=\left[\begin{array}{c}O\\ -x\left(t_l-{\tau }_j\right)+x\left(-{\tau }_j\right)\\ O\\ -u\left(t_l-{\tau }_j\right)+u\left(-{\tau }_j\right)\\ O\end{array}\right].$

$\square$

定理3.3 设 $\hat{\Theta }\left(\bar{\tau }\right)\in {ℝ}^{n\times \left(N+1\right)\left(n+p\right)}$，则

$\begin{array}{c}\frac{\partial }{\partial {\tau }_j}\hat{\Theta }\left(\bar{\tau }\right)=\tilde{X}{\left(\frac{\partial }{\partial {\tau }_j}\Phi \left(\bar{\tau }\right)\right)}^{\text{T}}{\left(\Phi \left(\bar{\tau }\right){\Phi }^{\text{T}}\left(\bar{\tau }\right)\right)}^{-1}\\ -\tilde{X}{\Phi }^{\text{T}}\left(\bar{\tau }\right){\left(\Phi \left(\bar{\tau }\right){\Phi }^{\text{T}}\left(\bar{\tau }\right)\right)}^{-1}[\left(\frac{\partial }{\partial {\tau }_j}\Phi \left(\bar{\tau }\right)\right){\Phi }^{\text{T}}\left(\bar{\tau }\right)\\ +\Phi \left(\bar{\tau }\right)\frac{\partial }{\partial {\tau }_j}{\Phi }^{\text{T}}\left(\bar{\tau }\right)]{\left(\Phi \left(\bar{\tau }\right){\Phi }^{\text{T}}\left(\bar{\tau }\right)\right)}^{-1}.\end{array}$

其中 $\frac{\partial }{\partial {\tau }_j}\Phi \left(\bar{\tau }\right)=\left[\begin{array}{ccc}\frac{\partial }{\partial {\tau }_j}\varphi \left(\bar{\tau }\mathrm{,1}\right)& \cdots & \frac{\partial }{\partial {\tau }_j}\varphi \left(\bar{\tau }\mathrm{,}m\right)\end{array}\right]$。

证明。

$\begin{array}{l}\frac{\partial }{\partial {\tau }_j}\tilde{X}{\Phi }^{\text{T}}\left(\bar{\tau }\right){\left(\Phi \left(\bar{\tau }\right){\Phi }^{\text{T}}\left(\bar{\tau }\right)\right)}^{-1}=\tilde{X}\left(\frac{\partial }{\partial {\tau }_j}{\Phi }^{\text{T}}\left(\bar{\tau }\right)\right){\left(\Phi \left(\bar{\tau }\right){\Phi }^{\text{T}}\left(\bar{\tau }\right)\right)}^{-1}\\ +\tilde{X}{\Phi }^{\text{T}}\left(\bar{\tau }\right)\left(\frac{\partial }{\partial {\tau }_j}{\left(\Phi \left(\bar{\tau }\right){\Phi }^{\text{T}}\left(\bar{\tau }\right)\right)}^{-1}\right),\end{array}$

由命题2.2得

$\frac{\partial }{\partial {\tau }_j}{\left(\Phi \left(\bar{\tau }\right){\Phi }^{\text{T}}\left(\bar{\tau }\right)\right)}^{-1}=-{\left(\Phi \left(\bar{\tau }\right){\Phi }^{\text{T}}\left(\bar{\tau }\right)\right)}^{-1}\left[\left(\frac{\partial }{\partial {\tau }_j}\Phi \left(\bar{\tau }\right)\right){\Phi }^{\text{T}}\left(\bar{\tau }\right)+\Phi \left(\bar{\tau }\right)\frac{\partial }{\partial {\tau }_j}{\Phi }^{\text{T}}\left(\bar{\tau }\right)\right]{\left(\Phi \left(\bar{\tau }\right){\Phi }^{\text{T}}\left(\bar{\tau }\right)\right)}^{-1}.$

$\square$

定理3.4 式(8)对延时量 ${\tau }_j$ 的偏导数如下

$\frac{\partial }{\partial {\tau }_j}{\Vert \tilde{X}-\hat{\Theta }\left(\bar{\tau }\right)\Phi \left(\bar{\tau }\right)\Vert }_F^2=-2tr\left({\left(\tilde{X}-\hat{\Theta }\left(\bar{\tau }\right)\Phi \left(\bar{\tau }\right)\right)}^{\text{T}}\left(\left(\frac{\partial }{\partial {\tau }_j}\hat{\Theta }\left(\bar{\tau }\right)\right)\Phi \left(\bar{\tau }\right)+\hat{\Theta }\left(\bar{\tau }\right)\left(\frac{\partial }{\partial {\tau }_j}\Phi \left(\bar{\tau }\right)\right)\right)\right)\mathrm{.}$ (9)

证明。

$\begin{array}{l}\frac{\partial }{\partial {\tau }_j}{\Vert \tilde{X}-\hat{\Theta }\left(\bar{\tau }\right)\Phi \left(\bar{\tau }\right)\Vert }_F^2\\ =2tr\left({\left(\tilde{X}-\hat{\Theta }\left(\bar{\tau }\right)\Phi \left(\bar{\tau }\right)\right)}^{\text{T}}\frac{\partial }{\partial {\tau }_j}\left(\tilde{X}-\hat{\Theta }\left(\bar{\tau }\right)\Phi \left(\bar{\tau }\right)\right)\right)\\ =-2tr\left({\left(\tilde{X}-\hat{\Theta }\left(\bar{\tau }\right)\Phi \left(\bar{\tau }\right)\right)}^{\text{T}}\left(\left(\frac{\partial }{\partial {\tau }_j}\hat{\Theta }\left(\bar{\tau }\right)\right)\Phi \left(\bar{\tau }\right)+\hat{\Theta }\left(\bar{\tau }\right)\left(\frac{\partial }{\partial {\tau }_j}\Phi \left(\bar{\tau }\right)\right)\right)\right)\mathrm{.}\end{array}$

$\square$

3.2. 算法描述

本文采用动量梯度法进行联合估计(见算法1)。

4. 数值算例

本节给出两个数值算例来验证结果的有效性

例4.1 考虑不稳定的线性单延时系统(1)，系统参数如下

$A_0=\left[\begin{array}{cc}-0.2511& -0.1511\\ -0.3194& -0.2511\end{array}\right],A_1=\left[\begin{array}{cc}0.3682& 0.3682\\ 0.1000& 0.3682\end{array}\right],$

$B_0=\left[\begin{array}{c}2.2000\\ 0.2000\end{array}\right],B_1=\left[\begin{array}{c}3.1000\\ 1.1000\end{array}\right],\bar{\tau }=\left[\begin{array}{c}0.8\end{array}\right],$

输入函数为 $u\left(t\right)=3\sin \left(3t\right)-\cos \left(t\right)$，初始函数为 $\phi \left(t\right)={\left[\begin{array}{cc}0& 0\end{array}\right]}^{\text{T}},T=5$。

初始点 ${\bar{\tau }}_0=\begin{array}{c}0.7\end{array}$，参数 $\alpha =0.7,\beta =0.2,t_l=1+0.02\left(l-1\right),l=1,\cdots ,201$。

利用算法1出参数

$A_0=\left[\begin{array}{cc}-0.2511& -0.1505\\ -0.3195& -0.2511\end{array}\right],A_1=\left[\begin{array}{cc}0.3684& 0.3672\\ 0.1001& 0.3682\end{array}\right],$

$B_0=\left[\begin{array}{c}2.2004\\ 0.2000\end{array}\right],B_1=\left[\begin{array}{c}3.0997\\ 1.1000\end{array}\right],\bar{\tau }=\begin{array}{c}0.8005\end{array}.$

用算法1估计例4.1中延时量的迭代过程如图1所示。

Figure 1. Iterative process of delay estimation

图1. 延时量估计的迭代过程

例4.2 考虑线性多延时系统(1)，系统参数如下

$A_0=\left[\begin{array}{cc}-15.5000& -15.5000\\ -15.0000& -15.5000\end{array}\right],A_1=\left[\begin{array}{cc}-10.4110& -10.4110\\ -1.0000& -10.4110\end{array}\right],A_2=\left[\begin{array}{cc}-5.9110& -5.9110\\ -6.5000& -5.9110\end{array}\right],$

$\begin{array}{l}{B}_0=\left[\begin{array}{c}0.5000\\ 1.0000\end{array}\right],B_1=\left[\begin{array}{c}3.1000\\ 1.1000\end{array}\right],\\ B_2=\left[\begin{array}{c}2.2000\\ 0.2000\end{array}\right],\bar{\tau }={\left[\begin{array}{cc}0.2& 0.8\end{array}\right]}^{\text{T}},\end{array}$

输入函数为 $u\left(t\right)=10+\sin \left(3t\right)-5\cos \left(t\right)$，初始函数为 $\phi \left(t\right)={\left[\begin{array}{cc}0& 0\end{array}\right]}^{\text{T}},T=10$ .

初始点 ${\bar{\tau }}_0={\left[\begin{array}{cc}0.25& 0.85\end{array}\right]}^{\text{T}}$，参数 $\alpha =0.9,\beta =20,t_l=1+0.5\left(l-1\right),l=1,\cdots ,19$ .

利用算法1出参数

$A_0=\left[\begin{array}{cc}-15.4998& -15.4998\\ -14.9999& -15.4999\end{array}\right],A_1=\left[\begin{array}{cc}-10.4109& -10.4109\\ -1.0000& -10.4109\end{array}\right],A_2=\left[\begin{array}{cc}-5.9109& -5.9109\\ -5.9109& -5.9109\end{array}\right],$

$\begin{array}{l}{B}_0=\left[\begin{array}{c}0.5000\\ 1.0000\end{array}\right],B_1=\left[\begin{array}{c}3.1000\\ 1.1000\end{array}\right],\\ B_2=\left[\begin{array}{c}2.2000\\ 0.2000\end{array}\right],\bar{\tau }=\left[\begin{array}{c}0.2001\\ 0.7998\end{array}\right].\end{array}$

用算法1估计例4.2中两个延时量的迭代过程通过等高线图方式呈现如图2所示。将图2中目标函数和迭代过程用三维曲面呈现后如图3所示。

Figure 2. Iterative process of two delay estimation

图2. 两个延时量估计的迭代过程

Figure 3. The iterative process and objective function of the delay amount

图3. 延时量的迭代过程及目标函数

5. 结论

针对线性多延时系统，给出了一种含积分的目标函数。通过最小二乘方法将联合参数的目标函数转化为仅依赖延时量的非线性函数。采用动量梯度方法求解关于延时量的优化问题。确定延时量后求矩阵方程的最小二乘解估计系数矩阵。通过数值仿真验证了算法的有效性。