用matlab最小二乘法编写程序求方程参数的值

Answer 1

function [sysd,sys,err] = ID(Y,U,Ts)
%
%基于递推最小二乘法的参数辨识程序
%仅针对二阶系统：)
%出处：http://blog.sina.com.cn/xianfa110
%=-=-=-=-=-=-=-=-=-=-==-=-=-=-=-=
%Inputs:
%=-=-=-=-=-=-=-=-=-=-==-=-=-=-=-=
%Y = nX1 vector of your model output
%U = nX1 vector of your model input
%Ts = sample time
%
%=-=-=-=-=-=-=-=-=-=-==-=-=-=-=-=
%Outputs:
%=-=-=-=-=-=-=-=-=-=-==-=-=-=-=-=
%sysd = discrete-time transfer function identified
%sys = continuous-time transfer function identified
%err = error
%
if nargin<3 || nargin>3
error('Must be three inputs!');
end
if length(Y)~=length(U)
error('length of inputs must be equal.');
end
n=length(U);
Y=reshape(Y,n,1);
U=reshape(U,n,1);
theta=[0.1;0.1;0.1;0.1;0.1];
P=2^25*eye(5);
R0=1;
for m=3:n
X=[Y(m-1) Y(m-2) U(m) U(m-1) U(m-2)]';
alfa=1/(R0+X'*P*X);
L=alfa*P*X;
theta(:,m-1)=theta(:,m-2)+L*(Y(m)-X'*theta(:,m-2));
P=P/R0-alfa*P*X*X'*P;
err=Y(m)-X'*theta(:,m-2);
if abs(err)<=1e-10
break;
end
end
m=length(theta(1,:));
result=[-theta(1:2,m);theta(3:5,m)];
t=1:m;
figure;
plot(t,theta(1,:),t,theta(2,:),t,theta(3,:),t,theta(4,:),t,theta(5,:));
legend('th1','th2','th3','th4','th5');
num=[result(3),result(4),result(5)];
den=[1,result(1),result(2)];
sysd=tf(num,den,Ts);
[n,d]=d2cm(num,den,Ts,'tustin');
sys=tf(n,d);
%%====================================================
exaple:
对于以下模型：

运行之后，数据通过Scope传递到工作空间，方法参见Simulink利用Scope输出及绘制仿真波形技巧。
输入以下代码：
Y=data(:,3);U=data(:,2);
[sysd,sys,e]=ID(Y,U,0.001);
得到结果如下：
Transfer function:
-4.314e-009 z^2 + 8.784e-005 z + 4.362e-005
-------------------------------------------
z^2 - 1.975 z + 0.9753

Sampling time: 0.001

Transfer function:
-1.12e-005 s^2 - 0.04417 s + 133.1
----------------------------------
s^2 + 25.02 s - 0.008906

e =
-6.8120e-011