应用性学习（一）：控制系统辨识基础知识

控制系统辨识

一、M序列

参数辨识的参考输入

全称Maximal Length Sequence

用于产生伪随机二进制序列，近似白噪声的统计特性

常用的反馈通道组

M序列参数选择

• 幅值$a$ ：根据被辨识系统的线性范围和允许的信噪比来确定。一方面，$a$必须对系统进行充分的激励;另一方面，$a$ 又不能太大，否则将引起系统非线性和较大的干扰噪声。

• 移位脉冲周期$\Delta$：根据被识系统的频带或截止频率这些重要的工作频率作为参考进行4的选择,设最高的参考工作频率为$w_{max}$为使 M 序列的有效频带能够覆盖被辨识系统的重要工作频率区，应满足
  $$\frac{2\pi }{3\Delta }>{\omega }_{\max }\text{或 }\Delta <\frac{2\pi }{3{\omega }_{\max }}$$

• 循环周期$N$：为使系统脉冲响应在 M 序列一个周期时间**$N\Delta$内近似衰减到零，应满足$T=N\Delta >t_s$**，其中t为系统的调节时间。如果循环周期取得太大，在实际应用中辨识用时间较长,一般取为 $N\Delta =(1.2-1.5)t_s$。

Example

$a=0.5mA,\Delta =0.6s,N=31$，使用四级移位寄存器实现M序列

n=4;
N=2^n-1;
a=0.5;
Delta=0.6;
A1=1;A2=1;A3=1;A4=0; %设置初始值
for i=1:3*N
    X1=xor(A3,A4);  %对第三级和第四级移位寄存器输出进行模二相加
    X2=A1;
    X3=A2;
    X4=A3;
    OUT(i)=A4;     %移位寄存器最后一级输出
    t(i)=Delta*i;
    if OUT(i)>0.5   
        u(i)=-a;   %确定电平幅值
    else u(i)=a;
    end
    A1=X1;A2=X2;A3=X3;A4=X4;    
end

仿真结果

二、最小二乘基础

2.1 一般最小二乘

• $u(k)$ 系统输入
• $y(k)$ 系统输出
• $v(k)$ 系统噪声
• $z(k)$ 实际可测输出

