最近在做关于过热气温的动态建模问题。有现场运行的历史数据,要找出导前区和惰性区的传递函数。
对这类算法了解不多,程序读起来比较吃力,所以就转来一篇完整的辨识程序,在原有基础上进行了简化,并稍加注解一下,也能让有共同需要的朋友少花点时间。
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;%就是一个赋值
form=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);%计算实际输出值和估值输出的误差
ifabs(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);