精华内容
下载资源
问答
  • 该代码为龚建伟《无人驾驶车辆模型预测控制》中第3章3.3.3的例子。 1 公式推导 这一部分的推导过程,比我上一篇文章的推导较简单一些,主要表现在三个方面:第1是控制量由前轮转角δ变成了角速度ω,第2是没有对状态...

    无人驾驶轨迹跟踪仿真——线性时变模型预测代码详解

    对照推导的公式,对代码进行一一详细注解,方便学习代码的同学。该代码为龚建伟《无人驾驶车辆模型预测控制》中第3章3.3.3的例子。

    1 公式推导

    这一部分的推导过程,比我上一篇文章的推导较简单一些,主要表现在三个方面:第1是控制量由前轮转角δ变成了角速度ω,第2是没有对状态量和控制量进行组合,第3是代价函数没有添加松弛变量。
    具体的推导过程,我放上了我的手稿。
    1
    推导公式中的矩阵形式,以及最终的二次规划形式在后边的代码中都有所体现,我会在代码部分作出解释。

    2 代码部分

    代码共分为了5部分。

    2.1 参考轨迹生成

    参考轨迹为Y=2的一条直线

    N=100;                                   %参考轨迹点的数量
    T=0.05;                                  %采样周期
    Xout=zeros(N,3);                         %生成N行3列的0矩阵,为参考点矩阵
    Tout=zeros(N,1);                         %生成N行1列的0矩阵
    for k=1:1:N
        Xout(k,1)=k*T;                       %1列表示x值
        Xout(k,2)=2;                         %2列表示y值
        Xout(k,3)=0;                         %3列表示φ值
        Tout(k,1)=(k-1)*T;                   %为时刻矩阵
    end
    

    2.2 常量参数赋值

    Nx=3;                                     %状态量个数
    Nu =2;                                    %控制量个数
    Tsim =20;                                 %文中注解仿真时间,我认为是预测时域Np
    X0 = [0 0 pi/3];                          %初始位置的状态
    [Nr,Nc] = size(Xout);                     %取Xout矩阵的行数和列数分别赋值给Nr和Nc
    L = 1;                                    %车辆轴距
    Rr = 1;                                   %轮胎直径
    w = 1;                                    %车轮转速
    vd1 = Rr*w;                               %参考系统的纵向速度
    vd2 = 0;                                  %参考系统的角速度
    

    2.3 定义矩阵

    x_real=zeros(Nr,Nc);                       %状态量矩阵   
    x_piao=zeros(Nr,Nc);                       %状态量误差矩阵
    u_real=zeros(Nr,2);                        %控制量矩阵 
    u_piao=zeros(Nr,2);                        %控制量误差矩阵
    x_real(1,:)=X0;                            %把初始状态赋值给状态量第一行
    x_piao(1,:)=x_real(1,:)-Xout(1,:);         %计算第一个状态量误差
    X_PIAO=zeros(Nr,Nx*Tsim);
    XXX=zeros(Nr,Nx*Tsim);                     %用于保持每个时刻预测的所有状态值
    q=[1 0 0;0 1 0;0 0 0.5];                   %状态量的权重系数
    Q_cell=cell(Tsim,Tsim);
    for i=1:1:Tsim
        for j=1:1:Tsim
            if i==j
                Q_cell{i,j}=q;
            else 
                Q_cell{i,j}=zeros(Nx,Nx);
            end 
        end
    end
    Q=cell2mat(Q_cell);                         %将多个矩阵合并成一个矩阵,得到权重系数矩阵Q
    R=0.1*eye(Nu*Tsim,Nu*Tsim);                 %权重系数矩阵R
    

    2.4 MPC的预测和优化

    for i=1:1:Nr
        t_d =Xout(i,3);                         %取航向角φ
        a=[1    0   -vd1*sin(t_d)*T;
           0    1   vd1*cos(t_d)*T;
           0    0   1;];                        %状态量系数矩阵
        b=[cos(t_d)*T   0;
           sin(t_d)*T   0;
           0            T;];                    %控制量系数矩阵   
        A_cell=cell(Tsim,1);
        B_cell=cell(Tsim,Tsim);
         for j=1:1:Tsim
            A_cell{j,1}=a^j;                    %A_cell为 ψ 矩阵的分量
            for k=1:1:Tsim
               if k<=j
                    B_cell{j,k}=(a^(j-k))*b;    %B_cell为 θ 矩阵的分量
               else
                    B_cell{j,k}=zeros(Nx,Nu);
               end
            end
        end
        A=cell2mat(A_cell);                     %合并A_cell生成 ψ 矩阵
        B=cell2mat(B_cell);                     %合并B_cell生成 θ 矩阵
        
        H=2*(B'*Q*B+R);                         %二次规划中的H
        f=(2*x_piao(i,:)*A'*Q*B)';              %二次规划中的f
       
        %% 约束条件
        A_cons=[];                              %不等式约束的系数矩阵
        b_cons=[];                              %不等式约束的值
        lb=[-1;-1];                             %自变量约束下界
        ub=[1;1];                               %自变量约束上界
        tic                                     %用来保存当前时间
        [X,fval(i,1),exitflag(i,1),output(i,1)]=quadprog(H,f,A_cons,b_cons,[],[],lb,ub);   %%二次规划函数求解
        toc                                     %用来记录程序完成时间
        X_PIAO(i,:)=(A*x_piao(i,:)'+B*X)';      %输出矩阵Y
          
          %% 输出每个时刻预测的所有的状态量
          %% i表示时刻k,预测时域值为20,则当时刻k=80开始后,预测到的时刻要超100了,所以之后都取100时刻的值。
        if i+j<Nr
             for j=1:1:Tsim
                 XXX(i,1+3*(j-1))=X_PIAO(i,1+3*(j-1))+Xout(i+j,1);
                 XXX(i,2+3*(j-1))=X_PIAO(i,2+3*(j-1))+Xout(i+j,2);
                 XXX(i,3+3*(j-1))=X_PIAO(i,3+3*(j-1))+Xout(i+j,3);
             end
            
        else
             for j=1:1:Tsim
                 XXX(i,1+3*(j-1))=X_PIAO(i,1+3*(j-1))+Xout(Nr,1);
                 XXX(i,2+3*(j-1))=X_PIAO(i,2+3*(j-1))+Xout(Nr,2);
                 XXX(i,3+3*(j-1))=X_PIAO(i,3+3*(j-1))+Xout(Nr,3);
             end
        end
        u_piao(i,1)=X(1,1);                     %将二次规划的第一个解的x值赋值给控制量矩阵当前时刻x值
        u_piao(i,2)=X(2,1);                     %将二次规划的第一个解的y值赋值给控制量矩阵当前时刻y值
        Tvec=[0:0.05:4];
        X00=x_real(i,:);                        %将当前时刻的状态量赋值给X00
        vd11=vd1+u_piao(i,1);                   %当前时刻的纵向速度
        vd22=vd2+u_piao(i,2);                   %当前时刻的角速度
        
        %% 基于控制量计算下一时刻的状态量
        XOUT=dsolve('Dx-vd11*cos(z)=0','Dy-vd11*sin(z)=0','Dz-vd22=0','x(0)=X00(1)','y(0)=X00(2)','z(0)=X00(3)');
         t=T; 
         %% 记录下一时刻的状态量
         x_real(i+1,1)=eval(XOUT.x);
         x_real(i+1,2)=eval(XOUT.y);
         x_real(i+1,3)=eval(XOUT.z);
         if(i<Nr)
             x_piao(i+1,:)=x_real(i+1,:)-Xout(i+1,:); %计算下一时刻的误差
         end
        u_real(i,1)=vd1+u_piao(i,1);  
        u_real(i,2)=vd2+u_piao(i,2);              %当前时刻的控制量矩阵
        
        figure(1);
        plot(Xout(1:Nr,1),Xout(1:Nr,2));          %作图参考轨迹点x-y
        hold on;
        plot(x_real(i,1),x_real(i,2),'r*');       %作图状态量x-y
        title('跟踪结果对比');
        xlabel('横向位置X');
        axis([-1 5 -1 3]);
        ylabel('纵向位置Y');
        hold on;
        for k=1:1:Tsim
             X(i,k+1)=XXX(i,1+3*(k-1));           %每个时刻的预测值x
             Y(i,k+1)=XXX(i,2+3*(k-1));           %每个时刻的预测值y
        end
        X(i,1)=x_real(i,1);   
        Y(i,1)=x_real(i,2);
        plot(X(i,:),Y(i,:),'y')                   %所有时刻的预测值x-y
        hold on;
        
    end
    

    注:代码每一个x_piao的分量在形式上是与推导中是转置的关系,所以在二次规划表达式中 f 的表达式可以看出来。如图1所示:
    在这里插入图片描述

    2.5 绘图

    figure(2)
    subplot(3,1,1);
    plot(Tout(1:Nr),Xout(1:Nr,1),'k--');                      %时刻-参考点 x 值
    hold on;
    plot(Tout(1:Nr),x_real(1:Nr,1),'k');                      %时刻-状态量 x 值
    %grid on;
    %title('状态量-横向坐标X对比');
    xlabel('采样时间T');
    ylabel('横向位置X')
    subplot(3,1,2);
    plot(Tout(1:Nr),Xout(1:Nr,2),'k--');                      %时刻-参考点 y 值
    hold on;
    plot(Tout(1:Nr),x_real(1:Nr,2),'k');                      %时刻-状态量 y 值
    %grid on;
    %title('状态量-横向坐标Y对比');
    xlabel('采样时间T');
    ylabel('纵向位置Y')
    subplot(3,1,3);
    plot(Tout(1:Nr),Xout(1:Nr,3),'k--');                      %时刻-参考点 φ 值
    hold on;
    plot(Tout(1:Nr),x_real(1:Nr,3),'k');                      %时刻-状态量 φ 值
    %grid on;
    hold on;
    %title('状态量-\theta对比');
    xlabel('采样时间T');
    ylabel('\theta')
    
    figure(3)
    subplot(2,1,1);
    plot(Tout(1:Nr),u_real(1:Nr,1),'k');                       %时刻-控制量 v 值
    %grid on;
    %title('控制量-纵向速度v对比');
    xlabel('采样时间T');
    ylabel('纵向速度')
    subplot(2,1,2)
    plot(Tout(1:Nr),u_real(1:Nr,2),'k');                       %时刻-控制量 w 值
    %grid on;
    %title('控制量-角加速度对比');
    xlabel('采样时间T');
    ylabel('角加速度')
    
    figure(4)
    subplot(3,1,1);
    plot(Tout(1:Nr),x_piao(1:Nr,1),'k');                       %时刻-控制量误差 x 值
    %grid on;
    xlabel('采样时间T');
    ylabel('e(x)');
    subplot(3,1,2);
    plot(Tout(1:Nr),x_piao(1:Nr,2),'k');                       %时刻-控制量误差 y 值
    %grid on;
    xlabel('采样时间T');
    ylabel('e(y)');
    subplot(3,1,3);
    plot(Tout(1:Nr),x_piao(1:Nr,3),'k');                       %时刻-控制量误差 φ 值
    %grid on;
    xlabel('采样时间T');
    ylabel('e(\theta)');
    

    所有图如下所示:
    图2
     图2
    图3
     图3
    图4
    图4
    本文主要对文中的代码做了较为详细的注释,并针对代码,先对公式进行了推导,其中公式中的矩阵和代码中是一一对应的,在学习代码时,需要仔细对照,希望对各位学习MPC仿真的同学有用,谢谢~
    参考1:龚建伟《无人驾驶车辆模型预测控制》

    展开全文
  • 代码中基本上每一句都有注释,由于第五章所采用的输出方程较第四章有所不同,所以本代码中参考相关论文进行详解
  • %预测时域 Nc=30;%控制时域 Row=10;%松弛因子 fprintf('Update start,t=%6.3f\n',t) t_d =u(3)*3.1415926/180; % t_d为横摆角,CarSim输出的为角度,角度转换为弧度 %直线路径 % r(1)=5*t; % r(2)=5; % r(3)=0; % ...

    可以先进行对MPC的推导,推导过程中注意矩阵的维度以及矩阵的构造,最后将我们的目标函数转化为二次规划问题,构造二次规划中所需要的矩阵,用matlab中二次规划求解函数quadprog进行求解。

    二次规划问题:https://blog.csdn.net/m0_50888394/article/details/115257655

    MPC推导过程:https://blog.csdn.net/m0_50888394/article/details/115556185

    % carsim中的输出端为X坐标(m)、Y坐标(m)、横摆角(°)、质心处的纵向速度(km/h)、方向盘转角(°)
    % Jeossiery 2021.4.13
    function [sys,x0,str,ts] = MY_MPCController3(t,x,u,flag)  %构建S-function函数,t为采样时间,x为状态变量,u是输入(simulink模块的输入),flag为仿真过程中的标志位
    switch flag  %sys输出根据flag的不同而不同
     case 0
      [sys,x0,str,ts] = mdlInitializeSizes; % Initialization初始化
     case 2
      sys = mdlUpdates(t,x,u); % Update discrete states 更新离散状态
     case 3
      sys = mdlOutputs(t,x,u); % Calculate outputs 计算输出
     case {1,4,9} % Unused flags 滞空
      sys = [];
     otherwise
      error(['unhandled flag = ',num2str(flag)]); % Error handling
    end
    % End of dsfunc.
    %==============================================================
    % Initialization
    %==============================================================
    
    function [sys,x0,str,ts] = mdlInitializeSizes
    sizes = simsizes;
    sizes.NumContStates  = 0; %连续状态量
    sizes.NumDiscStates  = 3; %离散状态量
    sizes.NumOutputs     = 2; %输出的个数
    sizes.NumInputs      = 3; %输入量的个数
    sizes.DirFeedthrough = 1; % Matrix D is non-empty.
    sizes.NumSampleTimes = 1; %采样时间的个数
    sys = simsizes(sizes); %设置完后赋值给sys输出
    x0 =[0;0;0]; %状态量初始化,横纵坐标以及横摆角都为0  
    global U; %设定U为全局变量
    U=[0;0]; %初始的控制量为0,0
    % Initialize the discrete states.
    str = [];             % Set str to an empty matrix.str为保留参数,mathworks公司还没想好怎么用它,一般在初始化中将其滞空即可
    ts  = [0.1 0];       % ts为1*2维度的向量,采样周期+偏移量sample time: [period, offset]
    %End of mdlInitializeSizes
    		      
    %==============================================================
    % Update the discrete states
    %==============================================================
    function sys = mdlUpdates(t,x,u) %flag=2表示此时要计算下一个离散状态,即x(k+1)
      
    sys = x;
    %End of mdlUpdate.
    
    %==============================================================
    % Calculate outputs
    %==============================================================
    function sys = mdlOutputs(t,x,u) %flag=3表示此时要计算输出,如果sys=[],则表示没有输出
        global a b u_piao;  %a,b和u_piao矩阵为全局矩阵
        global U;  %U全局控制量,2*1 维度矩阵
        global kesi; %kesi为全局,新的状态向量,[k时刻的状态量误差;k-1时刻的控制量误差],5*1矩阵
        tic
        Nx=3;%状态量的个数
        Nu =2;%控制量的个数
        Np =60;%预测时域
        Nc=30;%控制时域
        Row=10;%松弛因子
        fprintf('Update start,t=%6.3f\n',t)
        t_d =u(3)*3.1415926/180; % t_d为横摆角,CarSim输出的为角度,角度转换为弧度
       %直线路径
    %     r(1)=5*t;
    %     r(2)=5;
    %     r(3)=0;
    %      vd1=5;
    %      vd2=0;
    
        %半径为25m的圆形轨迹,速度为5m/s
        r(1)=25*sin(0.2*t);   %参考的X
        r(2)=25+10-25*cos(0.2*t); %参考的Y
        r(3)=0.2*t; %  参考的横摆角,0.2怎么算出来的? 答案可以根据书中p106 公式推出,根据车辆运动学的公式带入速度,前轮转角和轴距即可算出
        vd1=5; %速度5m/s
        vd2=0.104;  %前轮偏角=tan(轴距/半径)=tan(2.6/25)=0.104
        
    %     %半径为25m的圆形轨迹,速度为3m/s
    %     r(1)=25*sin(0.12*t);
    %     r(2)=25+10-25*cos(0.12*t);
    %     r(3)=0.12*t;
    %     vd1=3;
    %     vd2=0.104;
    
    	%半径为25m的圆形轨迹,速度为10m/s
    %      r(1)=25*sin(0.4*t);
    %      r(2)=25+10-25*cos(0.4*t);
    %      r(3)=0.4*t;
    %      vd1=10;
    %      vd2=0.104;
    
    %     %半径为25m的圆形轨迹,速度为4m/s
    %      r(1)=25*sin(0.16*t);
    %      r(2)=25+10-25*cos(0.16*t);
    %      r(3)=0.16*t;
    %      vd1=4;
    %      vd2=0.104;
        kesi=zeros(Nx+Nu,1); %构造新的状态量矩阵 ,5*1维矩阵
        kesi(1)=u(1)-r(1);%u(1)==X(1)  将状态量误差放到相应的位置上,第一行X的误差
        kesi(2)=u(2)-r(2);%u(2)==X(2)  第二行Y的误差
        kesi(3)=t_d-r(3); %u(3)==X(3)  第三行是横摆角的误差
        kesi(4)=U(1); % 速度误差放到第五行,这里其实要减去参考的控制量,由于我们设置的参考控制量为0所以等价
        kesi(5)=U(2); %前轮偏角误差放到第五行
        fprintf('Update start, u(1)=%4.2f\n',U(1))
        fprintf('Update start, u(2)=%4.2f\n',U(2))
    
        T=0.1; %采样时间为0.1s,即100毫秒
        T_all=40;%临时设定,总的仿真时间
        L = 2.6; %轴距为2.6m
           
    %矩阵初始化   
        u_piao=zeros(Nx,Nu); %构造3*2维度的矩阵来存放控制量误差
        Q=100*eye(Nx*Np,Nx*Np);  %权重矩阵Q为180*180的单位矩阵
        R=5*eye(Nu*Nc); %权重矩阵R为60*60的单位矩阵
        a=[1    0   -vd1*sin(t_d)*T;
           0    1   vd1*cos(t_d)*T;
           0    0   1;]; %a为我们线性离散化后的第一个系数矩阵
        b=[cos(t_d)*T   0;
           sin(t_d)*T   0;
           tan(vd2)*T/L      vd1*T/((cos(vd2)^2)*L);]; %b为我们线性离散化后的第二个系数矩阵
        A_cell=cell(2,2); % 构建2*2的元胞数组
        B_cell=cell(2,1); %构建2*1的元胞数组
        A_cell{1,1}=a;  %将a矩阵放到A_cell的第一行第一个位置
        A_cell{1,2}=b;  %将b矩阵放到A_cell的第一行第二个位置
        A_cell{2,1}=zeros(Nu,Nx); %将2*3的零矩阵放到A_cell第二行的第一个位置
        A_cell{2,2}=eye(Nu); %将2*2的单位阵放到A_cell第二行的第二个位置
        B_cell{1,1}=b;%将b矩阵放到B_cell的第一行
        B_cell{2,1}=eye(Nu);%将2*2的单位阵放到B_cell第二行
        A=cell2mat(A_cell); %这里的A就是我们在推导下一时刻的状态空间时候的A,详见CSDN推导过程
        B=cell2mat(B_cell); %这里的B就是我们在推导下一时刻的状态空间时候的B
        C=[1 0 0 0 0;0 1 0 0 0;0 0 1 0 0;];%这个C矩阵是我们输出方程yita的系数矩阵,因为我们可能不需要把每个状态量都输出所以就可以通过设置C矩阵来输出我们想要的状态量,在这里我们输出的是X的误差、Y的误差以及横摆角的误差
        PHI_cell=cell(Np,1);%这个PHI是我们通过总结规律得到的等式右边的第一个系数矩阵,60*1维度
        THETA_cell=cell(Np,Nc);%这里的THETA为我们通过总结规律的到的等式右边的第二个系数矩阵,60*30维度,具体请详见我们CSDN的推导过程
        for j=1:1:Np
            PHI_cell{j,1}=C*A^j; %通过循环来给第一个系数矩阵赋值
            for k=1:1:Nc
                if k<=j
                    THETA_cell{j,k}=C*A^(j-k)*B; %C为3*5矩阵;A为5*5;B为5*2,所以C*A*B为3*2矩阵
                else 
                    THETA_cell{j,k}=zeros(Nx,Nu); %详见CSDN推导过程
                end
            end
        end
        PHI=cell2mat(PHI_cell);%size(PHI)=[Nx*Np Nx+Nu],180*5维度
        THETA=cell2mat(THETA_cell);%size(THETA)=[Nx*Np Nu*Nc+1] 180*
        H_cell=cell(2,2); %这里的H为我们二次规划中的H矩阵,以下来构造二次规划中的H矩阵
        H_cell{1,1}=THETA'*Q*THETA+R;
        H_cell{1,2}=zeros(Nu*Nc,1); %60*1维度
        H_cell{2,1}=zeros(1,Nu*Nc); %1*60维度
        H_cell{2,2}=Row;  %H矩阵的右下角的元素就只有一个就是我们的松弛因子
        H=cell2mat(H_cell); %由于松弛因子的影响,最终的H矩阵为61*61
    
         error=PHI*kesi;%这里的error就是我们所设的E矩阵
        f_cell=cell(1,2); %f为二次规划的第二个向量,下面我们来构造它
        f_cell{1,1}=2*error'*Q*THETA; 
        f_cell{1,2}=0; %详见CSDN推导过程
    %     f=(cell2mat(f_cell))';
        f=cell2mat(f_cell);%将元胞数组转化为矩阵
        
     %% 以下为约束生成区域
     %不等式约束
        A_t=zeros(Nc,Nc);%见falcone论文 P181
        for p=1:1:Nc
            for q=1:1:Nc
                if q<=p 
                    A_t(p,q)=1;
                else 
                    A_t(p,q)=0;
                end
            end 
        end 
        A_I=kron(A_t,eye(Nu));%对应于falcone论文约束处理的矩阵A,求克罗内克积
        Ut=kron(ones(Nc,1),U);%此处感觉论文里的克罗内科积有问题,暂时交换顺序
        umin=[-0.2;-0.54;];%维数与控制变量的个数相同
        umax=[0.2;0.332];
        delta_umin=[-0.05;-0.0082;]; %源代码有错,速度下界没加负号
        delta_umax=[0.05;0.0082];
        Umin=kron(ones(Nc,1),umin);
        Umax=kron(ones(Nc,1),umax);
        A_cons_cell={A_I zeros(Nu*Nc,1);-A_I zeros(Nu*Nc,1)};
        b_cons_cell={Umax-Ut;-Umin+Ut};
        A_cons=cell2mat(A_cons_cell);%(求解方程)状态量不等式约束增益矩阵,转换为绝对值的取值范围
        b_cons=cell2mat(b_cons_cell);%(求解方程)状态量不等式约束的取值
       % 状态量约束
        M=10; 
        delta_Umin=kron(ones(Nc,1),delta_umin);
        delta_Umax=kron(ones(Nc,1),delta_umax);
        lb=[delta_Umin;0];%(求解方程)状态量下界,包含控制时域内控制增量和松弛因子
        ub=[delta_Umax;M];%(求解方程)状态量上界,包含控制时域内控制增量和松弛因子
        
        %% 开始求解过程
        %options = optimset('Algorithm','active-set'); %新版quadprog不能用有效集法,这里选用内点法
        options = optimset('Algorithm','interior-point-convex'); 
        [X,fval,exitflag]=quadprog(H,f,A_cons,b_cons,[],[],lb,ub,[],options);
        %% 计算输出
        u_piao(1)=X(1); % X1为二次规划求出来的是控制时域内最优的速度控制增量
        u_piao(2)=X(2); % X2为二次规划求出来的最优前轮转角控制增量
        U(1)=kesi(4)+u_piao(1); %kesi4为速度误差与速度控制增量相加再与下面的参考速度相加才是真正的速度
        U(2)=kesi(5)+u_piao(2); %kesi5为前轮偏角误差,道理同上
        u_real(1)=U(1)+vd1;
        u_real(2)=U(2)+vd2;
        sys= u_real;
        toc
    % End of mdlOutputs.

     

    展开全文
  • 北京理工无人驾驶车辆模型预测控制代码全集。 Matlab m文件:跟踪直线代码详解,跟踪弯道代码 联合仿真:跟踪圆形轨迹代码详解,双移线轨迹的代码详解+参考论文,局部路径规划+轨迹跟踪代码+参考论文,考虑轮胎滑移...
  • 卡尔曼滤波原理详解及系统模型建立(simulink)

    万次阅读 多人点赞 2019-12-06 11:27:56
    卡尔曼滤波原理详解及系统模型建立(simulink) 卡尔曼滤波器 卡尔曼滤波器是在上个世界五六十年代的时候提出的,到今天已经有六十年左右的时间,但卡尔曼滤波算法不管在控制、制导、导航或者通讯方面对数据的预测...

    卡尔曼滤波原理详解及系统模型建立(simulink)

    卡尔曼滤波器

    卡尔曼滤波器是在上个世纪五六十年代的时候提出的,到今天已经有六十年左右的时间,但卡尔曼滤波算法不管在控制、制导、导航或者通讯方面对数据的预测能力依然处在一个不可撼动的位置上,可是很多人对于其算法内部的工作原理究竟是怎么运作的依然不理解,所以在工程上很多人都只是把卡尔曼滤波当成是一种“黑箱”预测算法,并不清楚内部原理。但实际上没有任何算法是“黑箱”,只是算法内部的运行规律并不直观,所以让人很难理解,现在也有很多对卡尔曼滤波的解释,但是我这篇文章里希望从原理入手,尽可能定性地对卡尔曼滤波的每一步都做出更加通俗的解释,最后对卡尔曼滤波的系统过程建立相对应的模型,对其进行各种响应的测试,这样也能够更深入地理解卡尔曼滤波。

    原理

    首先我们要明白的就是,卡尔曼滤波不是传统意义上我们理解的滤波,这是一种算法,在变化的数据中去除噪声对系统未来输出做出预测的算法,是基于概率统计原理的预测算法,是一整个系统,是系统就有输入和输出,所以第一步,我们要明确卡尔曼滤波的输入输出到底是在干啥玩意儿。

    很简单看下面这个图:
    卡尔曼滤波干的事
    KF就是卡尔曼滤波,算法的输入值是一个可测的量,这个量可以是任何量反正得你能测量出来,而且还知道这个测量值的精度大概在多少,有了这个测量值我们就根据测量值来估计这个系统的真实输出,并同时给出我新估计的这个值的精度大概在什么范围内,这就是卡尔曼滤波做的工作,但这个工作是不断进行的,对系统不断测量,然后不断估计,这样持续一段时间之后就能估计出系统一个非常准确的输出值。这里要明确的一点是,测量值可能非常不准确,估计值也非常不准确,这符合工程中的很多工作状况,但仅仅根据这两个不准确的值最后就可以估计出一个相对准确的系统输出值,这也就是卡尔曼滤波的作用。

    下面就来具体讲卡尔曼滤波是咋运行的

    网络上有很多例子,小车啊室内温度的估计还是估计小兔兔体重啥的,但这里为了形象我还是要借用一下这个小车的例子,方便我们讲原理:

    有一辆小车车,在一个水平轴上向右行驶,它的初始位置我是用脚测量的,非常不准,我也只能确定这个位置在一个高斯分布的范围内,如图:
    在这里插入图片描述
    当然了我还有眼睛,大概能看出来这车的速度,现在我就可以根据我刚刚脚测得的位置和我目测的速度,大概估计下一秒他的位置,它大概运动到了这个位置:在这里插入图片描述
    这图非常合理,为啥呢,你会看到高斯分布的方差变大了,因为这是我估计出来的位置,脚测量和眼睛测量的误差叠加在一起所以导致我估计的这个位置的精度更加不准确,所以它是一个更加“肥胖”的高斯分布,但没关系,我再用脚测量一下不就好了很简单,好了我又用脚测量了一下,如图:

    在这里插入图片描述
    这里黄色的是我估计的位置分布,蓝色的是我此时重新测量的,那车车真实的位置在哪里呢?不知道,我永远不会知道,但是我可以同时相信这两个位置分布,然后估计出一个最优的位置分布,这里蓝色的高斯分布更加瘦,所以就更加值得我信赖,所以我会信它更多一点,那我们最终确定车车的位置就在这里:
    在这里插入图片描述
    图中绿色的部分就是我现在能信任的最准确的小车的位置,它是由我估计出来和测量出来的值共同估计出的,也就是卡尔曼滤波的输出值,而下一步我会根据这个最优估计值进行估计,(也就是我前面讲的“脚测”,但是这个时候脚的估计值精度已经明显提高了)然后再用脚测量一次,然后根据两个值的“可信赖度”去决定我信谁更加多一点,然后估计出一个最优估计值,那你就会问了,我怎么知道我该信谁多少呢?具体要信多大的量呢?这也就是卡尔曼最tm牛逼的地方,他给出了卡尔曼增益,告诉了大家每次估计完成之后应该更偏向于哪个值具体多少。

    前方劝退警告

    如果只想对卡尔曼滤波有一个感性的认识可以直接跳过这一部分去看我最后的系统建模,但是我觉得如果想对卡尔曼滤波有个深入的理解就跟我一起耐心地把整个过程推导一遍(主要是因为老子写得很辛苦,最好看一下),虽然看起来非常多,但是慢慢啃下来肯定会有收获的,这里需要有一些矩阵求导和控制原理相关的基础,如果不太懂可以去查一下,我也会尽量用通俗的语言来讲每一步:

    正式开始推导

    首先要问问自己我们到底想要啥,是不是就想要每一个时刻(k时刻)的系统真实状态,这里就用X(k)来表示k时刻系统的真实状态矩阵:
    在这里插入图片描述
    这里面A是状态转移矩阵,w(k)是噪声矩阵,认为是高斯分布的,X就是我们想要的每个时刻的系统状态的真实值,这个值永远不能得到,U是系统的输入向量,这个式子不是卡尔曼滤波的内容,这其实是马尔科夫提出的理论:一个系统的某些因素在转移中,第n次结果只受到第n-1次的结果影响,这里很好理解,如果知道前一时刻的系统状态,那我们可以根据这个系统的变化方式计算出这一时刻的系统状态。这里怎么去理解状态转移矩阵A,简单来说就是这个系统的变化方式,比如如果系统里变化量我们认为是室温,那X(k)就表示k时刻的室温,而状态转移矩阵这个时候就是1,因为我们认为温度是不变化的,如果X(k)是小车的运动速度、加速度、位移,那A就是速度和加速度的变化关系,位移和速度还有加速度的变化关系,就是那些高中的运动公式。

    下面我们在k时刻对系统的真实值进行测量,也就是下面公式里的Z(k),这时候的测量值就是我们可以得到的了,H是测量矩阵,比如我们要测的是小车的位移,那给状态变量矩阵左乘一个测量矩阵就得到了我们要的位移。v(k)是测量误差,也是高斯分布。
    在这里插入图片描述
    到这里为止,上面的两个式子都不是卡尔曼滤波的内容,上面是两个一定成立的等式,下面就要开始正式卡尔曼了:

    我们先看他的第一个等式:
    在这里插入图片描述
    这里面X(k-1|k-1)表示前一时刻系统的卡尔曼滤波最优估计值,咦,你会问,这还没开始预测呢哪来的最优估计值,因为卡尔曼滤波本身就是一个循环迭代的过程,每次估计都要用到上一次的最优估计值,(第一次用的是测量值,随便估计一个就可以了),所以就先假装已经有一个了(后面再讲怎么来的),而X(k|k-1)就是我们根据上一刻的最优估计值左乘一个状态转移矩阵再加上控制量得到的,这里因为我们已经估计值了所以没有噪声项。
    卡尔曼滤波的第二个等式:
    在这里插入图片描述
    这个实际上是一个假象的式子,先宏观来看,卡尔曼这个式子实际上就是用我们根据上一时刻估计的系统状态来修正一下,最终得到此刻系统输出的最准确的估计值,怎么估计呢,就是最开始那个小车的例子了,先假设有一个卡尔曼增益K(k),这个增益会在此刻衡量人为估计值和测量值哪个更可信,然后做一个决策,得到最终最优卡尔曼估计值,在这里我们先不要纠结这个值到底该怎么取,这一步要先理解整个式子是怎么假设出来的就可以了,简单来说就是假设了一个参数去衡量我该信哪个更多得到最终的估计值。

    在引出卡尔曼滤波的后三个关键等式之前:

    先来看看卡尔曼究竟是如何取到这个卡尔曼增益,最终让估计值更加精确的,要明白我们在干啥,看回最开始的小车的例子,我们每次希望得到的是最后这张图:
    在这里插入图片描述
    也就是说我们想要的是绿色的这个分布,我们不仅要得到位置,还要求位置的精度足够高,就相当于我们希望绿色高斯分布的方差尽量小,这里记P(k|k)是k时刻我们估计出的最优值的协方差,好了,我们的目标就是:最小化P(k|k)。

    那P(k|k)哪里来?

    所谓卡尔曼最优估计协方差就是最优估计与真实值的误差,那们记e(k|k)为k时刻最优估计的误差(后面e(k|k-1)是根据前时刻的估计值误差),显然该误差等于最优估计和系统真实输出的差即:e(k|k) = X(k)-X(k|k)

    但X(k|k)可以用卡尔曼滤波的第二个等式表示,所以对e(k|k)进行化简:在这里插入图片描述

    化到这一步以后,你会发现K(k)也就是卡尔曼增益依然存在,还是没有解决问题,没有告诉大家这个值该怎么取,但是这里出现了一个有意思的东西了:(X(k)-X(k|k-1)),这是啥玩意?刚刚我们说X(k)-X(k|k)是k时刻最优估计的误差,那这里显然就是人为估计的误差,也就是我们说的修正之前的误差,也就是说,这一步可以看出来,可以用修正前的估计误差和卡尔曼增益得到最优估计的误差。

    那下面就来计算最优估计协方差:
    在这里插入图片描述
    可以看到从第四个等号开始用了P(k|k-1)来代替E[(X(k)-X(k|k-1))(X(k)-X(k|k-1))^T],直接用修正前的估计协方差来替代这一坨乱七八糟的。那现在目标就很明确了,现在就要来最小化这个最优估计协方差,问题就转化成:将卡尔曼增益视为自变量,使得这个矩阵的迹最小时的卡尔曼增益就是我们要的K(k)。(最后一个等号把白噪声的协方差用R表示了,化简过程中要注意白噪声的期望是0,所以中间两项消去了)

    那我们就来取最优估计协方差矩阵的迹:
    在这里插入图片描述
    求最小值就顺便对K(k)求个导:
    在这里插入图片描述再顺便让导数等于零:
    在这里插入图片描述
    然后顺便把K(k)拿出来:
    在这里插入图片描述
    好了,到这里终于得到了我们朝思暮想的卡尔曼增益,可以发现是由未修正估计的协方差和测量误差的协方差决定的,以上也是卡尔曼滤波后三个等式中的一个,到这里,我们只要想个办法求出未修正估计值的协方差就大功告成了。同样的我们记e(k|k-1)为未修正估计值的误差,而且e(k|k-1) = X(k)-X(k|k-1),这步很简单就不讲了。

    同样的,求该协方差值:

    在这里插入图片描述
    第二行开始把k时刻系统真实输出用k-1时刻的真实输出来表示,同样把噪声的期望认为是0,这样就化简到了倒数第二个等号,可以看到又出现了一个神奇的东西!X(k-1)-X(k-1|k-1),这不就是e(k-1|k-1)吗,这个东西我们刚刚算过是不,那就替换过来,就有了红色的等式,噪声协方差用Q表示,这时候就很神奇了,估计协方差和前时刻的最优估计协方差是一个转移矩阵平方的关系,顺便加上噪声协方差,所以到这里为止我们可以用前时刻最优估计协方差来表示现在的估计协方差,也可以用现在的估计协方差表示现在的最优估计协方差,所以很明显,协方差的估计可以用来迭代了,红色的等式也就是卡尔曼滤波的第四个等式,就快要胜利了。

    整个完整的卡尔曼滤波的迭代过程这里已经可以写出来了,但你问为啥这里只有四个式子还少一个呢?其实我们最开始已经有第五个式子了,只是还不够好看,接下来就把刚刚没做完的做完:
    在这里插入图片描述
    这里就是把括号打开了哈,仔细化一下很简单的,然后把右边的卡尔曼增益项用刚刚的式子替代:
    在这里插入图片描述
    就化简出了一个非常简单的结果了,这就是卡尔曼最优估计值的协方差,只用估计协方差和卡尔曼增益就可以表示出来了,这里R被消去不是说跟估计值的协方差噪声没有关系了,而是把R放进卡尔曼增益里表示更加简洁,这也就是卡尔曼滤波的第五个式子了!到这里完整的卡尔曼滤波已经被推导出来了,喜大普奔!!!完结撒花!!!

    总结一下
    卡尔曼滤波的五个等式:
    在这里插入图片描述
    前两步是根据上一个最优估计值得出此刻的估计值和估计值的协方差,紧接着就可以得到此刻的最优估计值和最优估计值的协方差,然后利用此刻的最优估计值和最优估计值的协方差进行下一个迭代,完美。

    恭喜跟我看到这里的童鞋们,我真是谢谢您嘞,说明我也没白忙这么半天。

    我们可以再回头看一下卡尔曼滤波的五个等式,看起来有点多,但是实际上非常紧凑,推理非常严谨,但是有一些我也没想清楚卡尔曼究竟是怎么想到要这么化简的,觉得是巧合但又不是,只能感慨数学的伟大,卡尔曼的牛逼。

    开始系统建模
    这一部分开始最好随时把刚刚推导出来的五个卡尔曼滤波的等式放在一个方便看的地方,配合我的模型来边看边理解会比较好。

    下面就根据刚刚推导出的五个卡尔曼滤波的公式来建立一个simulink模型,模型如下:
    在这里插入图片描述

    这个模型要怎么看,就分两部分看,上面半部分实际上就是卡尔曼滤波的前两个公式,也就是预测值更新的过程,下面一大坨都是卡尔曼增益更新的迭代过程,每个节点我都把具体的变量标清楚了,仔细看是可以看出卡尔曼滤波五个公式的影子的,这里我默认是不加入控制量的,也方便后面讲整个系统的响应,所以我这里的状态转移矩阵和测量矩阵都是1,也就是用卡尔曼滤波来预测像温度变化这种简单的物理量变化,温度变化一般就认为是不变化的,所以状态转移矩阵也就是1,而温度也是直接可测的,所以测量矩阵也就是1。

    所以有了模型以后就很清楚了,卡尔曼滤波的五个等式看起来参数非常多,但大多数参数都在不断的循环更新中,所以整个系统真正能调整的参数只有三个,就是系统开始循环的初始值,这个要人为定一个。还有就是最最重要的两个参数,过程误差和测量误差的协方差,这就比较奇怪了,能调整的参数居然是两个噪声的协方差?!但是要知道过程误差和测量误差很多时候是不可得的,是系统本身就有的,所以我感觉这就是卡尔曼滤波很神奇的地方,即便是对误差的协方差有一个人为的不准确的判断,也能对系统的真实状态有个比较好的预测。

    下面我们就来对卡尔曼滤波的完整系统进行一个深入的研究,首先还是要把系统建立得更完整一点,我们给定一个真实系统,加上一个噪声,再对真实系统的输出进行测量,再加上测量噪声,将总的测量信号接入卡尔曼滤波,然后同时检测系统真实输出和测量值还有卡尔曼最优估计值,完整系统如下:
    在这里插入图片描述

    然后加入一个比较经典的阶跃输入信号,先去除所有噪声,系统过程噪声和测量噪声,包括卡尔曼滤波内部的所有协方差值都设置为0,看看卡尔曼滤波的输出:
    在这里插入图片描述
    这时候只有测量值可以检测到,没有卡尔曼输出,原因灰常简单,当R等于0的时候,卡尔曼增益等于测量矩阵的逆,这样P(k|k)=[I-K(k|k)H]P(k|k-1)就永远等于0,卡尔曼输出没有任何意义,所以要给R赋一个很小的初值才能认为传感器的精度非常高。同样的,先验误差如果等于0,也就是Q等于0,此时卡尔曼增益等于0,也就是说不考虑测量值,卡尔曼滤波只相信自己的估计值,但是你想想哈,卡尔曼滤波的输入就是测量值,你又不信它,太过分了吧,您凭空想象吗,所以这个时候的所有估计值都是没有任何意义的,即便我们希望某个误差很小,也必须是一个不等于0的值。

    那首先我们将先验误差和后验误差的协方差设同时设置为0.001,同样输入一个阶跃信号,看卡尔曼滤波的输出:
    在这里插入图片描述
    你问这跟刚那个图不是一毛一样么,不是的哈不信你翻回去瞅瞅,刚那个是蓝色的,这时候卡尔曼滤波已经有输出了,只不过卡尔曼滤波的输出和测量值重合了,为啥呢,也好理解,看卡尔曼增益:
    在这里插入图片描述
    其实这个时候先验误差影响是不大的,在传感器误差协方差接近于0的时候,就是我们认为测量误差非常非常小的时候,也就是认为传感器是一个神仙传感器,它测量出来的值是绝对真实可信的,这是一个最准的传感器,于是卡尔曼增益这里分母的R是一个高阶小量,分子约掉,最后卡尔曼增益等于测量矩阵的逆,最终看回卡尔曼滤波的第二个式子,会发现卡尔曼滤波的最优估计值永远等于测量值,这很好理解,算法都认为测量值准确无误了,那必然每次估计出的值只要相信传感器就好了呀,所以在这种情况下,它会复制传感器的值。

    下面我们固定先验误差的值不变,保持0.001,然后改变测量误差协方差的值,改变幅度很小的时候其实变化也并不大,我这里直接展示改变到10和1000的变化较大的时候。

    当R为10时:
    在这里插入图片描述
    会发现最优估计值和测量值开始有一点偏移,这里实际上也是比较好理解的,增大了R值,系统在预测过程中对测量值的信赖度降低,所以接近测量值的速度会降低,同样的把R值增大结果类似,下面R变为1000:
    在这里插入图片描述
    可不就是么,收敛速度会降低,如果固定R不动,增大先验误差的协方差道理也是一样的,会使收敛速度加快,因为增大Q的同时也就是减小了R,卡尔曼增益就更大了,卡尔曼增益变大,测量余量影响的比例就越大,所以测量值起到的作用也就越大,到这里Q和R的互相关系也就差不多描述清楚了。总结一下就是Q和R都不能为0,卡尔曼滤波的输出和先验误差、测量误差的协方差的具体值是没有关系的,和它们的比值高度相关,当然这都是从原理出发,实际在使用的时候还是要对传感器的漂移有个实际的估计。

    但是到这里为止,系统和传感器都是不真实带有误差的,我只不过在卡尔曼滤波器的内部主观加上了误差的协方差,下面就来测试一个真实带有噪声的系统。现在给系统过程误差加上0.01的方差,传感器本身也带有0.01的方差。系统输出和测量值如图所示:
    在这里插入图片描述红色线是卡尔曼滤波的值,还没加进去,蓝色是测量值,橙色是系统的真实输出,可以看到测量值围绕真实值上下飘动,下面开始用卡尔曼滤波预测,首先假设卡尔曼滤波内部对先验误差和测量误差的估计都是准确的,这里设置Q和R都是0.01,卡尔曼估计结果如下:

    在这里插入图片描述
    红色线围绕真实输出上下飘动,放大点可以看出来卡尔曼滤波可以得到一个比测量值更加接近真实值的输出:
    在这里插入图片描述

    这里我直接讲一个我的实验结论,就是刚刚我提到卡尔曼增益和单独的Q和R取值关系不大,只和它们的比值有关,所以这里如果同时给R和Q增大相同的倍数,卡尔曼滤波的估计结果都差不多,和真实输出的残差范围也基本没有变化,这个结果就不在这里展示了,直接看上面的图就可以,结果都类似于上面这个实验,这里我就展示一下输出结果有明显不同的实验。

    (下面的仿真对系统真实的过程误差和传感器的真实测量误差都没有改变,只改变卡尔曼滤波里的方差)

    下面固定Q不变,增大或者减小R值观察卡尔曼滤波的结果。

    首先固定Q为0.01不变,将R增大到1:
    在这里插入图片描述
    放大来看:
    在这里插入图片描述

    卡尔曼滤波的估计值收敛速度明显变慢了,但收敛后的值也比较平缓,原因是R增大,卡尔曼增益减小,对传感器的输入值不能足够相信,但更信任自己的估计值,所以每次卡尔曼估计出的结果都不会有太大的振动,因为每次都是根据上次自己的估计值来估计的,而且自己很信任自己的估计,所以很平稳,只是收敛速度慢。

    下面固定Q为0.01不变,将R减小为0.0001
    在这里插入图片描述
    放大看:
    在这里插入图片描述
    会发现看不到测量值了,其实这是收敛速度过快的表现,此时卡尔曼增益非常大,所以卡尔曼滤波系统内部对测量值的置信度是非常高的,几乎与测量值完全重合。

    做到这里就基本上可以概括卡尔曼滤波里两个重要参数对滤波器的影响了,最后来看一个正弦波输入卡尔曼滤波时R值逐渐增大的效果,不详细解释了,如果看到这里应该可以对照前面的平衡式自己看看哈,直接上图加深理解好了:R值分别为0.0001,0.01,1,10
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    最后再来总结一下卡尔曼滤波里先验误差和传感器误差的参数调节

    ①Q和R都不能取0,Q取0时意味着完全不相信测量值的输入,所以系统输出值衡为初值;R取0时意味着无条件相信测量值,所以卡尔曼滤波的输出与测量值完全重合,这两种情况下卡尔曼滤波毫无意义。

    ②讨论Q和R的单独取值同样意义不大,因为卡尔曼增益是受R和Q的比值影响的,在调节参数的时候要同时比较两个协方差的大小。

    ③卡尔曼增益的值越大,意味着越相信测量值的输出,收敛速度越快,最优估计值震荡越明显;卡尔曼增益的值越小,意味着测量值的可信度越低,越相信系统本身的估计值,最优估计值收敛速度越慢,输出越平稳。

    好了,终于完结,再撒花。

    有啥不太对的地方欢迎讨论哈,谢谢大噶!

    只收藏不点赞,生娃没有屁股蛋!

    展开全文
  • 无人驾驶车辆模型预测控制跟踪圆形轨迹代码详解!每一行都有注释!与我们的推导过程完全一致!可以让你对MPC的理解有极大的帮助!
  • 时间序列分析是一种动态数据处理的统计方法。该方法基于随机过程理论和数理统计学方法,研究...它包括一般统计分析(如自相关分析,谱分析等),统计模型的建立与推断,以及关于时间序列的最优预测控制与滤波等内容

    在时间序列分析中最常使用的一系列模型:AR、MR、ARMA,一直到ARIMA,都源于乔治·博克斯和格威利姆·詹金斯等人的一系列工作(他们的有关成果后汇集成该领域的权威经典著作【1】)。乔治·博克斯被认为是二十世纪的一代统计学大师,他有一句广为人知的名言:所有的模型都是错误的,但有一些是有用的(“All models are wrong, but some are useful”)。为了让统计模型发挥作用,深入理解并正确使用它们尤为重要。很多网上的民科或者错漏百出的教程往往把人引入歧途。时间序列分析在计量经济学研究中已经建立起了一套比较标准的流程。只有恰当地运用各种检验方法才能保证模型的使用是正确的。这其中就不得不提很多以克莱夫·格兰杰为代表的经济学家所做的贡献(经济时间序列分析大师,被认为是世界上最伟大的计量经济学家之一,2003年诺贝尔经济学奖获奖者)——例如本文将会讨论到的协整(cointegration)。

     

    一、ARIMA模型

    系列文章的前面几篇已经讨论了ARMA模型。如果要对时间序列数据应用ARMA模型,需要先保证数据是平稳的。如果数据非平稳,一般对其进行差分处理,可以得到平稳的时间序列。这也就引出了ARIMA(Auto-Regressive Integrated Moving Average)模型,它是将非平稳数据通过差分运算变成平稳后建立的模型:

    \begin{align*} ARIMA(p,d,q):\; \; Y_t^*=\beta_0&+\beta_1Y^*_{t-1}+\beta_2Y^*_{t-2}+\cdots+\beta_pY^*_{t-p} \\ &+\epsilon_t+\phi_1\epsilon_{t-1}+\phi_2\epsilon_{t-2}+\cdots+\phi_q\epsilon_{t-q} \end{align*}

    其中,Y^*_tYd阶差分。

     

    数据非平稳,可以通过差分的方法来避免“伪回归”。但差分回归会改变经济变量本身的意义。计量经济模型的建立必须依赖已有的经济学理论,才能发挥模型的解释力。例如,我们想建立,金融发展水平对经济增长(GDP)的影响模型,那么如何描述金融发展水平这个变量呢?戈德史密斯(1994)出版了《金融结构与金融发展》一书,他在书中指出“金融结构即各种金融工具和金融机构的相对规模”。他同时认为,金融结构的存在形式是多种多样的,而这些不同的形式就体现了金融发展的不同水平。注意到,根据戈德史密斯的理论,各种金融结构的规模表征了金融发展的水平。但假如我们发现金融的规模时间序列非平稳,所以对其做了差分,那么得到的就是金融增长率。注意,金融增长率和金融的规模是两个不同的意义。如果你要用增长率来建模,就首先要论证为什么可以用金融增长率来表示金融的发展水平。

     

    如何在避免“伪回归”的同时,还能分析变量本身的经济含义,有一种计量经济分析方法——“协整分析法”。如果两组序列是非平稳的,但它们的线性组合可以得到一个平稳序列,那么就说这两组时间序列数据具有协整的性质,我们同样可以把统计性质用到这个组合的序列上来。

     

    二、协整的概念与协整校验

    如果Y_tX_t是两个I(1)的过程,其中t=0,1,2,\cdots;通常来说,对于任何的\betaY_t-\beta X_t都是一个I(1)的过程。但是,对某些\beta\neq 0来说,Y_t-\beta X_t有可能是一个I(0)过程。也就是说,它有常数均值和固定方差,自相关仅取决于序列中任意两个变量之间的时间间隔。如果存在这样的\beta,就说YX是协整的,并称参数\beta是协整参数。

     

    Granger协整的定义:如果n维时间序列Y_{1t},Y_{2t}, \cdots, Y_{nt}都是d阶单整,即I(d)。存在一个向量\beta =(\beta_1, \beta_2\cdots, \beta_n),使得\beta Y'_t\sim I(d-b),其中d-b<d, 则称序列是(d-b)阶协整,记为CI(d, b)协整,\beta是协整向量。例如,d=1、b=1,d-b=0,即n维的I(1)序列,经过协整之后,\beta Y'_t\sim I(0)

     

    协整的存在(或者说非平稳时间序列经过线性组合后能够变成平稳时间序列)背后是尤其经济学解释的。例如,三月期国债利率和六月期国债利率往往都具有相同的趋势,二者之间存在协整关系;如果不这样的话,就可能引起套利行为,而套利行为又会将分道扬镳的二者拉回到一起,最后二者的利差总是维持在一个很小的正均值(长期投资者因为损失流动性,所以相对于短期投资者会获得更高一些的回报)。可见,这种协整关系的存在是有经济学原理作为基础的。但很多时候,变量之间的关系可能没有明确的经济理论支撑;此时,我们就需要按照一定的步骤去检验变量之间是否存在协整关系。

     

    Engle-Granger协整检验——是由Robert F. Engle和Clive Granger提出的协整检验方法。首先,建立1阶单整序列Y_tX_t的回归方程:Y_t=\alpha+\beta X_t+\mu_t,并利用最小二次估计得到参数\hat{\alpha}\hat{\beta},进而得到残差序列 \hat{\mu}_t=Y_t-\hat{\alpha}-\hat{\beta}X_t;第二步,用平稳性检验方法(例如ADF)检验残差序列\hat{\mu}_t是否是平稳的。经检验,如果序列\hat{\mu}_t是平稳的,则Y_tX_t之间存在协整关系,且协整向量为(1, \hat{\beta})';否则若序列\hat{\mu}_t是非平稳的,Y_tX_t之间就不存在协整关系。

     

    三、误差修正模型

    协整使得我们对两个序列之间潜在的长期关系有所理解,利用这种协整关系,我们可以丰富处理动态模型的种类,用于经济理论的检验或预测。如果Y_tX_t是I(1)过程,但不是协整的,估计一个一阶差分形式的动态(自回归分布滞后ADL)模型:

    \Delta Y_t=\alpha_0+\alpha_1\Delta Y_{t-1}+\gamma_0 \Delta X_t +\gamma_1 \Delta X_{t-1} +\epsilon_t

    如果Y_tX_t是协整的,协整参数是\beta,便又有一个I(0)变量(组合后产生的信息)加入上面的方程中,得到误差修正模型:

    \begin{align*} \Delta Y_t &= \alpha_0+\alpha_1\Delta Y_{t-1}+\gamma_0 \Delta X_t +\gamma_1 \Delta X_{t-1} +\delta \mu_{t-1}+\epsilon_t \\ &= \alpha_0+\alpha_1\Delta Y_{t-1}+\gamma_0 \Delta X_t +\gamma_1 \Delta X_{t-1} +\delta ( Y_{t-1} -\beta X_{t-1})+\epsilon_t \end{align*}

    根据需要,误差修正模型中可以包含当前期的\Delta X_t项,也可不包含,即:

    \Delta Y_t = \alpha_0+\alpha_1\Delta Y_{t-1}+\gamma_1 \Delta X_{t-1} +\delta ( Y_{t-1} -\beta X_{t-1})+\epsilon_t

    其中,\delta ( Y_{t-1} -\beta X_{t-1})称为误差修正项。

     

    误差修正模型的优点包括:

    • 克服了伪回归的问题;
    • 保留了数据本身的经济意义,同时刻画了变量之间的长期均衡关系;
    • 模型中包含了短期动态特征。

     

     

    【本文完】

     


    参考文献与推荐阅读材料

     

    【1】Time Series Analysis: Forecasting and Control (5th Edition), George E. P. Box, Gwilym M. Jenkins, Gregory C. Reinsel, Greta M. Ljung

     

    展开全文
  • 时间序列分析之AR、MA、ARMA、ARIMA详解(4)

    万次阅读 多人点赞 2015-12-04 21:09:12
    时间序列分析是一种动态数据处理的统计方法。该方法基于随机过程理论和数理统计学方法,研究...它包括一般统计分析(如自相关分析,谱分析等),统计模型的建立与推断,以及关于时间序列的最优预测控制与滤波等内容
  • R3功能详解-物料管理

    2014-01-26 20:11:19
    1.13 R/3-数据模型 2 组织结构 2.1 集团公司 2.2 采购部门 3 基础数据 3.1 供应商 3.2 物料 3.3 采购信息记录 3.4 物料清单 3.5 分类 3.6 条件 4 物料需求计划 4.1 计划过程 4.2 批量过程 4.3 物料...
  • 1、两个最常见的衡量指标是“准确率(precision)”(你给出的结果有多少...AUC可以综合衡量一个预测模型的好坏,这一个指标综合了precision和recall两个指标。 但AUC计算很麻烦,有人用简单的F-score来代替。F-score计
  • TCPIP详解--共三卷

    2015-11-30 17:17:21
    1.8 客户-服务器模型 8 1.9 端口号 9 1.10 标准化过程 10 1.11 RFC 10 1.12 标准的简单服务 11 1.13 互联网 12 1.14 实现 12 1.15 应用编程接口 12 1.16 测试网络 13 1.17 小结 13 第2章 链路层 15 2.1 引言 15 2.2 ...
  • 第23章 Elman神经网络的数据预测----电力负荷预测模型研究 第24章 概率神经网络的分类预测--基于PNN的变压器故障诊断 第25章 基于MIV的神经网络变量筛选----基于BP神经网络的变量筛选 第26章 LVQ神经网络的分类——...
  • TCP_IP详解卷1

    热门讨论 2010-12-29 10:53:54
    《TCP/IP详解,卷1:协议》是一本完整而详细的TCP/IP协议指南。描述了属于每一层的各个协议以及它们如何在不同操作系统中运行。作者用Lawrence Berkeley实验室的tcpdump程序来捕获不同操作系统和TCP/IP实现之间传输...
  • 时间序列分析是一种动态数据处理的统计方法。该方法基于随机过程理论和数理统计学方法,研究...它包括一般统计分析(如自相关分析,谱分析等),统计模型的建立与推断,以及关于时间序列的最优预测控制与滤波等内容
  • 第一章.概述 Internet的成功 Internet体系结构被设计成支持现有网络互联,同时提供了广泛的服务与协议操作。 选用数据包的分组交换主要是因为它的鲁棒性与效率,而相对来说数据安全性...差错控制与流量控制 总体...
  • 这项新技术可以实时预测人类即将说出的内容,实时生成回应,并控制对话节奏,从而使长程语音交互成为可能。而采用该技术的智能硬件设备不需要用户在每轮交互时都说出唤醒词,仅需一次唤醒,就可以轻松实现连续对话,...
  • TCP/IP详解part_2

    2010-12-29 10:58:48
    《TCP/IP详解,卷1:协议》是一本完整而详细的TCP/IP协议指南。描述了属于每一层的各个协议以及它们如何在不同操作系统中运行。作者用Lawrence Berkeley实验室的tcpdump程序来捕获不同操作系统和TCP/IP实现之间传输...
  • TCPIP详解卷[1].part10

    2010-12-29 11:08:06
    1.8 客户-服务器模型 8 1.9 端口号 9 1.10 标准化过程 10 1.11 RFC 10 1.12 标准的简单服务 11 1.13 互联网 12 1.14 实现 12 1.15 应用编程接口 12 1.16 测试网络 13 1.17 小结 13 第2章 链路层 15 2.1 引言 15 2.2 ...
  • TCPIP详解卷[1].part07

    2010-12-29 11:04:06
    1.8 客户-服务器模型 8 1.9 端口号 9 1.10 标准化过程 10 1.11 RFC 10 1.12 标准的简单服务 11 1.13 互联网 12 1.14 实现 12 1.15 应用编程接口 12 1.16 测试网络 13 1.17 小结 13 第2章 链路层 15 2.1 引言 15 2.2 ...
  • TCPIP详解卷[1].part03

    2010-12-29 10:59:26
    1.8 客户-服务器模型 8 1.9 端口号 9 1.10 标准化过程 10 1.11 RFC 10 1.12 标准的简单服务 11 1.13 互联网 12 1.14 实现 12 1.15 应用编程接口 12 1.16 测试网络 13 1.17 小结 13 第2章 链路层 15 2.1 引言 15 2.2 ...
  • TCPIP详解卷[1].part12

    2010-12-29 11:12:15
    《TCP/IP详解,卷1:协议》是一本完整而详细的TCP/IP协议指南。描述了属于每一层的各个协议以及它们如何在不同操作系统中运行。作者用Lawrence Berkeley实验室的tcpdump程序来捕获不同操作系统和TCP/IP实现之间传输...
  • TCPIP详解卷[1].part11

    2010-12-29 11:09:56
    《TCP/IP详解,卷1:协议》是一本完整而详细的TCP/IP协议指南。描述了属于每一层的各个协议以及它们如何在不同操作系统中运行。作者用Lawrence Berkeley实验室的tcpdump程序来捕获不同操作系统和TCP/IP实现之间传输...
  • TCPIP详解卷[1].part09

    2010-12-29 11:05:16
    《TCP/IP详解,卷1:协议》是一本完整而详细的TCP/IP协议指南。描述了属于每一层的各个协议以及它们如何在不同操作系统中运行。作者用Lawrence Berkeley实验室的tcpdump程序来捕获不同操作系统和TCP/IP实现之间传输...
  • TCPIP详解卷[1].part08

    2010-12-29 11:04:39
    《TCP/IP详解,卷1:协议》是一本完整而详细的TCP/IP协议指南。描述了属于每一层的各个协议以及它们如何在不同操作系统中运行。作者用Lawrence Berkeley实验室的tcpdump程序来捕获不同操作系统和TCP/IP实现之间传输...
  • TCPIP详解卷[1].part06

    2010-12-29 11:03:10
    《TCP/IP详解,卷1:协议》是一本完整而详细的TCP/IP协议指南。描述了属于每一层的各个协议以及它们如何在不同操作系统中运行。作者用Lawrence Berkeley实验室的tcpdump程序来捕获不同操作系统和TCP/IP实现之间传输...
  • TCPIP详解卷[1].part05

    2010-12-29 11:02:35
    《TCP/IP详解,卷1:协议》是一本完整而详细的TCP/IP协议指南。描述了属于每一层的各个协议以及它们如何在不同操作系统中运行。作者用Lawrence Berkeley实验室的tcpdump程序来捕获不同操作系统和TCP/IP实现之间传输...
  • TCPIP详解卷[1].part04

    2010-12-29 11:01:49
    《TCP/IP详解,卷1:协议》是一本完整而详细的TCP/IP协议指南。描述了属于每一层的各个协议以及它们如何在不同操作系统中运行。作者用Lawrence Berkeley实验室的tcpdump程序来捕获不同操作系统和TCP/IP实现之间传输...
  • 卡尔曼滤波器是一种由卡尔曼(Kalman)提出的用于时变线性系统的递归滤波器。这个系统可用包含正交状态变量的微分方程模型来描述,这种滤波器是将过去的测量估计误差...运动预测方面,起因是博主在通过视觉算法控制PI
  • 这种描述了如何表征和控制模型预测的不确定性的概率框架,在科学数据分析、机器学习、机器人技术、认知科学以及人工智能领域中扮演着中心角色。这篇评论介绍了这种框架,并讨论了该领域的最新进展——即概率编程、...
  • 第23章 Elman神经网络的数据预测----电力负荷预测模型研究 第24章 概率神经网络的分类预测--基于PNN的变压器故障诊断 第25章 基于MIV的神经网络变量筛选----基于BP神经网络的变量筛选 第26章 LVQ神经网络的分类——...

空空如也

空空如也

1 2 3
收藏数 46
精华内容 18
关键字:

模型预测控制详解