被辨识模型为 $G(z)$
$$G(z)=\frac{y(z)}{u(z)}=\frac{b_1{z}^{-1}+b_2{z}^{-2}+\cdot \cdot \cdot +b_n{z}^{-n}}{1+a_1{z}^{-1}+a_2{z}^{-2}+\cdot \cdot \cdot +a_n{z}^{-n}}$$
对应差分方程为
$$\begin{array}{rl}y(k)=& -\sum _{i=1}^n{a}_{i}y(k-i)+\sum _{i=1}^n{b}_{i}u(k-i)\end{array}$$
改写成
$$z(k)=-\sum _{i=1}^n{a}_{i}y(k-i)+\sum _{i=1}^n{b}_{i}u(k-i)+v(k)$$
定义待辨识参数矩阵和系统输入输出矩阵
$$\begin{gathered} h(k)=\left[-y(k-1),-y(k-2),\cdots ,-y(k-n),u(k-1),u(k-2),\cdots ,u(k-n)\right] \\ \theta ={\left[a_1,a_2,\cdots ,a_n,b_1,b_2,\cdots ,b_n\right]}^{⊺} \end{gathered}$$
则被辨识模型式
$$\mathbf{z}(k)=\mathbf{h}(k)\mathbf{\theta }+\mathbf{v}(k)$$
写成矩阵形式，如下
$${\mathbf{Z}}_m={\mathbf{H}}_m\mathbf{\theta }+{\mathbf{V}}_m$$
其中
$$\begin{array}{rl}{\mathbf{Z}}_m& =\left[\begin{array}{c}z(1)\\ z(2)\\ ⋮\\ z(m)\end{array}\right]\\ {\mathbf{H}}_m& =\left[\begin{array}{c}h(1)\\ h(2)\\ ⋮\\ h(m)\end{array}\right]=\left[\begin{array}{cccccc}-y(0)& \cdots & -y(1-n)& u(0)& \cdots & u(1-n)\\ -y(1)& \cdots & -y(2-n)& u(1)& \cdots & u(2-n)\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ -y(m-1)& \cdots & -y(m-n)& u(m-1)& \cdots & u(m-n)\end{array}\right]\\ \mathbf{\theta }& ={\left[\begin{array}{c}{a}_1,a_2,\cdots ,a_n,b_1,b_2,\cdots ,b_n\end{array}\right]}^{⊺},{\mathbf{V}}_m={\left[\begin{array}{cccc}v(1)& v(2)& \cdots & v(m)\end{array}\right]}^{⊺}\end{array}$$
估计的损失函数
$$J(\hat{\mathbf{\theta }})=({\mathbf{Z}}_m-{\mathbf{H}}_m\hat{\mathbf{\theta }}{)}^{\mathrm{T}}({\mathbf{Z}}_m-{\mathbf{H}}_m\hat{\mathbf{\theta }})$$
为使其损失函数最小，根据极值定理
$${\frac{\partial J}{\partial \mathbf{\theta }}|}_{\mathbf{\theta }=\hat{\mathbf{\theta }}}=-2{\mathbf{H}}_m^{\mathrm{T}}({\mathbf{Z}}_m-{\mathbf{H}}_m\hat{\mathbf{\theta }})=0$$
则
$${\mathbf{H}}_m^{\mathrm{T}}{\mathbf{H}}_m\hat{\mathbf{\theta }}={\mathbf{H}}_m^{\mathrm{T}}{\mathbf{Z}}_m$$
若${\mathbf{H}}_m^{\mathrm{T}}{\mathbf{H}}_m$ 满秩，则最小二乘估计为
$$\hat{\mathbf{\theta }}=({\mathbf{H}}_m^{\mathrm{T}}{\mathbf{H}}_m{)}^{-1}{\mathbf{H}}_m^{\mathrm{T}}{\mathbf{Z}}_m$$

2.2 加权最小二乘

其损失函数为
$$J(\hat{\mathbf{\theta }})=({\mathbf{Z}}_m-{\mathbf{H}}_m\hat{\mathbf{\theta }}{)}^{\mathrm{T}}{\mathbf{W}}_m({\mathbf{Z}}_m-{\mathbf{H}}_m\hat{\mathbf{\theta }})$$
其中加权矩阵为
$${\mathbf{W}}_m=\mathrm{diag}[w(1),w(2),\cdots ,w(m)]$$
加权最小二乘估计为
$$\hat{\mathbf{\theta }}=({\mathbf{H}}_m^{\mathrm{T}}{\mathbf{W}}_m{\mathbf{H}}_m{)}^{-1}{\mathbf{H}}_m^{\mathrm{T}}{\mathbf{W}}_m{\mathbf{Z}}_m$$
若${\mathbf{W}}_m={\mathbf{R}}^{-1}$ ，R为噪声协方差矩阵，则为马尔科夫估计，是所有加权最小二乘估计的最优者。

2.3 递推最小二乘

不用每次对整体的数据进行最小二乘的估计运算，利用上次计算的值从和下一次到来的数据进行融合计算，增加计算效率，修正旧得计算值从而达到新的计算效果。

