最近在做關于過熱氣溫的動態建模問題。有現場運行的歷史數據，要找出導前區和惰性區的傳遞函數。

對這類算法了解不多，程序讀起來比較吃力，所以就轉來一篇完整的辨識程序，在原有基礎上進行了簡化，并稍加注解一下，也能讓有共同需要的朋友少花點時間。

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

%基于遞推最小二乘法的參數辨識程序

%僅針對二階系統(工業對象大都是高階系統，用起來很復雜，所以最終還是要降階成為二階近似系統)

%Inputs: Y是模型輸出值，U是模型輸入。二者都作為辨識程序的輸入

%Outputs: sysd=離散傳函，sys=連續傳函，err=error

n=length(U);

Y=reshape(Y,n,1);?%輸出和輸入的維數肯定是要對齊的，無須多言

U=reshape(U,n,1);

theta=[0.1;0.1;0.1;0.1;0.1]?%注意這是一個5*1矩陣

P=2^25*eye(5);?%P，和theta都是遞推參數初始值，每一步遞推都會改變

R0=1;?%就是一個賦值

for?m=3:n?%遞推過程

X=[Y(m-1)

Y(m-2) U(m) U(m-1)

U(m-2)]';?%X就是遞推最小二乘法的算法公式中的Φ(k)

alfa=1/(R0+X'*P*X);?%計算因子

L=alfa*P*X;

theta(:,m-1)=theta(:,m-2)+L*(Y(m)-X'*theta(:,m-2));?%theta=θ?%每迭代一次，theta就增加一列，新增加的那一列就是最新的參數估值。

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);

輸入以下代碼即可：

Y=data(:,3);U=data(:,2);

[sysd,sys,e]=ID(Y,U,0.001);

總結

以上是生活随笔為你收集整理的函数的递推matlab,关于递推最小二乘法辨识参数的matlab编程（含注释）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。