自从Gauss在1775年提出LS算法以来, 对其进行改进和完善的报导几乎没有停止过^[1,2], 但在实际工程中大量使用的仍是普通LS算法。这主要是因为它计算简单, 便于理解和掌握。当然LS算法确实存在诸如不稳健、易出现病态估计等缺陷, 许多文献给出了值得称道的改进算法, 但多数是停留在理论或仿真研究上, 很少有应用于工程实践的。其主要原因是改进算法复杂难以理解, 且研究人员受实践条件的制约。另外, 这些算法成立的前提是系统要满足持续激励条件, 持续激励条件不仅影响系统的可估计性, 而且影响自适应控制系统的鲁棒性。而实际过程中许多系统的输入信号变化不大, 持续激励条件难以满足, 普通LS算法及其相关改进算法不能直接使用, 必须外加伪随机信号, 这将给工程实现带来许多麻烦。事实上, 许多工业过程的先验知识是可以利用的, 根据先验知识重新修正目标函数是克服LS算法病态估计的有效途径之一。

笔者立足于工程实践, 提出了基于修正目标函数的时变参数估计算法, 有效地克服了病态估计, 巧妙地避免使用伪随机信号。同时考虑到某些生产过程参数的高度时变性, 在新算法中引入加权矩阵, 强化当前观测数据对参数估计的作用, 提高跟踪时变参数的能力。新算法一方面摆脱了伪随机信号带来的麻烦;另一方面避免了数据饱和现象, 提高了自适应控制系统的鲁棒性和控制精度。它已用于实际工业生产, 取得了良好效果。仿真结果及实际应用结果都证明:新型估计算法简单易行、收敛速度快和估计精度高。

时变参数估计是工程实践中经常遇到的问题, 为了使所估计的参数能跟踪实际参数的变化, 需要人为地削弱旧数据的影响, 以增强新数据的权。常用的办法主要有两类:a. 采用矩形窗函数, 舍弃旧的数据, 只用一段新的数据进行参数估计的限定记忆法;b. 是通过遗忘因子对历史数据进行指数加权, 以逐渐削除旧数据影响的渐消记法, 本文的方法只涉及第二类。

考虑单变量离散线性时变系统

$y{}_{(t)}=\mathit{\lambda }{}_{(t)}^{\text{Τ}}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}+\tilde{\mathit{\omega }}{}_{(t)}‚(1)$

式中t=1, 2, 3, …,

$\begin{array}{l}\mathit{\lambda }{}_{(t)}^{\text{Τ}}=[-y(t-1),-y(t-2),\cdots ,-y(t-n{}_{\text{a}}),\\ u(t-1),u(t-2),\cdots ,u(t-n{}_{\text{b}})]‚(1\text{a})\\ \mathit{\Xi }{}_{(t)}^{\text{Τ}}=[a{}_{1(t)},a{}_{2(t)},\cdots ,a{}_{\text{n}\text{a}(t)},\\ b{}_{1(t)},b{}_{2(t)},\cdots ,b{}_{\text{n}\text{b}(t)}],(1\text{b})\end{array}$

其中{u _(i) }, {y _(j) } 分别是系统的输入和输出测量序列, $\tilde{\omega }$ _(t) 为随机噪声, n_a, n_b为过程的结构参数, a_{i (t)} , b_{j (t)} 是待识时变参数。

为跟踪参数Ξ _(t) 的变化应引入遗忘因子。

定义1 设遗忘因子为δ, (0<δ<1) 且定义加权矩阵为

$\mathit{W}{}_{(t)}=\left[\begin{array}{ccccc}\delta {}^{2(t-1)}& & & & \\ & \delta {}^{2(t-2)}& & 0& \\ & & \ddots & & \\ & 0& & \delta {}^2& \\ & & & & 1\end{array}\right]‚(2)$

式中t为总测量次数。在LS算法目标函数中引入W _(t) 是准确跟踪时变参数实现实时估计的重要途径。