原m次数据的加权最小二乘估计为
$${\hat{\mathbf{\theta }}}_m=({\mathbf{H}}_m^{\mathrm{T}}{\mathbf{W}}_m{\mathbf{H}}_m{)}^{-1}{\mathbf{H}}_m^{\mathrm{T}}{\mathbf{W}}_m{\mathbf{Z}}_m$$
获得新输入输出数据
$$z(m+1)=h(m+1)\theta +v(m+1)$$
加权最小二乘估计的递推算法为（一般的同理）
$$\begin{array}{rl}{\hat{\mathbf{\theta }}}_{m+1}& ={\hat{\mathbf{\theta }}}_m+{\mathbf{K}}_{m+1}\left[z(m+1)-\mathbf{h}(m+1){\hat{\mathbf{\theta }}}_m\right]\\ {\mathbf{P}}_{m+1}& ={\mathbf{P}}_m-{\mathbf{P}}_m{\mathbf{h}}^T(m+1)[w^{-1}(m+1)+\mathbf{h}(m+1){\mathbf{P}}_m{\mathbf{h}}^T(m+1){]}^{-1}\mathbf{h}(m+1){\mathbf{P}}_m\\ & ={\mathbf{P}}_m[\mathbf{I}-{\mathbf{K}}_{m+1}\mathbf{h}(m+1)]\\ {\mathbf{K}}_{m+1}& ={\mathbf{P}}_m{\mathbf{h}}^T(m+1)[w^{-1}(m+1)+\mathbf{h}(m+1){\mathbf{P}}_m{\mathbf{h}}^T(m+1){]}^{-1}\end{array}$$
其中${\mathbf{P}}_m=({\mathbf{H}}_m^{\mathrm{T}}{\mathbf{W}}_m{\mathbf{H}}_m{)}^{-1}$ ，$w(m+1)$ 为新输入输出数据的噪声方差

精度相关的停机准则
$$\underset{\forall i}{\max }|\frac{{\hat{\theta }}_i(m+1)-{\hat{\theta }}_i(m)}{{\hat{\theta }}_i(m)}|<\epsilon$$

对于没有遗忘因子的最小二乘递推算法，若对于时变系统和观测数据累计的情况下会出现误差越来越大的情况，甚至会有数据饱和的情况。

2.4 时变最小二乘

（1）矩形窗法

矩形窗法的特点是在时刻 k 的估计只依靠近k时刻的有限个数据，此前的数据全部被剔除掉，长度为 m 的矩形窗，每个时刻增加一个新数据，就要剔除一个旧数据，保持数据窗内的数据数目始终为 m。

（2）指数窗法

和上述差不多，但是带遗忘因子的矩形窗法

2.5 递推阻尼最小二乘

随着协方差矩阵的减小，参数易产生爆发现象。为防止辨识参数产生爆发现象，在最小二乘参数辨识的目标函数中增加一个对参数变化量的阻尼项，来增加算法的稳定性。

递推阻尼最小二乘法的目标函数为
$$J=\left[\begin{array}{c}(Y_m-{\mathbf{H}}_m{\hat{\mathbf{\theta }}}_m{)}^{\mathrm{T}}{\mathbf{W}}_m({\mathbf{Y}}_m-{\mathbf{H}}_m{\hat{\mathbf{\theta }}}_m)\end{array}\right]+\mu {\left[\begin{array}{c}{\hat{\mathbf{\theta }}}_m-{\hat{\mathbf{\theta }}}_{m-1}\end{array}\right]}^{\mathrm{T}}\left[\begin{array}{c}{\hat{\mathbf{\theta }}}_m-{\hat{\mathbf{\theta }}}_{m-1}\end{array}\right]$$
公式为
$$\begin{array}{rl}{\hat{\mathbf{\theta }}}_m=& {\hat{\mathbf{\theta }}}_{m-1}+{\lambda }_f\mu {\mathbf{P}}_m[{\hat{\mathbf{\theta }}}_{m-1}-{\hat{\mathbf{\theta }}}_{m-2}]+{\mathbf{P}}_m{\mathbf{h}}_m^{\mathrm{T}}[{\mathbf{y}}_m-{\mathbf{h}}_m{\hat{\mathbf{\theta }}}_{m-1}]\\ {\mathbf{P}}_m=& [\mu \mathbf{I}+\lambda {\mathbf{H}}_{m-1}^{\mathrm{T}}{\mathbf{W}}_{m-1}{\mathbf{H}}_{m-1}+{\mathbf{h}}_m^{\mathrm{T}}{\mathbf{h}}_m{]}^{-1}\end{array}$$

