自适应控制-系统辨识部分-Part1

自适应控制

00-自适应控制引言（绪论1）

• 参考书目
  • Adaptive Control(2-edition),K.J.Astrom
  • System Identification-Theory for the User(2-Edition),Lennart Ljung
  • 过程辨识，方崇智

01-自适应控制的发展概况（绪论2）

一、什么是自适应控制器

• 自适应控制器定义：1）具有“控制器结构或参数整定机构”的控制器；2）相对于定常参数控制器而言，“结构或者参数”能够在线调整的控制器。

二、为什么要改变“控制器的参数或结构”

1. 不改不行：被控对象或者环境发生变化，控制律无法使系统稳定；
2. 参数或结构存在不确定性，需要进一步修正控制律以提升控制性能。

三、自适应控制适用的被控对象

1. 参数、结构不确定系统（线性、非线性）：如 模型参考自适应系统
2. 参数未知或随时间“缓慢”变化的被动对象：如 随机自适应系统
3. 非线性系统：如 增益调度自适应系统，T-S模糊系统

四、自适应控制的发展概况

五、自适应控制理论的发展动力

六、自适应控制理论的发展方向

七、自适应控制系统的分类

1. 自整定调节器（Self-tuning Regulators,STR）

   
   • 多种自适应方案：不同设计方法+不同辨识方法
   • 没有基于稳定性的考虑，缺乏“从顶向下”的整体设计方案

2. 双重控制（Dual Control）

3. 模型参考自适应控制系统（Model-Reference Adaptive System）

   

4. 增益调度系统（Gain Scheduling）

   

5. 无模型自适应控制系统

   • 定义：控制器的设计仅利用受控系统的输入输出数据，控制器中不包含受控过程数学模型的任何信息的控制理论与方法。
   • 方法：伪梯度向量、迭代无模型控制优化方法（迭代反馈自校正方法）、去伪控制、迭代学习控制和重复控制、强化学习自适应控制

02-系统辨识的基本概念（参数在线估计）

一、什么是系统辨识（参数估计、学习）

• 系统辨识是根据系统的输入输出时间函数来确定描述系统行为的数学模型。即按照一定的准则，从某一类模型中找出一个与输入输出数据拟合得最好的模型

1. 机理分析建模方法（白箱法）

   • 问题：效率低，物理参数需进一步确定，不方便计算机在线决策

2. 系统辨识法（黑箱法）

   1. 输入信号为什么要选M序列？

      M序列的自相关性较好,具有伪随机性,容易产生和复制。正余弦函数的信息不够丰富。阶跃信号不易产生。

   2. 预定的模型结构、阶次如何选定？

   3. 怎么确定具体参数？

   4. 如何通过迭代的方法在线估计系统的参数？

3. 机理分析法+系统辨识法（灰箱法）

二、系统辨识方法的基本分类