某些工业过程的先验知识可用来克服输入信号不满足持续激励条件引起的病态估计。这里假设先验知识向量为

$\mathit{{\rm M}}=[\overline{a}{}_1‚\overline{a}{}_2‚\cdots ,\overline{a}{}_{\text{n}\text{a}}‚\overline{b}{}_1‚\overline{b}{}_2‚\cdots ,\overline{b}{}_{\text{n}\text{b}}]‚(3)$

式中{$\overline{a}$_i}, {$\overline{b}$_j}是待估计过程的时变参数均值, 参数均值这一先验知识可以通过日常生产数据分析获得, 这个先验知识是指导选样的重要参考。

为克服病态估计, 在新算法中应尽量利用先验知识向量M, 以降低过程对持续激励条件的要求;为适应时变过程的要求, 新算法还必须重视W _(t) 的作用。考虑到上述两个要求, 根据最小二乘原理可定义下列二次型目标函数。

定义2 基于修正目标函数的时变参数LS估计算法的目标函数为

$\begin{array}{l}J=[\mathit{Y}{}_{(t)}-\mathit{\Lambda }{}_{(t)}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}]{}^{\text{Τ}}\mathit{W}{}_{(t)}[\mathit{Y}{}_{(t)}-\mathit{\Lambda }{}_{(t)}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}]+\\ \alpha [\mathit{{\rm M}}{}_{(t)}-\frac{1}{t}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}]{}^{\text{Τ}}[\mathit{{\rm M}}{}_{(t)}-\frac{1}{t}[\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}]‚(4)\end{array}$

式中,

$\begin{array}{l}Y{}_{(t)}=[y{}_{(1)},y{}_{(2)},\cdots ,y{}_{(t)}]{}^{\text{Τ}},(4\text{a})\\ \mathit{\Lambda }{}_{(t)}=\left[\begin{array}{l}\mathit{\lambda }{}_{(1)}^{\text{Τ}}\\ \mathit{\lambda }{}_{(2)}^{\text{Τ}}\\ ⋮\\ \mathit{\lambda }{}_{(t)}^{\text{Τ}}\end{array}\right](4\text{b})\\ \mathit{{\rm M}}{}_{(t)}=\mathit{{\rm M}}-\frac{1}{t}{\displaystyle \sum _{i=1}^{t-1}\widehat{{\displaystyle \mathit{\Xi }}}}{}_{(i)}‚(4\text{c})\end{array}$

其中, α为权因子, 它决定了对先验知识的重视程度, $\widehat{{\displaystyle \mathit{\Xi }}}{}_{(i)}$是由i个采样值得到的参数估计值组成的向量。

定理1 使目标函数式 (4) 达到最小的参数估计值$\widehat{{\displaystyle \mathit{\Xi }}}{}_{(i)}$应满足下列正规方程:

$\begin{array}{l}\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{Y}{}_{(t)}+\alpha \mathit{{\rm M}}{}_{(t)}=\\ \left(\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{\Lambda }{}_{(t)}+\frac{\alpha }{t}\mathit{{\rm I}}\right)\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}。(5)\end{array}$

如果$\left(\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{\Lambda }{}_{(t)}+\frac{\alpha }{t}\mathit{{\rm I}}\right)$满秩可得唯一解

$\begin{array}{l}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}=(\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{Y}{}_{(t)}+\alpha \mathit{{\rm M}}{}_{(t)})\cdot \\ \left(\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{\Lambda }{}_{(t)}+\frac{\alpha }{t}\mathit{{\rm I}}\right){}^{-1}‚(6)\end{array}$

式 (6) 便是基于修正目标函数的时变参数LS估计算法。

定理1证明: 对目标函数式 (4) 求导数得。