2.6 增广最小二乘

当噪声均值为0时，最小二乘参数估计算法为无偏估计；当噪声的均值不为 0时，最小乘参数估计算法为有偏估计。为了解决最小二乘参数估计的有偏性，将噪声模型的辨识同时考虑进去，因此称为增广最小二乘法。

辨识对象为
$$z(k)=-\sum _{i=1}^n{a}_{i}y(k-i)+\sum _{i=1}^n{b}_{i}u(k-i)+\sum _{i=1}^n{c}_{i}v(k-i)+v(k)$$
对辨识对象进行扩充（考虑对有色噪声的辨识）
$$\begin{array}{rl}\mathbf{h}(k)=& [-y(k-1),-y(k-2),\cdots ,-y(k-n),u(k-1),u(k-2),\cdots ,\\ & u(k-n),\hat{v}(k-1),\cdots ,\hat{v}(k-n)]\end{array}$$

$$\mathbf{\theta }=[a_1,a_2,\cdot \cdot \cdot ,a_n,b_1,b_2,\cdot \cdot \cdot ,b_n,c_1,\cdot \cdot \cdot ,c_n{]}^T$$

则简化方程为
$$z(k)=\mathbf{h}(k)\mathbf{\theta }+\upsilon (k)$$
此种增广递推算法为
$$\begin{array}{c}{\hat{\mathbf{\theta }}}_{m+1}={\hat{\mathbf{\theta }}}_m+{\mathbf{K}}_{m+1}\left[z(m+1)-\mathbf{h}(m+1){\hat{\mathbf{\theta }}}_m\right]\\ {\mathbf{P}}_{m+1}={\mathbf{P}}_m-{\mathbf{P}}_m{\mathbf{h}}^{\mathrm{T}}(m+1){\left[w^{-1}(m+1)+\mathbf{h}(m+1){\mathbf{P}}_m{\mathbf{h}}^{\mathrm{T}}(m+1)\right]}^{-1}\mathbf{h}(m+1){\mathbf{P}}_m\\ {\mathbf{K}}_{m+1}={\mathbf{P}}_m{\mathbf{h}}^{\mathrm{T}}(m+1){\left[w^{-1}(m+1)+\mathbf{h}(m+1){\mathbf{P}}_m{\mathbf{h}}^{\mathrm{T}}(m+1)\right]}^{-1}\end{array}$$

2.7 多变量系统最小二乘辨识

多变量系统表示方法
$$\begin{array}{rl}\mathbf{Y}(k)+{\mathbf{A}}_1\mathbf{Y}(k-1)+\cdots +{\mathbf{A}}_n\mathbf{Y}(k-n)& ={\mathbf{B}}_0\mathbf{U}(k)+{\mathbf{B}}_1\mathbf{U}(k-1)+\cdot \cdot \cdot +\\ & {\mathbf{B}}_n\mathbf{U}(k-n)+\mathbf{V}(k)\end{array}$$

• $\mathbf{Y}(k)$ 为m维输出
• $\mathbf{U}(k)$ 为r维输入
• $\mathbf{V}(k)$ 为m维噪声
• ${\mathbf{A}}_i$ 为待辨识m×n维矩阵
• ${\mathbf{B}}_i$ 为待辨识m×r维矩阵

即
$$\mathbf{A}(z^{-1})\mathbf{Y}(k)=\mathbf{B}(z^{-1})\mathbf{U}(k)+\mathbf{V}(k)$$
其中
$$\begin{gathered} \mathbf{A}(z^{-1})=\mathbf{I}+{\mathbf{A}}_1{z}^{-1}+\cdots +{\mathbf{A}}_n{z}^{-n}=\mathbf{I}+\sum _{i=1}^n{\mathbf{A}}_i{z}^{-i} \\ \mathbf{B}(z^{-1})={\mathbf{B}}_0+{\mathbf{B}}_1{z}^{-1}+\cdot \cdot \cdot +{\mathbf{B}}_n{z}^{-n}=\sum _{i=0}^n{\mathbf{B}}_i{z}^{-i} \end{gathered}$$
待辨识总数目：$n{m}^2+(n+1)mr$

