基于子空间辨识的PEMFC电特性状态空间模型

clear all close all %% 输入输出数据 load u1_1000; load u2_1000; N=1000; % 采样数 I=u1; % 负载电流 f_h2=u2; % 氢气流量 u=[u1,u2]; [Vout,Pout]=pemfc_model(f_h2,I,N); y=[Vout,Pout]; u=u'; y=y'; m=2; l=2; %% Hankel矩阵 i=5; % Hankel矩阵行数 j=N-2*i+1; % Hankel矩阵列数 % Up for cnt1=1:i for cnt2=1:j Up1(:,cnt2)=u(:,cnt1+cnt2-1); end Up(2*cnt1-1:2*cnt1,:)=Up1; end % Yp for cnt1=1:i for cnt2=1:j Yp1(:,cnt2)=y(:,cnt1+cnt2-1); end Yp(2*cnt1-1:2*cnt1,:)=Yp1; end % Uf for cnt1=1:i for cnt2=1:j Uf1(:,cnt2)=u(:,cnt1+cnt2+i-1); end Uf(2*cnt1-1:2*cnt1,:)=Uf1; end % Yf for cnt1=1:i for cnt2=1:j Yf1(:,cnt2)=y(:,cnt1+cnt2+i-1); end Yf(2*cnt1-1:2*cnt1,:)=Yf1; end % Upa+ for cnt1=1:i+1 for cnt2=1:j Up1(:,cnt2)=u(:,cnt1+cnt2-1); end Upa(2*cnt1-1:2*cnt1,:)=Up1; end % Ypa+ for cnt1=1:i+1 for cnt2=1:j Yp1(:,cnt2)=y(:,cnt1+cnt2-1); end Ypa(2*cnt1-1:2*cnt1,:)=Yp1; end % Ufs- for cnt1=1:i-1 for cnt2=1:j Uf1(:,cnt2)=u(:,cnt1+cnt2+i); end Ufs(2*cnt1-1:2*cnt1,:)=Uf1; end % Yfs- for cnt1=1:i-1 for cnt2=1:j Yf1(:,cnt2)=y(:,cnt1+cnt2+i); end Yfs(2*cnt1-1:2*cnt1,:)=Yf1; end Wp=[Up;Yp]; Wpa=[Upa;Ypa]; Ob=(Yf-Yf*Uf'*pinv(Uf*Uf')*Uf)*pinv(Wp-Wp*Uf'*pinv(Uf*Uf')*Uf)*Wp; Obs=(Yfs-Yfs*Ufs'*pinv(Ufs*Ufs')*Ufs)*pinv(Wpa-Wpa*Ufs'*pinv(Ufs*Ufs')*Ufs)*Wpa; %% SVD分解 W1=eye(l*i); W2=eye(j)-Uf'*pinv(Uf*Uf')*Uf; [U,S,V]=svd(W1*Ob*W2); ss = diag(S); % figure(1); % bar(ss); % 绘制直方图 n=2; S1=S(1:n,1:n); U1=U(:,1:n); %% 求A B C D T=pinv(W1)*U1*S1^0.5; a=size(T,1); Ts=T(1:a-l,:); X=pinv(T)*Ob; Xa=pinv(Ts)*Obs; U=u(:,1:j); Y=y(:,1:j); M=[X;U]; Ma=[Xa;Y]; E=Ma*pinv(M); [sr,sc]=size(E); A=E(1:n,1:n); B=E(1:n,n+1:sc); C=E(n+1:sr,1:n); D=E(n+1:sr,n+1:sc); %% 模型验证 eplot=zeros(441,1); cnt=1; cnta=zeros(10,1); cntb=zeros(10,1); k=1; % e=0 u=u'; y=y'; for cnt1=-3.92:0.001:-3.9 for cnt2=-2.91:0.001:-2.89 x0=[cnt1,cnt2]; yi=dlsim(A,B,C,D,u,x0); e=sumsqr(yi-y); if e<=229.633 cnta(k)=cnt1; cntb(k)=cnt2; k=k+1; end eplot(cnt)=e; cnt=cnt+1; end end x0=[-3.905,-2.89]; yi=dlsim(A,B,C,D,u,x0); figure(2) plot(y,'r') hold on plot(yi,'b') title('subspace identification coefficient matrix'); xlabel('time'); ylabel('y'); % figure(3) % plot(eplot)