$\begin{array}{l}\frac{\partial J}{\partial \widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}}=2\frac{\partial (\mathit{Y}{}_{(t)}-\mathit{\Lambda }{}_{(t)}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}){}^{\text{Τ}}}{\partial \widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}}\cdot \\ \mathit{W}{}_{(t)}(\mathit{Y}{}_{(t)}-\mathit{\Lambda }{}_{(t)}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)})+\\ 2\alpha t\frac{\partial \left(\mathit{{\rm M}}{}_{(t)}-\frac{1}{t}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}\right)}{\partial \widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}}\left(\mathit{{\rm M}}{}_{(t)}-\frac{1}{t}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}\right)=\\ -2\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}(\mathit{Y}{}_{(t)}-\mathit{\Lambda }{}_{(t)}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)})-\\ 2\alpha t\frac{1}{t}\left(\mathit{{\rm M}}{}_{(t)}-\frac{1}{t}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}\right)=\\ -2(\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{Y}{}_{(t)}-\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{\Lambda }{}_{(t)}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}+\\ \alpha \mathit{{\rm M}}{}_{(t)}-\frac{\alpha }{t}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)})=\\ -2[(\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{Y}{}_{(t)}+\alpha \mathit{{\rm M}}{}_{(t)})-\\ (\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{\Lambda }{}_{(t)}+\frac{\alpha }{t}\mathit{{\rm I}})\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}]。(7)\end{array}$

根据极值定理, 使目标函数式 (4) 最小的必要条件是$\partial J/\partial \widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}=0$, 即式 (5) , 又Heese矩阵为

J极小化的充分条件为

$\frac{\partial }{\partial \widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}}\left[\frac{\partial J}{\partial \widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}}\right]{}^{\text{Τ}}=2\left[\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{\Lambda }{}_{(t)}+\frac{\alpha }{t}\mathit{{\rm I}}\right]>0。(9)$

由于权因子α可调整, 故即使矩阵Λ^T _(t) 不满秩式 (9) 亦可成立且式 (6) 有唯一解。定理1得证。

推论1 在式 (6) 中, 当先验知识权因子α=0时, 则式 (6) 便是加权LS算法, 即:

$\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}=(\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{\Lambda }{}_{(t)}){}^{-1}\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{Y}{}_{(t)}‚(10)$

对应的目标函数为

$J=[\mathit{Y}{}_{(t)}-\mathit{\Lambda }{}_{(t)}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}]{}^{\text{Τ}}\mathit{W}{}_{(t)}[\mathit{Y}{}_{(t)}-\mathit{\Lambda }{}_{(t)}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}]。(11)$

推论2 在式 (6) 中, 当α=0且当W _(t) =I时, 则式 (6) 便是普通LS算法, 即

$\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}=(\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{\Lambda }{}_{(t)}){}^{-1}\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{Y}{}_{(t)}‚(12)$

对应的目标函数为

$J=[\mathit{Y}{}_{(t)}-\mathit{\Lambda }{}_{(t)}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}]{}^{\text{Τ}}[\mathit{Y}{}_{(t)}-\mathit{\Lambda }{}_{(t)}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}]。(13)$

从定理1中不难看出, 引入加权矩阵W _(t) 和权因子α后, 为获得时变参数的估计值和克服病态估计提供了更多的灵活性。

定理2 即使在过程式 (1) 持续激励条件不满足的情况下, 通过调整权因子α仍能保证式 (6) 有唯一解。

定理2说明:从理论上讲, 若过程式 (1) 输入信号不满足持续激励条件, 在计算中表现为矩阵Λ^T _(t) W _(t) Λ _(t) 不满秩, 而通过选择α, 则能保证矩阵$\left[\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{\Lambda }{}_{(t)}+\frac{\alpha }{t}\mathit{{\rm I}}\right]$满秩, 仍能使式 (6) 有唯一的正确的解。显而易见, 若用普通算法就会出现多解而得不到正确的估算值。

性质1 先验知识向量M对基于式 (4) 的最小二乘解的影响, 随着观测次数的增大而减小, 当时t→∞时, 式 (6) 趋向式 (10) , 即新算法收敛于加权LS算法。

性质1说明:从 (6) 式可以看出, 先验知识向量M对最小二乘解的影响, 不仅与权因子α有关, 而且与观测次数t有关。t较小时影响明显, 可有效地避免由于输入信号不持续激励而带来的问题;t太大时, 先验知识向量几乎不起作用, 保证收敛性。