将待辨识的矩阵化为向量形式，如下
$$\begin{array}{c}{\mathbf{\theta }}_j={\left[\begin{array}{cccccccccccccc}{a}_{11}& \cdots & a_{1m}^j& \cdots & a_{n1}^j& \cdots & a_{nm}^j& b_{01}^j& \cdots & b_{0r}& \cdots & b_{n1}^j& \cdots & b_{nr}^j\end{array}\right]}^T\\ {\mathbf{H}}_j=[\begin{array}{cccccc}-{\mathbf{Y}}_j^{\mathrm{T}}(0)& \cdots & -{\mathbf{Y}}_j^{\mathrm{T}}(1-n)& {\mathbf{U}}^{\mathrm{T}}(1)& \cdots & {\mathbf{U}}^{\mathrm{T}}(1-n)\\ -{\mathbf{Y}}_j^{\mathrm{T}}(1)& \cdots & -{\mathbf{Y}}_j^{\mathrm{T}}(2-n)& {\mathbf{U}}^{\mathrm{T}}(2)& \cdots & {\mathbf{U}}^{\mathrm{T}}(2-n)\\ ⋮& ⋮& ⋮& ⋮& ⋮& \\ -{\mathbf{Y}}_j^{\mathrm{T}}(N-1)& \cdots & -{\mathbf{Y}}_j^{\mathrm{T}}(N-n)& {\mathbf{U}}^{\mathrm{T}}(N)& \cdots & {\mathbf{U}}^{\mathrm{T}}(N-n)\end{array}]\end{array}$$
于是简化为
$${\mathbf{Y}}_j={\mathbf{H}}_j{\mathbf{\theta }}_j+{\mathbf{V}}_j$$
（1）离线辨识（直接辨识）
$${\hat{\mathbf{\theta }}}_j=({\mathbf{H}}_j^{\mathrm{T}}{\mathbf{H}}_j{)}^{-1}{\mathbf{H}}_j^{\mathrm{T}}{\mathbf{Y}}_j$$
（2）在线辨识（递推辨识）

前N次结果之后得到新的观测值 ${\mathbf{y}}_j(N+1)$ 和 $\mathbf{U}(N+1)$
$${\hat{\mathbf{\theta }}}_{jN}=({\mathbf{H}}_{jN}^{\mathrm{T}}{\mathbf{H}}_{jN}{)}^{-1}{\mathbf{H}}_{jN}^{\mathrm{T}}{\mathbf{Y}}_{jN}$$

$$\begin{array}{c}{\hat{\mathbf{\theta }}}_{j(N+1)}={\hat{\mathbf{\theta }}}_{jN}+{\mathbf{K}}_{j(N+1)}{\text{[y}}_j\left(N+1\right)-{\mathbf{h}}_{j(N+1)}^{\mathrm{T}}{\hat{\mathbf{\theta }}}_{jN}]\\ {\mathbf{K}}_{j(N+1)}={\mathbf{P}}_{jN}{\mathbf{h}}_{j(N+1)}[1+{\mathbf{h}}_{j(N+1)}^{\top }{\mathbf{P}}_{jN}{\mathbf{h}}_{j(N+1)}{]}^{-1}\\ {\mathbf{P}}_{j(N+1)}={\mathbf{P}}_{jN}-{\mathbf{P}}_{jN}{\mathbf{h}}_{j(N+1)}[1+{\mathbf{h}}_{j(N+1)}^{\top }{\mathbf{P}}_{jN}{\mathbf{h}}_{j(N+1)}{]}^{-1}{\mathbf{h}}_{j(N+1)}{\mathbf{P}}_{jN}\\ {\mathbf{P}}_{jN}=({\mathbf{H}}_{jN}^{\mathrm{T}}{\mathbf{H}}_{jN}{)}^{-1}\end{array}$$