$$\{\begin{array}{l}参数辨识方法\{\begin{array}{l}a.经典辨识方法\\ \text{b.最小二乘类参数辨识方法}\\ \text{c.基于梯度信息的参数辨识方法}\\ d.极大似然法和预报误差方法\\ ...\end{array}\\ 结构辨识方法\{\begin{array}{l}\text{a.根据Hankel矩阵的秩估计模型的阶次}\\ b.行列式比\\ c.残差的方差\\ ...\end{array}\end{array}$$

三、辨识的基本要素

$$\{\begin{array}{l}\text{1.输入输出数据（辨识的基础）}\{\begin{array}{l}必须包含有关系统特性的足够信息\\ \text{时域的角度：}信号变化强烈，且呈现非周期性\\ \text{频域的角度：}频谱宽\end{array}\\ \text{2.模型类}\\ \text{3.准则：}评判“辨识得到的模型”是否满足“实际需要”的一个“准则”\end{array}$$

辨识就是按照一定的准则从某一类模型中找出一个与输入输出数据拟合得最好的模型。

四、模型简介（Book.1 Page.81）

1. ARX模型

   • 有源自回归模型

   • 分量形式
     $$\begin{gathered} y(k)+a_{1}y(k-1)+\cdots +a_{na}y(k-na)=b_{1}u(k-1)+\cdots +b_{nb}u(k-nb)+e(k) \\ e(k)为白噪声 \end{gathered}$$

   • 矩阵形式
     $$\begin{gathered} A(z)y(k)=B(z)u(k)+e(k) \\ A(z)=1+a_1{z}^{-1}+...+a_{na}{z}^{-na} \\ B(z)=b_1{z}^{-1}+...+b_{nb}{z}^{-nb} \end{gathered}$$

   • Figure

     

2. ARMAX模型

   • MA：滑动平均

   • 矩阵形式
     $$\begin{gathered} A(z)y(k)=B(z)u(k)+C(z)e(k) \\ C(z)=1+c_1{z}^{-1}+...+c_{nc}{z}^{-nc} \end{gathered}$$

   • Figure

     

3. ARARX模型

   • 矩阵形式
     $$\begin{gathered} A(z)y(k)=B(z)u(k)+\frac{1}{D(z)}e(k) \\ D(z)=1+d_1{z}^{-1}+...+d_{nd}{z}^{-nd} \end{gathered}$$

   • Figure

     

4. ARARMAX模型

   • 矩阵形式
     $$A(z)y(k)=B(z)u(k)+\frac{C(z)}{D(z)}e(k)$$

5. Output Error Model Structure（OE model structure）

6. Box-Jenkins Model Structure

7. A General Family of Model Structures

8. 神经网络模型（Artificial Neural Networks）

9. 基于规则的模型（Rule based Model）

10. 基于模糊推理系统的自适应网络（Adaptive-Network based Fuzzy Inference System,ANFIS）

03-最小二乘参数辨识标准算法（最小二乘类1）

一、引例

• 总结：
  1. 以“测量值与估计值残差的平方和最小”指标，在“线性函数”这一类函数中找到最优解；
  2. 如果函数类扩大到“非线性”方程，有可能找到一个使“偏差的平方和”更小的方程；
  3. 多元函数优化问题

二、基本定义和基本计算公式

1. 最小二乘模型类和标准格式(ARX模型)

   • 分量形式:
     $$y(k)+a_{1}y(k-1)+...+a_{n}y(k-n)=b_{1}u(k-1)+...+b_{n}u(k-n)+e(k)$$

   • 向量形式:

   $$y(k)=\left[\begin{array}{cccccc}-y(k-1)& \cdots & -y(k-n)& u(k-1)& \cdots & u(k-n)\end{array}\right]\left[\begin{array}{c}{a}_1\\ ⋮\\ a_n\\ b_1\\ ⋮\\ b_n\end{array}\right]$$

2. 最小二乘准则（最小二乘估计原理）和正规方程组

   • 准则：$$J=\sum _{i=1}^N[y(n+i)-\hat{y}(n+i){]}^2$$

   • 正规方程组：$${\Phi }^T(N)Y(N)-{\Phi }^T(N)\Phi (N)\theta (N)=0$$

     其中：
     $$\begin{gathered} \Phi (N)=\left[\begin{array}{cccccccc}-y(n)& -y(n-1)& \cdots & -y(1)& u(n)& u(n-1)& \cdots & u(1)\\ -y(n+1)& -y(n)& \cdots & -y(2)& u(n+1)& u(n)& \cdots & u(2)\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮\\ -y(n+N-1)& -y(n+N-2)& \cdots & -y(N)& u(n+N-1)& u(n+N-2)& \cdots & u(N)\end{array}\right] \\ \theta (N)={\left[\begin{array}{cccccc}{a}_n& \cdots & a_1& u_n& \cdots & u_1\end{array}\right]}^T \\ Y(N)={\left[\begin{array}{ccc}y(n+1)& \cdots & y(n+N)\end{array}\right]}^T \end{gathered}$$

3. 参数解的表达式：

$$\hat{\theta }(N)=[{\Phi }^T(N)\Phi (N){]}^{-1}{\Phi }^T(N)Y(N)$$

4. $[{\Phi }^T(N)\Phi (N){]}_{2n\times 2n}$ 的非奇异性

   • 从矩阵理论上讲，非奇异要求 $\Phi (N)$ 的行数 $N\ge 2n$；从物理背景上看，$N>2n$；
   • 当 $N > 2n $ 时，$[{\Phi }^T(N)\Phi (N)]$ 不一定非奇异。

5. 开环可辨识性条件

   • 充分必要条件：输入信号必须是 $2n$ 阶持续激励信号，即要求：
     $$\begin{gathered} {\overline{U}}_L^T{\overline{U}}_L>0 \\ \{\begin{array}{l}{\overline{U}}_L=[F{u}_L,F^2{u}_L,\dots ,F^{2n}{u}_L]\\ u_L=[u(1),u(2),\dots ,u(L){]}^T\\ F=\left[\begin{array}{ccc}0& & \mathrm{O}\\ 1& \ddots & \\ \mathrm{O}& 1& 0\end{array}\right]\\ n=\mathbf{max}(n_a,n_b)\end{array} \end{gathered}$$

   • 常用的信号：白噪声序列；M序列或逆M序列；n种频率的正弦信号组合成的信号，且其频率不成整数倍关系

三、算法程序示例

四、算法演示

1. 例1：稳定的被控对象

   1. 系统函数：
      $$\begin{gathered} G(z)=\frac{0.45z+0.23}{z^2-1.81z+0.82} \\ y(k)-1.81y(k-1)+0.82y(k-2)=0.45u(k-1)+0.23u(k-2) \end{gathered}$$

   2. MATLAB程序：C03P19_example1.m：

      1. 向量 $\phi$：
         $$\phi (1)=[\begin{array}{cccc}-y(2)& -y(1)& u(2)& u(1)\end{array}]$$

      2. 矩阵 $\Phi$ 的赋值：

         Phi =zeros(N,2*n);
         for i=1:N
             for j=1:n
                 Phi(i,j)=-y(n+i-j);
             end
             for j=1:n
                 Phi(i,n+j)=u(n+i-j);
             end
         end
         

2. 例2：临界稳定的被控对象

   1. 系统函数：
      $$\begin{gathered} G(z)=\frac{0.45z+0.23}{z^2-1.82z+0.82} \\ y(k)-1.82y(k-1)+0.82y(k-2)=0.45u(k-1)+0.23u(k-2) \end{gathered}$$

3. 例3：不稳定的被控对象

   1. 系统函数
      $$\begin{gathered} G(z)=\frac{0.45z+0.23}{z^2-1.83z+0.82} \\ y(k)-1.83y(k-1)+0.82y(k-2)=0.45u(k-1)+0.23u(k-2) \end{gathered}$$

4. 仿真结论：

   1. 最小二乘既可以用到“稳定对象”，也可以用到”不稳定对象“；辨识自身不稳定对象，应注意保持系统在平衡位置附近；

   2. 噪声的方差对辨识精度影响很大，其影响程度对不同的被控对象有所不同。

五、加权最小二乘

1. 指标：$$J=\sum _{i=1}^N\lambda (i){\left[y(n+i)-\hat{y}(n+i)\right]}^2$$
2. 计算公式：

$$\begin{gathered} \hat{\theta }(N)={\left[{\Phi }^T(N)\Lambda (N)\Phi (N)\right]}^{-1}{\Phi }^T(N)\Lambda (N)y(N) \\ \Lambda (N)=diag([\lambda (1),\lambda (2),\dots ,\lambda (N)])>0 \end{gathered}$$

3. 在噪声信号为平稳随机序列的前提下，最小二乘指标可以看成是对最小方差指标的一种近似逼近：

$$J=\sum _{i=1}^N{\left[y(n+i)-\hat{y}(n+i)\right]}^2\to J=E{\left[y(n+i)-\hat{y}(n+i)\right]}^2$$

六、带“初始估计”的最小二乘

1. 指标：$$J=\sum _{i=1}^N{\left[y(n+i)-\hat{y}(n+i)\right]}^2+(\theta -{\theta }_0{)}^{T}R(\theta -{\theta }_0)，其中{\theta }_0是对\theta 的初步猜测；R>0$$

2. 计算公式：

$$\theta (N)={\left[{\Phi }^T(N)\Phi (N)+R\right]}^{-1}({\Phi }^T(N)Y(N)+R{\theta }_0)$$

3. 问题：

   1. 假设被控对象是一个线性定常系统，什么情况下能通过系统辨识获得系统的真实参数？

   2. 辨识的精度与加到系统中的噪声干扰有何联系？如何从理论上说明这一点？

   3. 实际数据不免受到噪声的影响，如何对付这些噪声对输入输出数据的影响？

   4. 如何针对不稳定的被控对象展开参数辨识？

      尽量保证在稳定点附近辨识。

04-递推最小二乘参数估计算法（最小二乘类2）

一、问题的提出

1. 基本思想：

$$新的估计值\hat{\theta }(k)=老的估计值\hat{\theta }(k-1)+修正项$$

2. 建立数学模型：$$x(k)=a+e(k)，x(k)为第k次测量，a表示真实距离，e(k)为高斯白噪声$$
   • 指标：$$J=\sum _{k=1}^N{\left[x(k)-\hat{x}(k)\right]}^2$$
   • 批处理求解：$$\hat{a}(k)=\frac{{\sum }_{k-1}^N[x(k)]}{N}$$
   • 递推求解：$$\hat{a}(k)=\hat{a}(k-1)+\frac{1}{k}\left[x(k)-\hat{a}(k-1)\right]$$
   • 对初值不敏感

二、递推算法结果

$$\begin{array}{l}K(N+1)=P(N)\phi (N+1){\left[1+{\phi }^T(N+1)P(N)\phi (N+1)\right]}^{-1}\\ \hat{\theta }(N+1)=\hat{\theta }(N)+K(N+1)\left[y(n+N+1)-{\phi }^T(N+1)\hat{\theta }(N)\right]\\ P(N+1)=P(N)-K(N+1)K^T(N+1){\left[1+{\phi }^T(N+1)P(N)\phi (N+1)\right]}^{-1}\end{array}$$

$$\{\begin{array}{l}\hat{\theta }(N+1):2n\times 1维向量\\ {\phi }^T(N+1)=\left[-y(n+N),-y(n+N-1),\dots ,-y(N+1),u(n+N),u(n+N-1),\dots ,u(N+1)\right]:1\times 2n维向量\\ P(N):2n\times 2n维矩阵\end{array}$$

三、递推算法计算流程

$$\hat{\theta }(0),P(0)\stackrel{y(1),u(1)}{⟶}K(1),\hat{\theta }(1),P(1)\stackrel{y(2),u(2)}{⟶}K(2),\hat{\theta }(2),P(2)\stackrel{y(3),u(3)}{⟶}\cdots \stackrel{y(n),u(n)}{⟶}K(n),\hat{\theta }(n),P(n)$$

四、初值的选取方法

1. 用最小二乘的封闭形式计算：$$P_0(N_0),{\hat{\theta }}_0(N_0)$$
2. 认为给定初值：$$P_0={\epsilon }^{2}I,\epsilon 足够大;{\hat{\theta }}_0=0$$

五、递推算法的仿真实验

1. 例1：稳定被控对象

   1. 系统函数：
      $$\begin{gathered} G(z)=\frac{1.0z+0.5}{z^2-1.5z+0.7} \\ y(k)-1.5y(k-1)+0.7y(k-2)=1.0u(k-1)+0.5u(k-2) \end{gathered}$$

   2. 初值选择：
      $$P(0)=10^{10}I;\hat{\theta }(0)=0;噪声方差\sigma =0.01$$

   3. MATLAB程序：C04P26_example1.m

      1. 向量 $\phi (N+1)$ ：
         $$\phi (1)=[\begin{array}{cccc}-y(2)& -y(1)& u(2)& u(1)\end{array}]$$

      2. 问题：前几个值和最后几个值怎么办（原因：由于出现 $y(n+N+1)$ 项）

2. 例2：不稳定被控对象

   1. 系统函数：
      $$\begin{gathered} G(z)=\frac{0.45z+0.23}{z^2-1.83z+0.82} \\ y(k)-1.83y(k-1)+0.82y(k-2)=0.45u(k-1)+0.23u(k-2) \end{gathered}$$

3. 问题：

   • 为什么进一步增大迭代次数时计算发散了？
     • （猜测答案：系统本身发散了）
   • 对于不稳定的被控对象进行参数辨识时，应该采取什么措施？
     • （猜测答案：采用适当的控制策略来稳定被控对象）

六、加权最小二乘递推算法

$$\begin{array}{l}K(N)=P(N-1)\phi (N){\left[{\lambda }^{-1}(N)+{\phi }^T(N)P(N-1)\phi (N)\right]}^{-1}\\ \hat{\Theta }(N)=\hat{\Theta }(N-1)+K(N)\left[y(n+N)-{\phi }^T(N)\hat{\Theta }(N-1)\right]\\ P(N)=P(N_1)-K(N){\phi }^T(N)P(N-1)\end{array}$$

七、带初值估计的最小二乘递推算法

$$\begin{array}{rl}\hat{\theta }(N+1)& =\hat{\theta }(N)-P(N+1)\phi (N+1)\left[y(n+N+1)-{\phi }^T(N+1)\hat{\theta }(N)\right]\\ P(N+1)& =P(N)-\frac{P(N)\phi (N+1){\phi }^T(N+1)P(N)}{1+{\phi }^T(N+1)P(N)\phi (N+1)}\\ \theta (0)& ={\theta }_0,P(0)=R^{-1}\end{array}$$

八、递推算法和卡尔曼滤波之间的关系

1. 辨识问题转化为滤波问题：

$$\begin{gathered} \theta (k+1)=A\theta (k) \\ z(k)={\phi }^T(k)\theta (k)+e(k) \end{gathered}$$

2. 滤波算法：

$$\begin{array}{l}\hat{\theta }(k)=\hat{\theta }(k-1)+K(k)\left[z(k)-{\phi }^T(k)\hat{\theta }(k-1)\right]\\ K(k)=P^{-}(k)\phi (k){\left[{\phi }^T(k)P^{-}(k)\phi (k)+{\sigma }_n^2\right]}^{-1}\\ P^{-}(k)=AP(k-1)A+Qc\\ P(k)=\left[I-K(k){\phi }^T(k)\right]P^{-}(k)\\ \hat{\theta }(0)=0,P(0)=P_0\end{array}$$

附录1、递推算法的推导过程

05-渐消记忆法与限定记忆法最小二乘（最小二乘类3）

一、数据饱和现象

• 直观上，随着采集到的数据越来越多，递推最小二乘法应该给出更精确的参数估计值；但实际上，随着迭代次数增加，“估计值”与“真实值”的偏差往往会越来越远。

二、数据饱和现象的原因

• 随着递推次数的增加，旧的数据会累积的越来越多，造成结果是把新的数据信息淹没，最终导致参数估计无法进行，算法失效。

三、解决的方案

1. 渐消记忆法：增加新数据在计算中的权重，减小老数据在计算中的权重；
2. 限定记忆法：去掉一部分老数据。

四、渐消记忆算法流程

$$\begin{array}{l}K(N+1)=P(N)\phi (N+1){\left[\lambda +{\phi }^T(N+1)P(N)\phi (N+1)\right]}^{-1}\\ \hat{\Theta }(N+1)=\hat{\Theta }(N)+K(N+1)\left[y(n+N+1)-{\phi }^T(N+1)\hat{\Theta }(N)\right]\\ P(N+1)=\frac{1}{\lambda }\left\{P(N)-K(N+1)K^T(N+1)\left[\lambda +{\phi }^T(N+1)P(N)\phi (N+1)\right]\right\}\\ 其中，\lambda \in (0,1]\end{array}$$

五、渐消记忆算法推导过程

六、仿真分析

1. 仿真结论：
   1. 渐消记忆法在一定程度上克服了“数据饱和问题”，该算法对于“慢时变系统的辨识问题”比一般递推算法更有效；
   2. 遗忘因子对辨识精度有显著影响，需要适当选取。
2. 问题：
   1. 能否从动态系统的角度考虑为什么遗忘因子会导致很大的震动？

七、限定记忆算法流程

$$\begin{array}{l}P(k+1)=Q(k)-Q(k)\phi (k+1){\left[Q(k)\phi (k+1)\right]}^T{\left[1+{\phi }^T(k+1)Q(k)\phi (k+1)\right]}^{-1}\\ Q(k+1)=P(k+1)+P(k+1)\phi (k+1-N){\left[P(k+1)\phi (k+1-N)\right]}^T{\left[1-{\phi }^T(k+1-N)P(k+1)\phi (k+1-N)\right]}^{-1}\\ \hat{\theta }(k+1)=\hat{\theta }(k)+Q(k+1)\left[\phi (k+1)\Delta y(n+k+1)-\phi (k+1-N)\Delta y(n+k+1-N)\right]\end{array}$$

$$\begin{array}{lc}其中，& \Delta y(n+k+1)=y(n+k+1)-{\phi }^T(k+1)\hat{\theta }(K)\\ & \Delta (n+k+1-N)=y(n+k+1-N)-{\phi }^T(k+1-N)\hat{\theta }(k)\end{array}$$

八、限定记忆算法推导过程

九、限定记忆最小二乘算法程序流程

1. 先用递推算法根据前 $n+N$ 个数据估计 $\hat{\theta }(N)$：

$$\hat{\theta }(0),Q(0)\stackrel{y(1),u(1)}{\to }K(0),\hat{\theta }(1),P(1)\stackrel{y(2),u(2)}{\to }K(1),\hat{\theta }(2),P(2)\stackrel{y(3),u(3)}{\to }\cdots \stackrel{y(n),u(n)}{\to }K(n-1),\hat{\theta }(n),P(n)$$

2. 将 $Q(N)=P(N)$ 作为限定记忆算法的初值，开始调用限定记忆算法：

$$\stackrel{y(n+N+1),u(n+N+1)}{\to }\Delta y(n+k+1),\Delta (n+k+1-N)⟶P(k+1),Q(k+1),\hat{\theta }(k+1)$$

十、仿真

1. 仿真结论：限定记忆法可以有效改善数据饱和现象，可在很大程度上改善递推最小二乘算法的性能。

2. 问题：如何确定 $N$？

附录1、数据饱和现象的原因

附录2、渐消记忆法证明

附录3、限定记忆法证明

06-最小二乘解的几何意义及其统计特性（最小二乘类4）

一、引例

• 正交投影的定义：设 $y$ 是具有前二阶矩的 $n$ 维随机向量，$X$ 是适当维数的随机矩阵。如果存在一个与 $y$ 同维的随机向量 $y^{\ast }$，并且具备以下三个性质：

  1. $y^{\ast }$ 可以由 $X$ 线性表示，即：$y^{\ast }=a+Xb$
  2. 无偏性，即：$E[y^{\ast }]=E[y]$
  3. $y-y^{\ast }$ 与 $X$ 相互正交，即：$E[(y-y^{\ast }{)}^{T}X]=0$
     则称 $y^{\ast }$ 是 $y$ 在空间 $X$ 上的正交投影。

二、最小二乘估计的几何解释

• 如果噪声向量 $E(N)$ 的均值为零且与 $\Phi (N)$ 统计独立，则输出估计向量 $\hat{Y}(N)$ 是输出测量向量 $Y(N)$ 在由 $\alpha (1),\alpha (2),\dots ,\alpha (2n)$ 张成的空间上的正交投影，或者说输出残差向量 $\xi (N)$ 垂直于由 $\alpha (1),\alpha (2),\dots ,\alpha (2n)$ 张成的空间。

三、最小二乘参数估计值的统计性质

1. 无偏性

   1. 物理意义：估计值是否围绕真值波动。

   2. 定理：对于如下模型
      $$Y(N)=\Phi (N){\theta }_0+E(N)$$
      如果噪声向量 $E(N)$ 的均值为零，并且 $E(N)$ 和 $\Phi (N)$ 是统计独立的，则最小二乘估计值 $\hat{\theta }(N)$=$ 是无偏估计量，即
      $$E[\hat{\theta }]={\theta }_0$$
      其中，${\theta }_0$ 表示真实参数。

   3. 无偏性证明

   4. 问题：什么情况下 $E(N)$ 和 $\Phi (N)$ 是统计独立的？

      • 静态线性系统的回归分析中；
      • 动态系统的参数辨识中，一般情况下 $E(N)$ 和 $\Phi (N)$ 是相关的，最小二乘估计值 $\hat{\theta }(N)$ 是有偏估计；对于高斯白噪声序列，最小二乘估计具有无偏性。

2. 参数估计偏差的协方差性质

   1. 定理：对于如下模型
      $$Y(N)=\Phi (N){\theta }_0+E(N)$$
      如果噪声向量 $E(N)$ 的均值为零，协方差矩阵为 ${\sigma }^{2}I$ ，并且 $E(N)$ 和 $\Phi (N)$ 是统计独立的，则最小二乘参数估计偏差的协方差阵为
      C o v { θ 0 − θ ^ ( N ) } = σ 2 E { ( Φ T ( N ) Φ ( n ) ) − 1 } Cov\{ \theta_0- \hat\theta(N) \}= \sigma^2 E\{ (\Phi^T(N)\Phi(n))^{-1} \} Cov{θ0​−θ^(N)}=σ2E{(ΦT(N)Φ(n))−1}
      其中，${\theta }_0$ 表示真实参数。

   2. 在以上条件下，如果方差 ${\sigma }^2$ 未知，则下式给出方差的一个无偏估计：
      $$S^2=\frac{1}{N-2n}(Y(N)-\Phi (N)\hat{\theta }{)}^T(Y(N)-\Phi (n)\hat{\theta })$$
      其中，$n$ 为系统的阶次。

3. 一致性

   1. 物理意义：估计值将以概率 $1$ 收敛于真值。当 $$N$ 很大时，工程上的计算值是可以接受的。

   2. 定理：对于如下模型
      $$y(k)=-a_{1}y(k-1)-a_{2}y(k-2)-\cdots -a_{na}y(k-na)+b_{1}u(k-1)+b_{2}u(k-2)+\cdots +b_{nb}u(k-nb)+e(k)$$
      当满足以下四个条件

      • 条件1，如果 $e(k)$ 是白噪声序列（四阶矩存在）；
      • 条件2，待辨识对象的特征值在单位圆内；
      • 条件3，输入序列 $u(k)$ 与 $e(k)$ 噪声序列相互独立；
      • 条件4，输入信号必须是 $nb$ 阶持续激励信号；

      则最小二乘参数估计是一致收敛的，即
      $$\underset{N\to \infty }{\lim }{\hat{\theta }}_{LS}={\theta }_0,W.P.1$$

四、非一致性的举例说明

• 考虑如下模型

$$y(k)=-ay(k-1)+bu(k-1)+e(k)$$

​ 其中， $|a|<1$ ， $e(k)$ 和 $u(k)$ 是统计独立的， $e(k)$ 是均值为零的各台遍历平稳噪声序列，假设 $R_e(l)=0,\forall l\ge 2$。
$${\hat{\theta }}_{LS}\stackrel{N\to \infty ,W.P.1}{\to }\left[\begin{array}{l}{a}_0-\frac{R_u(0)R_e(1)}{\Delta }\\ b_0-\frac{R_{u}y(0)R_e(1)}{\Delta }\end{array}\right]$$

附录1、关于噪声方差估计的证明

附录2、最小二乘几何解释的证明过程

07-最小二乘参数估计的缺陷以及改进算法（最小二乘类5）

一、最小二乘算法的缺陷

• 噪声序列一般是相关的，所以“一致性”中的假设“噪声项不相干”一般并不成立。

二、增广矩阵最小二乘方法

1. Figure(ARMAX模型):

   

2. 矩阵形式系统模型：

$$A(z^{-1})y(k)=B(z^{-1})u(k)+C(z^{-1})\epsilon (k)$$

$$\begin{array}{l}A(z^{-1})=1+a_1{z}^{-1}+\cdots +a_n{z}^{-n}\\ B(z^{-1})=b_1{z}^{-1}+\cdots +b_n{z}^{-n}\\ C(z^{-1})=1+c_1{z}^{-1}+\cdots +c_n{z}^{-n}\end{array}$$

2. 分量形式系统模型：

$$\begin{array}{rl}y(k)=& -a_{1}y(k-1)-a_{2}y(k-2)-\cdots -a_{n}y(k-n)\\ & +b_{1}u(k-1)+b_{2}u(k-2)+\cdots +b_{n}u(k-n)\\ & +\epsilon (k)+c_1\epsilon (k-1)+\cdots +c_n\epsilon (k-n)\end{array}$$

3. 向量形式系统模型：

$$\begin{gathered} y(k)={\phi }_e(k{)}^T{\theta }_e+\epsilon (k) \\ {\phi }_e(k{)}^T=\left[-y(k-1),\cdots ,-y(k-n),u(k-1),\cdots ,u(k-n),\epsilon (k-1),\cdots ,\epsilon (k-n)\right] \\ {\theta }_e^T=\left[a_1,\cdots ,a_n,b_1,\cdots ,n_n,c_1,\cdots ,c_n\right] \\ \epsilon (k+1)=y(k+1)-{\phi }_e(k+1{)}^T{\hat{\theta }}_e(k) \\ \hat{\epsilon }(1)=\hat{\epsilon }(2)=\cdots =\hat{\epsilon }(n)=0 \end{gathered}$$

4. 例1：

   1. 系统方程：$$y(k)-1.5y(k-1)+0.7y(k-2)=u(k-1)+0.5u(k-2)+e(k)+0.2e(k-1)$$

   2. MATLAB程序：C07P13_example1.m

   3. 先估计噪声：
      $$\begin{gathered} \hat{\epsilon }(1)=\hat{\epsilon }(2)=\cdots =\hat{\epsilon }(n)=0 \\ fo{r}_{(i=n:end)}：\hat{\epsilon }(i+1)=y(i+1)-{\phi }_e(k+1{)}^T{\hat{\theta }}_e(k) \\ 对n=2的情况，\hat{\epsilon }(1)=\hat{\epsilon }(2)=0，\hat{\epsilon }(3)之后不为0 \\ {\phi }_e(i)=[-y(i-1),-y(i-2),u(i-1),u(i-2),\epsilon (i-1),\epsilon (i-2)] \end{gathered}$$
      再计算递推最小二乘的 $\phi$ 向量：
      $$\begin{gathered} {\phi }^T(i+1)=[-y(n+i),-y(n+i-1),\dots ,-y(i+1),\dots \\ u(n+i),u(n+i-1),\dots ,u(i+1),\dots \\ \epsilon (n+i),\epsilon (n+i-1),\dots ,\epsilon (i+1)] \\ 对n=2的情况，fo{r}_{(i=0:end)}：{\phi }^T(i+1)=[-y(i+2),-y(i+1),u(i+2),u(i+1),\epsilon (i+2),\epsilon (i+1)] \end{gathered}$$
      最后进行递推最小二乘：
      $$\begin{array}{l}fo{r}_{(i=0:end-n-1)}：\\ K(i+1)=P(i)\phi (i+1){\left[1+{\phi }^T(i+1)P(N)\phi (i+1)\right]}^{-1}\\ \hat{\theta }(i+1)=\hat{\theta }(N)+K(i+1)\left[y(n+i+1)-{\phi }^T(i+1)\hat{\theta }(i)\right]\\ P(i+1)=P(i)-K(i+1)K^T(i+1){\left[1+{\phi }^T(i+1)P(i)\phi (i+1)\right]}^{-1}\end{array}$$

   4. 问题：

      1. 一个向前推，一个向后推，差着步数，怎么解决？
      2. 噪声模型只有一阶，与输入和模型阶数不一致怎么办？

5. 注释：相比普通最小二乘法，同时考虑了噪声的模型，阶次扩大了，因此称之为增广矩阵法。

三、广义最小二乘方法的思想以及推导

1. Figure:

   

2. 矩阵形式系统模型：

$$\begin{array}{rl}y(k)& =\frac{B(z^{-1})}{A(z^{-1})}u(k)+\frac{G{(}^{-1})}{F(z^{-1})}\epsilon (K)\\ A(z^{-1})y(k)& =B(z^{-1})u(k)+\frac{A(z^{-1})G{(}^{-1})}{F(z^{-1})}\epsilon (K)\\ & \underset{e(k)=H(z^{-1})\epsilon (k)}{\overset{H(z^{-1})=\frac{A(z^{-1})G{(}^{-1})}{F(z^{-1})}}{\to }}\\ H^{-1}(z^{-1})[A(z^{-1})y(k)]& =H^{-1}(z^{-1})[B(z^{-1})u(k)]+\epsilon (K)\\ A(z^{-1})[H^{-1}(z^{-1})y(k)]& =B(z^{-1})[H^{-1}(z^{-1})u(k)]+\epsilon (K)\\ & ⟶\\ A(z^{-1})y_f(k)& =B(z^{-1})u_f(k)+\epsilon (k)\\ & \{\begin{array}{l}{y}_f(k)=H(z^{-1})y(k)\\ u_f(k)=H(z^{-1})u(k)\end{array}\end{array}$$

2. 广义最小二乘法的思想：
   1. 对噪声项用“线性滤波器”尽量好的白色化处理；
   2. 利用滤波之后的数据重新进行参数估计。
3. 对于简化的情况：

$$H(z^{-1})=\frac{1}{1+C_1{z}^{-1}+\dots +C_Y{z}^{-Y}}=\frac{1}{C(z^{-1})}$$

​ 则：
$$\{\begin{array}{l}{y}_f(k)=C(z^{-1})y(k)\\ u_f(k)=C(z^{-1})u(k)\\ \epsilon (k)=C(z^{-1})e(k)\end{array}$$

4. 例2

   1. 系统函数：
      $$\begin{gathered} y(k)-1.5y(k-1)+0.7y(k-2)=u(k-1)+0.5u(k-2)+e(k) \\ e(k)-1.2e(k-1)+0.5e(k-2)=\epsilon (k) \\ e(k)=\frac{1}{1-1.2z^{-1}+0.5z^{-2}}\epsilon (k) \\ H(z^{-1})=\frac{1}{1-1.2z^{-1}+0.5z^{-2}} \\ C(z^{-1})=1-1.2z^{-1}+0.5z^{-2} \end{gathered}$$

   2. MATLAB程序：C07P_example2.m

5. 问题：如何估计噪声模型 $C(Z)$ ，能否用迭代的方式？

四、广义最小二乘离线算法流程

$$\begin{gathered} 令C(z^{-1})=1⟶计算\{\begin{array}{l}{y}_f(k)=C(z^{-1})y(k)\\ u_f(k)=C(z^{-1})u(k)\end{array}⟶\begin{array}{c}依据A(z^{-1})y_f(k)=B(z^{-1})u_f(k)\\ 利用最小二乘估计A(z^{-1}),B(z^{-1})得到\hat{A}(z^{-1}),\hat{B}(z^{-1})\end{array} \\ ⟶计算\{\begin{array}{l}e(k)=\hat{A}(z^{-1})y(k)-\hat{B}(z^{-1})u(k)\\ e^{\ast }(k)=\hat{A}(z^{-1})y_f(k)-\hat{B}{(}^{-1})u_f(k)\end{array}⟶\begin{array}{c}依据C(z^{-1})e(k)=\epsilon (k)\\ 利用最小二乘估计\hat{C}(z^{-1})\end{array}⟶\begin{array}{c}如果收敛并满足精度要求则停止运算；\\ 否则转至第二步。\\ 所谓收敛：{\lim }_{j\to \infty }{\hat{C}}_j^{\ast }(z^{-1})=1\end{array} \end{gathered}$$

五、广义最小二乘程序部分代码

六、广义最小二乘在线算法流程

附录1、广义最小二乘递推推导过程