性质1非常重要, 因为t较小时, 持续激励条件一般不易满足, 增大M的作用显然对避免病态估计十分有利;t太大时, 保证收敛性则是新算法成立的理论基础。

值得注意的是, 遗忘因子δ和权因子α要经过适当的选择。在进行估计的初期, 为了得到一个好的收敛速度, δ可选得小些, 而α可选得大些;在估计的后期, 为得到好的精度 (小方差) , δ可趋近于1, 而α可趋近于0。

利用新算法对时变过程式 (1) 进行参数估计, 显然希望得到其递推形式。而新算法不能直接给出递推形式, 需作近似处理。当t相对于α较大时, 数量矩阵$\frac{\alpha }{t}\mathit{{\rm I}}$相对于Λ^T _(t) W _(t) Λ _(t) 很小, 可以忽略。这时, 式 (6) 演变为

$\begin{array}{l}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}=(\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{Y}{}_{(t)}+\\ \alpha \mathit{{\rm M}}{}_{(t)})(\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{\Lambda }{}_{(t)}){}^{-1}。(14)\end{array}$

在推导递推公式时要频繁引用矩阵求逆引理。

引理1 设A+BC^T和I+C^TB^-1B满秩, 则下面的矩阵恒等式成立

$\begin{array}{l}(\mathit{A}+\mathit{B}\mathit{C}{}^{\text{Τ}}){}^{-1}=\\ \mathit{A}{}^{-1}-\mathit{A}{}^{-1}\mathit{B}(\mathit{{\rm I}}+\mathit{C}{}^{\text{Τ}}\mathit{A}{}^{-1}\mathit{B}){}^{-1}\mathit{C}{}^{\text{Τ}}\mathit{A}{}^{-1}。(15)\end{array}$

定理3 未知变量Ξ的基于修正目标函数的加权最小二乘估计$\widehat{{\displaystyle \mathit{\Xi }}}$的递推公式为

$\begin{array}{l}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t+1)}=\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}+\mathit{\Psi }{}_{(t+1)}(y{}_{(t+1)}-\mathit{\lambda }{}_{(t+1)}^{\text{Τ}}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)})+\\ \alpha \mathit{\Pi }{}_{(t+1)}(\mathit{{\rm M}}{}_{(t+1)}-\delta {}^2\mathit{{\rm M}}{}_{(t)})‚(16)\end{array}$

$\Psi {}_{(t+1)}=\frac{\mathit{\Pi }{}_{(t)}\mathit{\lambda }{}_{(t+1)}}{\delta {}^2+\mathit{\lambda }{}_{(t+1)}^{\text{Τ}}\mathit{\Pi }{}_{(t)}\mathit{\lambda }{}_{(t+1)}}$, (17)

$\mathit{\Pi }{}_{(t+1)}=\frac{1}{\delta {}^2}(\mathit{{\rm I}}-\mathit{\Psi }{}_{(t+1)}\mathit{\lambda }{}_{(t+1)}^{\text{Τ}})\mathit{\Pi }{}_{(t)}$。 (18)

定理3证明: 令

$\begin{array}{l}\mathit{\Lambda }{}_{(t+1)}=\left[\begin{array}{c}\mathit{\Lambda }{}_{(t)}\\ \mathit{\lambda }{}_{(t+1)}^{\text{Τ}}\end{array}\right]‚(19)\\ \mathit{Y}{}_{(t+1)}=\left[\begin{array}{c}\mathit{Y}{}_{(t)}\\ y{}_{(t+1)}\end{array}\right]‚(20)\\ \mathit{\Pi }{}_{(t)}=(\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{\Lambda }{}_{(t)}){}^{-1}‚(21)\\ \text{则}\mathit{\Pi }{}_{(t+1)}=(\mathit{\Lambda }{}_{(t+1)}^{\text{Τ}}\mathit{W}{}_{(t+1)}\mathit{\Lambda }{}_{(t+1)}){}^{-1}。(22)\end{array}$

