增广矩阵法实现最小二乘的系统辨识

《系统辨识》上机实验报告

北京工商大学

《系统建模与辨识》课程

上机实验报告

（2014年秋季学期）

专业名称 ： 控制工程 上机题目 ： 增广矩阵法实现最小二乘的系统辨识 专业班级 ：

2015 年 1 月

一、实验目的

通过仿真实验掌握增广矩阵法进行参数估计的原理和方法。

二、实验原理

考虑CARMA模型

《系统辨识》上机实验报告

其中： 还可表示为 ：

其中：

以上两向量均由2n维扩展为（2n+r）维。

在式（2）中若但是｛是已知量，则可直接用最小二乘法估计出，｝是不可测的未知量。一个简单可行的方法是用计算的代替，｛｝由下式递推得出：

其中：

（计算残差）

开始对前100组数据进行增广矩阵最小二乘进入一个新的观测数据yi，对101组数据进行参数估计剔除最早的一个数据，对余下的数据进行参数估计参数收敛情况是否满足要求NY显示被辨识参数结束开始对前100组数据进行增广矩阵最小二乘进入一个新的观测数据yi，对101组数据进行参数估计剔除最早的一个数据，对余下的数据进行参数估计参数收敛情况是否满足要求NY显示被辨识参数结束

三、实验内容

按附表9-11、表9-12给出的二阶线性离散系统的输入输出数据，用增广矩阵法进行参数估计。对任务进行方案设计，画出实验流程图，用MATLAB编程实现。撰写实验报告。

《系统辨识》上机实验报告

四、流程图

开始对前100组数据进行增广矩阵最小二乘进入一个新的观测数据yi，对101组数据进行参数估计剔除最早的一个数据，对余下的数据进行参数估计参数收敛情况是否满足要求NY显示被辨识参数结束 五、代码实现

clc

clear all; close all;

u=[1.147,0.201 ,-0.787 ,-1.584 -1.052, 0.866 ,1.152 ,1.573 ,0.626, 0.433 ...

-0.958,0.810, -0.044, 0.947, -1.474, -0.719, -0.086, 1.099, 1.450, 1.151...

0.485,1.633, 0.043 ,1.326, 1.706, -0.340, 0.890,0.433, -1.177, -0.390...

-0.982,1.435 ,-0.119 ,-0.769, -0.899, 0.882,-1.008, -0.844, 0.628,-0.679 ...

1.541,1.375, -0.984, -0.582, 1.609,0.090, -0.813, -0.428 ,-0.848,-0.410...

0.048,-1.099 ,-1.108 ,0.259,-1.627, -0.528, 0.203, 1.204,1.691, -1.235...

-1.228,-1.267, 0.309,0.043, 0.043, 1.461, 1.585, 0.552,-0.601,

《系统辨识》上机实验报告

-0.319...

0.744 0.829,-1.626, -0.127, -1.578, -0.822, 1.469, -0.379,-0.212, 0.178...

0.493 -0.056, -0.1294, 1.228, -1.606, -0.382, -0.229, 0.313,-0.161, -0.810...

-0.277 0.983, -0.288, 0.846, 1.325, 0.723, 0.713, 0.643 0.463,0.786...

1.161,0.850, -1.349, -0.596, 1.512, 0.795, -0.713, 0.453,-1.604, 0.889...

-0.938, 0.056, 0.829,-0.981, -1.232, 1.327, -0.681,0.114,-1.135, 1.284...

-1.201 0.758, 0.590, -1.007, 0.390, 0.836,-1.52, -1.053, -0.083, 0.619...

0.840 -1.258, -0.354, 0.629, -0.242,1.680, -1.236, -0.803,0.537, -1.100...

1.417,-1.024, 0.671, 0.688,-0.123, -0.952, 0.232, -0.793,-1.138, 1.154...

0.206,1.196, 1.013,1.518, -0.553, -0.987, 0.167, -1.445,0.630, 1.255...

0.311,-1.726,0.975, 1.718, 1.360, 1.667, 1.111, 1.018, 0.078,-1.665...

-0.760, 1.184, -0.614, 0.994, -0.089, 0.947, 1.706, -0.395, 1.222, -1.351...

0.231,1.425, 0.114, -0.689, -0.704, 1.070, 0.262, 1.610 ,1.489,-1.602...

0.020, -0.601, -0.020, -0.601, -0.235, 1.245, 1.226, -0.204, 0.926, -1.297];

