整体最小二乘参数辨识（MATLAB代码）

"整体最小二乘"通常指的是在统计学和数据分析中使用的一种方法，用于估计线性回归模型中的参数。在简单线性回归中，整体最小二乘方法旨在找到一条直线，使得观察到的数据点到这条直线的垂直距离之和最小。在多元线性回归中，该方法扩展到了多个自变量的情况。

利用整体最小二乘算法辨识模型，改变噪声大小的情况下观察辨识结果，整体最小二乘算法的辨识效果好于最小二乘算法。

。

close all
clear all
clc

J=5;
D=100;
w=10;
k=50;
L=1;
num = [0 0 k/J/w];
den = [1 D/J k/J/w];
sys  = tf(num,den)% 真实的传递函数
sysdiscrete = c2d(sys,0.02)

% 生成输入信号和输出信号的样本数据
N = 1000; % 数据点数量
Ts = 0.02; % 采样时间间隔
t = (0:N-1)*Ts; % 时间向量
u = idinput(N);  % 输入信号

figure(1)
plot(u)
title('输入随机序列')

y = lsim(sys, u, t); % 通过真实的传递函数生成输出信号
y = y + 0.00001*randn(N,1);  %假设信号采集有噪声，0.0001为调节噪声大小的参数，噪声太大则采集的数据不可靠，辨识的结果会很差

%根据差分后的方程编写最小二乘 p(k) = h1*u(k) + h2*u(k-1) + h3*u(k-2) - f1*p(k-1) - f2*p(k-2);
uk = u(3:N);
uk_1 = u(2:N-1);
uk_2 = u(1:N-2);
pk = y(3:N);
pk_1 = y(2:N-1);
pk_2 = y(1:N-2);

% 普通最小二乘

AA = [uk uk_1 uk_2 pk_1 pk_2];
BB = pk;
the = pinv(AA' * AA )*AA'*BB;

% 输出结果
disp('普通参数估计结果:');
disp('----------------------------------');
disp(['h1 = ', num2str(the(1))]);
disp(['h2 = ', num2str(the(2))]);
disp(['h3 = ', num2str(the(3))]);
disp(['f1 = ', num2str(the(4))]);
disp(['f2 = ', num2str(the(5))]);

%% 整体最小二乘

A_augmented = [AA pk];

[U,S,V] = svd(A_augmented);

S_hat = S;

S_hat(5,5)=0;

A_augmented_hat = U*S_hat*V';

A_delta = A_augmented_hat - A_augmented;

% x_hat
xx_hat_TLS = reshape(V(:,6),6,1);
xx_hat_TLS = xx_hat_TLS./(-xx_hat_TLS(6));

% 输出结果
disp('整体最小二乘参数估计结果:');
disp('----------------------------------');
disp(['h1 = ', num2str(xx_hat_TLS(1))]);
disp(['h2 = ', num2str(xx_hat_TLS(2))]);
disp(['h3 = ', num2str(xx_hat_TLS(3))]);
disp(['f1 = ', num2str(xx_hat_TLS(4))]);
disp(['f2 = ', num2str(xx_hat_TLS(5))]);


%% 绘图 
% 定义传递函数的系数

h_ls = [the(1),the(2), the(3)];
f_ls = [1, -the(4), -the(5)];
% 创建离散时间系统对象
sys1 = tf(h_ls, f_ls, 0.02, 'Variable', 'z');

h_tls = [xx_hat_TLS(1),xx_hat_TLS(2), xx_hat_TLS(3)];
f_tls = [1, -xx_hat_TLS(4), -xx_hat_TLS(5)];
% 创建离散时间系统对象
sys2= tf(h_tls, f_tls, 0.02, 'Variable', 'z');
figure()
xx = 1:N;
plot(xx,lsim(sys,u,t),'b-',xx,lsim(sys1,u,t),'r-.',xx,lsim(sys2,u,t),'k--','linewidth',1)
legend('真实系统','普通最小二乘','整体最小二乘')
xlabel('时间')
ylabel('输出')
title('带有噪声情况下辨识结果对比')

%将离散传递函数转换为连续的计算J/W等参数
c_sys1 = d2c(sys1);
[num1, den1] = tfdata(c_sys1);
c_sys2 = d2c(sys2);
[num2, den2] = tfdata(c_sys2);
% 4个参数，2个方程，需要已知2个参数，假设已知w=10;K=50;
J1 = k / den1{1}(3)/w;
D1 = den1{1}(2)*J1;

disp('电机参数估计结果:');
disp('----------------------------------');
disp('普通最小二乘算法电机参数估计结果:');
disp(['J1=',num2str(J1),'    D1=',num2str(D1)])

J2 = k / den2{1}(3)/w;
D2 = den2{1}(2)*J2;
disp('整体最小二乘算法电机参数估计结果:');
disp(['J2=',num2str(J2),'    D2=',num2str(D2)])