由式 (2) 得

$\begin{array}{l}\mathit{W}{}_{(t+1)}=\left[\begin{array}{cc}\delta {}^2\left[\begin{array}{ccccc}\delta {}^{2(t-1)}& & & & \\ & \delta {}^{2(t-2)}& & 0& \\ & & \ddots & & \\ & 0& & \delta {}^2& \\ & & & & 1\end{array}\right]& 0\\ 0& 1\end{array}\right]=\\ \left[\begin{array}{cc}\delta {}^2\mathit{W}{}_{(t)}& 0\\ 0& 1\end{array}\right]。(23)\end{array}$

由式 (19) 、式 (20) 及式 (23) , 式 (22) 可变为

$\begin{array}{l}\mathit{\Pi }{}_{(t+1)}=\{\left[\mathit{\Lambda }{}_{(t)}^{\text{Τ}}┊\mathit{\lambda }{}_{(t+1)}\right]\cdot \\ \left[\begin{array}{cc}\delta {}^2\mathit{{\rm M}}{}_{(t)}& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}\mathit{\Lambda }{}_{(t)}\\ \mathit{\lambda }{}_{(t+1)}^{\text{Τ}}\end{array}\right]\}{}^{-1}=\\ (\delta {}^2\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{\Lambda }{}_{(t)}+\mathit{\lambda }{}_{(t+1)}\mathit{\lambda }{}_{(t+1)}^{\text{Τ}}){}^{-1}=\\ [\delta {}^2\mathit{\Pi }{}_{(t)}^{-1}+\mathit{\lambda }{}_{(t+1)}\mathit{\lambda }{}_{(t+1)}^{\text{Τ}}]{}^{-1}。(24)\end{array}$

由引理1, 式 (24) 可化为

$\begin{array}{l}\mathit{\Pi }{}_{(t+1)}=\frac{1}{\delta {}^2}(\mathit{{\rm I}}-\beta {}_{(t+1)}\mathit{\Pi }{}_{(t)}\mathit{\lambda }{}_{(t+1)}\mathit{\lambda }{}_{(t+1)}^{\text{Τ}})\mathit{\Pi }{}_{(t)}‚(25)\\ \text{式}\text{中}\beta {}_{(t+1)}=(\delta {}^2+\mathit{\lambda }{}_{(t+1)}^{\text{Τ}}\mathit{\Pi }{}_{(t)}\mathit{\lambda }{}_{(t+1)}){}^{-1}‚(26)\end{array}$

β (t+1) 是一个标量。

$\text{令}:\mathit{\Psi }{}_{(t+1)}=\mathit{\Pi }{}_{(t+1)}\mathit{\lambda }{}_{(t+1)},(27)$

将式 (25) 代入式 (26) 可得

$\begin{array}{l}\mathit{\Psi }{}_{(t+1)}=\frac{1}{\delta {}^2}\mathit{\Pi }{}_{(t)}\mathit{\lambda }{}_{(t+1)}\cdot \\ (\mathit{{\rm I}}-\beta {}_{(t+1)}\mathit{\Pi }{}_{(t)}\mathit{\lambda }{}_{(t+1)}^{\text{Τ}}\mathit{\Pi }{}_{(t)}\mathit{\lambda }{}_{(t+1)})=\\ \frac{\mathit{\Pi }{}_{(t+1)}\mathit{\lambda }{}_{(t+1)}}{\delta {}^2+\mathit{\lambda }{}_{(t+1)}^{\text{Τ}}\mathit{\Pi }{}_{(t)}\mathit{\lambda }{}_{(t+1)}}。(28)\end{array}$

式 (17) 得证。

由式 (28) , 式 (25) 可变为Π _(t+1) = (1/δ²) (I-Ψ _(t+1) λ^T _(t+1) ) Π _(t) , 式 (18) 得证。

另外, 由式 (19) 、式 (20) 及式 (23) 可得

