【自适应控制】最小二乘法离线辨识

前言

想要学习自适应控制的话，需要先学习系统辨识，系统辨识中比较常见的是最小二乘法，因此后续会先介绍最小二乘法，再去完成自适应控制

理论分析：

假设一个单输入单输出的系统如下：

其离散传递函数为：

将传递函数写成差分方程：

则：

上式是一组输入与输出的关系，扩展到N组输入输出，则矩阵变为：

小知识：矩阵的转置乘以矩阵本身即为最小二乘法

最小二乘法的缺陷：

最小二乘法归根结底是时域的辨识方法，当输入信号频率增高时，在高频段，最小二乘法的辨识结果往往是错误的，但在低频段，辨识结果还是可靠的。

根据理论实际进行Matlab编程：

假定我们的传递函数Z变换之后为：

 %   Y     b1 * z^-1 + b2 * z^-2 + b3 * z^-3
 %   -  =  --------------------------------------
 %   U     1 + a1 * z^-1 + a2 * z^-2 + a3 * z^-3

对于离散传递函数，B0一般为1，转换成差分方程为：

 %   Y = - a1 * Y(k-1) - a2 * Y(k-2) - a3 * Y(k-3) +
 %         b1 * U(k-1) + b2 * U(k-2) + b3 * U(k-3) 

后续的思想就是使用原始传递函数进行扫频，之后采集数据，之后假定系统为不同的阶数，然后进行辨识，之后就可以对比不同的阶数，最终选择一个好一点的阶数：

%**********************************************假定系统为2阶**************************************************%
H2 = zeros(Sum - 2,4);   %定义一个0矩阵，因为待求参数是6个，因此为6列 (B0 = 1)
Y2 = zeros(Sum - 2,1);   %定义一个0矩阵，保存输出值  用于后续矩阵计算
%Theta2 = zeros(4,1);   %定义一个0矩阵，保存待求参数
for i = 3:1:Sum
    H2(i - 2,:) = [-TFunOut(i-1),-TFunOut(i-2),TFunIn(i-1),TFunIn(i-2)];
    Y2(i - 2) = TFunOut(i);
end
Theta2 = lsqminnorm(H2' * H2 , H2') * Y2;
IdfySys2 = tf([0,Theta2(3),Theta2(4)],[1,Theta2(1),Theta2(2)],0.001);
Out2 = zeros(1,Sum);  %计算拟合传递函数的输出曲线
In2 = zeros(1,Sum);  %计算正弦输入信号
u2_1 = 0;    %上时刻控制量
u2_2 = 0;    %上上时刻控制量
y2_1 = 0;    %上时刻输出量
y2_2 = 0;    %上上时刻输出量
for i = 1:1:Sum
    In2(i) = SinSignal(i);      %正弦输入信号
    Out2(i) = - Theta2(1)*y2_1 - Theta2(2)*y2_2 + + Theta2(3)*u2_1 + Theta2(4)*u2_2;
    u2_2 = u2_1;
    u2_1 = In2(i);    	        %前一个的控制器输出值
    y2_2 = y2_1;
    y2_1 = Out2(i);    	        %前一个的系统响应输出值    
end

%**********************************************假定系统为3阶**************************************************%
H3 = zeros(Sum - 3,6);   %定义一个0矩阵，因为待求参数是6个，因此为6列 (B0 = 1)
Y3 = zeros(Sum - 3,1);   %定义一个0矩阵，保存输出值  用于后续矩阵计算
%Theta3 = zeros(6,1);   %定义一个0矩阵，保存待求参数
for i = 4:1:Sum
    H3(i - 3,:) = [-TFunOut(i-1),-TFunOut(i-2),-TFunOut(i-3),TFunIn(i-1),TFunIn(i-2),TFunIn(i-3)];
    Y3(i - 3) = TFunOut(i);
