本人在学习柴天佑老师的自适应控制一书时，在2.3.2参数估计小节遇到了最小二乘法的参数估计，虽然书中给出了详细的推导过程，但是随后的仿真实验中并没有给出详细的代码，亦没有给出很好的解释，只是给出了仿真结果，让人摸不着头脑。在参考了大量资料后，我算是有一点点的理解，趁热打铁，赶忙写下来，也算是自己复习了。 测试使用的时间序列函数为： y(k)+1.5y(k-1)+0.6y(k-2) = 2u(k-3)-1.4u(k-4)+V(k) 1.构造出输入输出序列y(k) 、u(k) 这种非病态算法 最后肯定会趋向于某个值 所以 k可以取大一些 书中取的是150 这里取400 输出序列u(k) 可以任意给 输出序列y(k)由输入序列决定 2.定义 P 数据向量序列seita 、未知参数序列H 根据下图的算法公式 我们需要先初始化一些向量和矩阵 未知参数H（0）=0 P（0）=10^6I 其中I是单位阵 3.根据公式开始计算 先得到 数据向量序列seita 这里包括 {-y(k-1),-y(k-2),u(k-3),u(k-4)} 因为k-i大于等于0，所以k应从5开始取值 构造出这种结果 过程基本就是这样 matlab代码我会附在结尾处
下面看结果 测试使用的时间序列函数为： y(k)+1.5y(k-1)+0.6y(k-2) = 2u(k-3)-1.4u(k-4)+V(k) 输出的估计值序列变化过程如下图：
以a0为例 图中红线 基本保持在了 1.5 调出最后一组估计值
a0 = 1.5037 a1 = 0.6052 b0 = 2.0709 b1 = -1.4479 对比实际值 a0 = 1.5 a1 = 0.6 b0 = 2.0 b1 = -1.4 看起来还不错_ 改变a0的值为1.0 看是否可以继续估计成功 红线是a0的估计值 看起来也还行 调出最后一组数据 a0 = 1.0286 a1 = 0.6305 b0 = 1.9105 b1 = -1.3208 对比实际值 a0 = 1.0 a1 = 0.6 b0 = 2.0 b1 = -1.4 也算是基本保持跟随了
好了 基本情况就是这样了 如有错误 请大家批评我！ 附代码：(我不知道如何插入matlab的代码 就用C的格式了)
```%递推最小二乘辨识算法仿真 %测试使用的时间序列函数为： y(k)+1.5y(k-1)+0.6y(k-2) = 2u(k-3)-1.4u(k-4)+V(k) %待估计量 1.5 0.6 2 -1.4 %y(k) = -1.5y(k-1)-0.6y(k-2)+2u(k-3)-1.4u(k-4)+V(k) %阶数 Na = 2 Ma = 4 %先产生出 y() u() u未知 是一个行向量 clear clc x = [0 1 0 1 1 0 1 1 1]; k = 400; num = 4; u = []; for i=1:k temp = xor(x(4),x(9)); u(i) = x(9); for j = 9:-1:2 x(j) = x(j-1); end x(1) = temp; end %此处的目的是产生一个无规则的向量 可以任意赋值 此处参考了网上的一种做法 v = randn(1,400); y = zeros(400,1); y(1) = -1; y(2) = 0; y(3) = -1; y(4) = 0; y(5) = 0; for i=5:k y(i) = -1.0*y(i-1)-0.6*y(i-2)+2*u(i-3)-1.4*u(i-4)+v(i); end %y(k) = -1.5y(k-1)-0.6y(k-2)+2u(k-3)-1.4u(k-4)+V(k) %y序列 和 u序列 完成后 开始写 %还需要3个向量 %P Seita Data P = (10^6)*eye(num); P_temp = zeros(num,k); P_temp(:,4) = 10^6;%P_temp = { P(1) P(2) P(3).....} Data = zeros(num,k);%估计值序列 = {Data(1) Data(2) Data(3) ......} seita = zeros(num,1); for i = 5:k for j = 1:2 seita(j) = -y(i-j); end for l = 3:4 seita(l) = u(i-l); end %h完毕 K = (P*seita)/(seita'*P*seita+1); Data(:,i) = Data(:,i-1) + K*[y(i)-seita'*Data(:,i-1)]; P = (eye(4)-K*seita')*P; %新的P % for kk = 1:4 % P_temp(kk,i-1) = P(kk,kk); % end end d1 = Data(1,:); d2 = Data(2,:); d3 = Data(3,:); d4 = Data(4,:); i = 1:k; plot(i,d1,'r'); title('辨识结果'); hold on; plot(i,d2,'b'); hold on; plot(i,d3,'g'); hold on; plot(i,d4,'y'); hold on; legend('a0','a1','b0','b1'); a0 = Data(1,k); a1 = Data(2,k); b0 = Data(3,k); b1 = Data(4,k); a0 a1 b0 b1