$\begin{array}{l}\mathit{\Lambda }{}_{(t+1)}\mathit{W}{}_{(t+1)}\mathit{Y}{}_{(t+1)}=\\ \{[\mathit{\Lambda }{}_{(t)}^{\text{Τ}}┊\mathit{\lambda }{}_{(t+1)}]\left[\begin{array}{cc}\delta {}^2\mathit{{\rm M}}{}_{(t)}& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}\mathit{Y}{}_{(t)}\\ y{}_{(t+1)}\end{array}\right]\}=\\ \delta {}^2\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{Y}{}_{(t)}+\mathit{\lambda }{}_{(t+1)}\mathit{y}{}_{(t+1)}。(29)\end{array}$

由式 (29) 及式 (14) 得

$\begin{array}{l}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t+1)}=[\mathit{\Lambda }{}_{(t+1)}^{\text{Τ}}\mathit{W}{}_{(t+1)}\mathit{\Lambda }{}_{(t+1)}]{}^{-1}\cdot \\ [\mathit{\Lambda }{}_{(t+1)}^{\text{Τ}}\mathit{W}{}_{(t+1)}\mathit{Y}{}_{(t+1)}+\alpha \mathit{{\rm M}}{}_{(t+1)}]=\\ \mathit{\Pi }{}_{(t+1)}[\delta {}^2\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{Y}{}_{(t)}+\mathit{\lambda }{}_{(t+1)}y{}_{(t+1)}]+\\ \alpha \mathit{\Pi }{}_{(t+1)}\mathit{{\rm M}}{}_{(t+1)}。(30)\end{array}$

由式 (14) 及式 (21) 得

$\begin{array}{l}\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{Y}{}_{(t)}=(\mathit{\Lambda }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\mathit{\Lambda }{}_{(t)})\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}-\\ \alpha \mathit{{\rm M}}{}_{(t)}=\mathit{\Pi }{}_{(t+1)}^{-1}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}-\alpha \mathit{{\rm M}}{}_{(t)}。(31)\end{array}$

将式 (31) 代入式 (30) 可得

$\begin{array}{l}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t+1)}=\mathit{\Pi }{}_{(t+1)}(\delta {}^2\mathit{\Pi }{}_{(t)}^{-1}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}+\mathit{\lambda }{}_{(t+1)}y{}_{(t+1)})+\\ \alpha \mathit{\Pi }{}_{(t+1)}(\mathit{{\rm M}}{}_{(t+1)}-\delta {}^2\mathit{{\rm M}}{}_{(t)})(32)\end{array}$

由式 (24) 得

$\delta {}^2\mathit{\Pi }{}_{(t)}^{-1}=\mathit{\Pi }{}_{(t+1)}^{-1}-\mathit{\lambda }{}_{(t+1)}\mathit{\lambda }{}_{(t+1)}^{\text{Τ}}。(33)$

将式 (33) 代入式 (32) 有

$\begin{array}{l}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t+1)}=\mathit{\Pi }{}_{(t+1)}[\mathit{\Pi }{}_{(t+1)}^{-1}-\mathit{\lambda }{}_{(t+1)}\mathit{\lambda }{}_{(t+1)}^{\text{Τ}}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}+\\ \mathit{\lambda }{}_{(t+1)}y{}_{(t+1)}]+\alpha \mathit{\Pi }{}_{(t+1)}(\mathit{{\rm M}}{}_{(t+1)}-\delta {}^2\mathit{{\rm M}}{}_{(t)})=\\ \widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}+\mathit{\Pi }{}_{(t+1)}\mathit{\lambda }{}_{(t+1)}(y{}_{(t+1)}-\mathit{\lambda }{}_{(t+1)}^{\text{Τ}}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)})+\\ \alpha \mathit{\Pi }{}_{(t+1)}(\mathit{{\rm M}}{}_{(t+1)}-\delta {}^2\mathit{{\rm M}}{}_{(t)})。(34)\end{array}$

将式 (27) 代入式 (34) 得

