matlab 最小二乘法辨识,[转载]自编最小二乘法的Matlab参数辨识程序（含实例）

乘法辨识,[转载]自编最小二乘法的Matlab参数辨识程序（含实例）"/>

matlab 最小二乘法辨识,[转载]自编最小二乘法的Matlab参数辨识程序（含实例）

function [sysd,sys,err] = ID(Y,U,Ts)

%

%基于递推最小二乘法的参数辨识程序

%仅针对二阶系统：)

%出处：

%=-=-=-=-=-=-=-=-=-=-==-=-=-=-=-=

%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

相关文章：