y=[ 0.086,2.210,0.486,-1.812,-3.705,-2.688,1.577,2.883,2.883,3.705, 1.642...

0.805,-2.088,0.946,-0.039,1.984,-2.545,-1.727,-0.231,2.440,3.583...

2.915,1.443,3.598,0.702,2.638,3.611,-0.168,1.732,0.666,2.377... -0.554,-2.088,2.698,0.189,-1.633,-2.010,1.716,-1.641,-1.885,1.061...

-0.968,2.911,3.088,-1.629,-1.533,3.030,0.614,-1.483,-1.029,-1.948...

-1.066,-0.113,-2.144,-2.626,0.134,-3.043,-1.341,0.338,2.702,3.813...

-1.924,-2.813,-1.795,3.002,1.027,1.027,2.755,3.584,1.737,-0.837...

-0.617,1.703,2.045,-2.886,-0.542,-2.991,-1.859,3.045,0.068,-0.375...

0.541,1.036,0.153,-0.474,2.512,-2.681,-0.954,-0.307,0.628,-0.270...

《系统辨识》上机实验报告

-0.277,0.983,-0.288,0.846,1.325,0.723,1.750,1.401,1.340,0.916... 1.396,2.446,2.103,2.432,-1.486,3.031,2.373,-0.763,-0.752,-3.207...

1.385,-1.642,-0.118,1.756,-1.613,-1.690,2.136,-1.136,-0.005,-2.210...

2.331,-2.204,0.983,1.347,-1.691,0.595,1.809,-2.204,-2.330,-0.454...

1.290,2.080,-1.990,-0.770,1.240,-0.252,3.137,-2.379,1.206,1.221...

-1.977,2.471,-1.680,1.148,1.816,0.055,-1.856,0.269,-1.323,-2.486...

1.958,0.823,2.481,2.209,3.167,-0.762,-2.225,-0.123,-2.786,1.026...

2.843,1.071,-3.317,1.514,3.807,3.388,3.683,-1.935,-1.935,0.309... -3.390,-2.124,2.192,-0.855,-1.656,0.016,1.804,3.774,-0.059,2.371...

-2.322,-0.032,2.632,0.565,-1.460,-1.839,1.917,0.865,3.180,3.261...

-2.755,-0.536,-1.171,-0.905,-3.303,-0.834,2.490,3.039,0.134,1.901 ];

v=randn(1,200);

P=100*eye(6);%估计方差? Pstore=zeros(6,200);

Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6)]; Theta=zeros(6,200);%参数的估计值，存放中间过程估值? Theta(:,1)=[3;3;3;3;3;3];%?K=zeros(4,400);??%增益矩阵?? K=[10;10;10;10;10;10]; for i=3:200

h=[-y(i-1);-y(i-2);u(i-1);u(i-2);v(i-1);v(i-2)]; K=P*h*inv(h'*P*h+1);

Theta(:,i-1)=Theta(:,i-2)+K*(y(i)-h'*Theta(:,i-2)); P=(eye(6)-K*h')*P;

Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6)]; end

disp('参数a1、a2、b1、b2、d1、d2估计结果：') Theta(:,199) i=1:200; %输入数据

figure(1),subplot(4,1,1) stem(u) grid on

title('图1 输入信号') %输入数据

figure(1),subplot(4,1,2)

《系统辨识》上机实验报告

stem(y) grid on

title('图2 输出信号') figure(1),subplot(4,1,3)

plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:),i,Theta(6,:))

title('待估参数过渡过程') figure(1),subplot(4,1,4)

plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Pstore(5,:),i,Pstore(6,:)) title('估计方差变化过程')

六、实验结果

利用matlab仿真得到的输入输出数据、模型参数的收敛值、方差变化图如下图所示

图1 输入信号20-250-550-5100500020406080100120140160180200

图1 实验结果

020406080100120图2 输出信号140160180200050100150待估参数过渡过程200250020406080100120估计方差变化过程140160180200模型参数的辨识结果为： a1 -0.1296

a2 0.0623 b1 0.0608 b2 1.7788 d1 0.0748 d2 0.0914 《系统辨识》上机实验报告

七、实验总结

此次实验使我对闭环系统的可辨识性的概念有了进一步的深入学习，通过使用MATLAB2008a软件进行完成了系统辨识的实验，从实验结果可以看出，增广矩阵法在数据比较多的时候，方差逐步趋近于零，也就是说，递推增广矩阵法具有很好的收敛性，因此它在许多领域都得到了广泛的应用。

本文来源：https://www.bwwdw.com/article/chxp.html