$\begin{array}{l}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t+1)}=\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}+\mathit{\Psi }{}_{(t+1)}(y{}_{(t+1)}-\mathit{\lambda }{}_{(t+1)}^{\text{Τ}}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)})+\\ \alpha \mathit{\Pi }{}_{(t+1)}(\mathit{{\rm M}}{}_{(t+1)}-\delta {}^2\mathit{{\rm M}}{}_{(t)})。(35)\end{array}$

定理3即式 (16) 证毕。

在定理3中, 新递推算法的初值一般可令$\widehat{{\displaystyle \mathit{\Xi }}}{}_{(0)}=0,\mathit{\Pi }{}_{(0)}=\sigma {}^2\mathit{{\rm I}}$, 其中σ为充分大的正数。

新递推算法的步骤如下:

Step 1 置初值$\widehat{{\displaystyle \mathit{\Xi }}}{}_{(0)}‚\mathit{\Pi }{}_{(0)}$, 输入初始数据;

Step 2 采样当前输入和输出;

Step 3 按式 (16) 、式 (17) 及式 (18) 计算$\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}‚\mathit{\Psi }{}_{(t)}‚\mathit{\Pi }{}_{(t)}$;

Step 4 返回Step 2 直到收敛满足要求为止。

根据数理统计及LS算法的知识, 不难发现, 当过程式 (1) 中的$\tilde{\mathit{\omega }}$ _(t) 为有色噪声时, 定理1及定理3所给出的算法同普通LS算法一样, 其估计值$\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}$不能保证是一致估计。为此, 这里仿照普通最小二乘辅助变量算法可以得到基于修正目标函数的时变参数辅助变量算法。

定理4 引入辅助变量Σ _(t) , 并选择加权矩阵W _(t) 使其满足

$\mathit{\Sigma }{}_{(t)}=\mathit{W}{}_{(t)}\mathit{\Lambda }{}_{(t)}‚(36)$

如果Σ _(t) 满足以下两个条件:

$\begin{array}{l}\mathit{E}\{\mathit{\Sigma }{}_{(t)}^{\text{Τ}}\mathit{W}{}_{(t)}\}=0‚(37)\\ \mathit{E}\{\mathit{\Sigma }{}_{(t)}^{\text{Τ}}\mathit{\Lambda }{}_{(t)}\}=\mathit{\Gamma }‚\mathit{\Gamma }\text{非}\text{奇}\text{异}‚(38)\end{array}$

那么, 基于修正目标函数的时变参数辅助变量算法可表示为

$\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}=\left[\mathit{\Sigma }{}_{(t)}^{\text{Τ}}\mathit{\Lambda }{}_{(t)}+\frac{\alpha }{t}\mathit{{\rm I}}\right]{}^{-1}[\mathit{\Sigma }{}_{(t)}^{\text{Τ}}\mathit{Y}{}_{(t)}+\alpha \mathit{{\rm M}}{}_{(t)}]。(39)$

定理4证明从略。

同样地将式 (39) 稍作近似处理便可得到基于修正目标函数的辅助变量递推算法。

定理5 未知变量Ξ的基于修正目标函数的辅助变量递推公式为

$\begin{array}{l}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t+1)}=\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}+\mathit{\Psi }{}_{(t+1)}[y{}_{(t+1)}-\mathit{\lambda }{}_{(t+1)}^{\text{Τ}}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}]+\\ \alpha \mathit{\Pi }{}_{(t+1)}[\mathit{{\rm M}}{}_{(t+1)}-\delta {}^2\mathit{{\rm M}}{}_{(t)}]‚(40)\end{array}$

$\mathit{\Psi }{}_{(t+1)}=\frac{\mathit{\Pi }{}_{(t)}\mathit{\zeta }{}_{(t+1)}}{\delta {}^2+\mathit{\lambda }{}_{(t+1)}^{\text{Τ}}\mathit{\Pi }{}_{(t)}\mathit{\zeta }{}_{(t+1)}}$, (41)