end
Theta3 = lsqminnorm(H3' * H3,H3') * Y3;
IdfySys3 = tf([0,Theta3(4),Theta3(5),Theta3(6)],[1,Theta3(1),Theta3(2),Theta3(3)],0.001);
Out3 = zeros(1,Sum);  %计算拟合传递函数的输出曲线
In3 = zeros(1,Sum);  %计算正弦输入信号
u3_1 = 0;    %上时刻控制量
u3_2 = 0;    %上上时刻控制量
u3_3 = 0;    %上上上时刻控制量
y3_1 = 0;    %上时刻输出量
y3_2 = 0;    %上上时刻输出量
y3_3 = 0;    %上上上时刻输出量
for i = 1:1:Sum
    In3(i) = SinSignal(i); %正弦输入信号     
    Out3(i) = - Theta3(1)*y3_1 - Theta3(2)*y3_2 - Theta3(3)*y3_3 + Theta3(4)*u3_1 + Theta3(5)*u3_2 + Theta3(6)*u3_3;%系统响应输出序列           
    u3_3 = u3_2;
    u3_2 = u3_1;
    u3_1 = In3(i);    	        %前一个的控制器输出值
    y3_3 = y3_2;
    y3_2 = y3_1;
    y3_1 = Out3(i);    	        %前一个的系统响应输出值
end

%**********************************************假定系统为4阶**************************************************%
H4 = zeros(Sum - 4,8);   %定义一个0矩阵，因为待求参数是6个，因此为6列 (B0 = 1)
Y4 = zeros(Sum - 4,1);   %定义一个0矩阵，保存输出值  用于后续矩阵计算
%Theta4 = zeros(8,1);   %定义一个0矩阵，保存待求参数
for i = 5:1:Sum
    H4(i - 4,:) = [-TFunOut(i-1) -TFunOut(i-2) -TFunOut(i-3) -TFunOut(i-4) TFunIn(i-1) TFunIn(i-2) TFunIn(i-3) TFunIn(i-4)];
    Y4(i - 4) = TFunOut(i);
end
Theta4 = lsqminnorm(H4' * H4,H4') * Y4;
IdfySys4 = tf([0,Theta4(5),Theta4(6),Theta4(7),Theta4(8)],[1,Theta4(1),Theta4(2),Theta4(3),Theta4(4)],0.001);
Out4 = zeros(1,Sum);  %计算拟合传递函数的输出曲线
In4 = zeros(1,Sum);  %计算正弦输入信号
u4_1 = 0;    %上时刻控制量
u4_2 = 0;    %上上时刻控制量
u4_3 = 0;    %上上上时刻控制量
u4_4 = 0;    %上上上上时刻控制量
y4_1 = 0;    %上时刻输出量
y4_2 = 0;    %上上时刻输出量
y4_3 = 0;    %上上上时刻输出量
y4_4 = 0;    %上上上上时刻输出量
for i = 1:1:Sum
    In4(i) = SinSignal(i); %正弦输入信号    
    Out4(i) = - Theta4(1)*y4_1 - Theta4(2)*y4_2 - Theta4(3)*y4_3 - Theta4(4)*y4_4  + Theta4(5)*u4_1  + Theta4(6)*u4_2  + Theta4(7)*u4_3 + Theta4(8)*u4_4;
    u4_4 = u4_3;
    u4_3 = u4_2;
    u4_2 = u4_1;
    u4_1 = In4(i);    	        %前一个的控制器输出值
    y4_4 = y4_3;
    y4_3 = y4_2;
    y4_2 = y4_1;
    y4_1 = Out4(i);    	        %前一个的系统响应输出值    
end

%*************************画出输入正弦信号与传递函数输出信号*********************%
% subplot(2,1,1),plot(time,TFunIn);   %画数据点不连线
% grid on;

subplot(4,1,1),plot(time,TFunOut);   %原始传递函数图像
grid on;

subplot(4,1,2),plot(time,Out2);   %2阶拟合传递函数图像
grid on;

subplot(4,1,3),plot(time,Out3);   %3阶拟合传递函数图像
grid on;

subplot(4,1,4),plot(time,Out4);   %4阶拟合传递函数图像
grid on;

效果如图：