$\mathit{\Pi }{}_{(t+1)}=\frac{1}{\delta {}^2}[\mathit{{\rm I}}-\mathit{\Psi }{}_{(t+1)}\mathit{\lambda }{}_{(t+1)}^{\text{Τ}}]\mathit{\Pi }{}_{(t)}$。 (42)

式中ζ _(t) =[-υ _(t-1) , -υ _(t-2) , …, -υ _(t-na) , u _(t-1) , u _(t-2) , …, u _(t-nb) ], 可取$\upsilon {}_{(t)}=\mathit{\Sigma }\text{Τ}{}_{(t)}\widehat{{\displaystyle \mathit{\Xi }}}{}_{(t)}$。

定理5证明从略。

另外, 新算法也可推广到广义LS算法, 在此暂不作研究。

分别采用LS算法和基于修正目标函数的LS算法来估计时变系统, 采用的时变系统是水泥生料混合过程。

水泥生料混合控制历来是水泥生产的关键问题。图1示出水泥混合系统的流程。计算机根据X荧光分析仪检测数据分析计算后, 与给定值比较, 控制电子皮带秤, 调节原料的配比。

如果假设系统的工作点稳定, 且混合过程理想, 没有死区时间, 那么水泥生粉中氧化物的成分与生料重量的关系如下:

$\left(\begin{array}{l}{\rm O}{}_{\text{C}\text{a}\text{Ο}}\\ {\rm O}{}_{\text{S}\text{i}\text{Ο}{}_2}\\ {\rm O}{}_{\text{A}\text{l}{}_2\text{Ο}{}_3}\\ {\rm O}{}_{\text{F}\text{e}{}_2\text{Ο}{}_3}\end{array}\right)=\left(\begin{array}{cccc}C{}_1& C{}_2& {\rm O}& {\rm O}\\ S{}_1& S{}_2& S{}_3& {\rm O}\\ A{}_1& A{}_2& {\rm O}& {\rm O}\\ F{}_1& F{}_2& {\rm O}& F{}_4\end{array}\right)\left(\begin{array}{l}\omega {}_1\\ \omega {}_2\\ \omega {}_3\\ \omega {}_4\end{array}\right)‚(43)$

或用下列向量形式表示:

$\mathit{{\rm O}}=\mathit{C}\mathit{\omega }(44)$

式中, 矩阵C为生料成分矩阵;O为生粉中氧化物的含量;ω表示生料中石灰石、粘土、石英砂和铁矿石的重量比。

由于水泥原料 (石灰石、粘土、石英砂和铁矿石等) 大多呈块状、分散性大, 直接实时测量其成分是不可能的, 只能根据磨机输出的矩阵O的生粉化学成分和原料进料量ω来估计生料成分C, 从而为控制算法提供重要依据。这里以成分波动较大的石灰石中的SiO₂含量为例给出了其对比仿真结果。

图2表明:基于修正目标函数的整匹估计算法的估计精度优于普通LS算法, 引入遗忘因子δ以后, 新算法对时变参数的跟踪能力大为提高, 精度也有所改善。图3表明:基于修正目标函数的递推算法比普通递推LS算法对时变参数的估计精度大为提高。

本算法是多年以来在水泥生料自适应控制系统研究的基础上经不断修正总结而来的。该算法结合自校正控制技术用于工程实际, 已将磨机输出生粉成分方差降低30 %左右。工业现场使用结果表明:所提出的系统具有较高的鲁棒性和控制精度, 可将磨机输出生粉成分方差降低到1.58, 而且长期运行稳定。

导出了一类基于修正目标函数的时变参数LS算法, 克服普通算法使病态程度加剧的缺陷, 并提高对时变参数的跟踪能力;避免数据饱和现象, 可提高自适应控制系统的鲁棒性。这种算法还具有运算量小、收敛性好的优点, 十分便于工程实现和应用。已取得良好的控制效果。此外, 以所提出的算法为基础, 还易于构成相应的递推辅助变量算法、广义LS算法等, 还可以构成对数据加权处理的限定记忆估计算法。