GPC11_广义预测控制代码_广义预测控制_模型预测控制_广义预测_gpc_

%close all %clear all %Parameter %Mode %mode=1 tao does not change, while a changes %mode=2 tao changes , while a does not change Mode=1; %System structure na=2; nb=1; %predictive control parameter N=5; %predictive time domain Nu=2; %control time domain %identification parameters po_thita=zeros(na+nb+1,1); mu=0.97; fai=zeros(na+nb+1,1); P=10^6*eye(na+nb+1); %set trace parameters T=1; ta=1.8; alpha=exp(-T/ta); c=2; %input and output thita1=zeros(5,1000); y_temp=zeros(1,na+1); y=zeros(1,1000); u=zeros(1,1000); u(1,4)=1; noise=0.01*randn(1,991); %noise=random('Normal', 0,0.01,1,991); %Start for k=10:1000 switch Mode case 1 tao=6; if k<=201 a=0.8854; elseif k<=301 a=0.6725; elseif k<=401 a=0.5968; elseif k<=501 a=0.7762; else a=0.8854; end case 2 a=0.8854; if k<=201 tao=6; elseif k<=301 tao=4; elseif k<=401 tao=8; elseif k<=501 tao=6; else tao=6; end end y_temp(1)=(a+0.264)*y_temp(2)-0.264*a*y_temp(3)+0.864*u(k-tao)-0.2731*u(k-tao-1); y(k)=y_temp(1)+noise(k-9); y_temp(3)=y_temp(2); y_temp(2)=y_temp(1); %identification %fai for i=1:na fai(i)=-(y(k-i)-y(k-i-1)); end for i=1:nb+1 fai(na+i)=u(k-tao+1-i)-u(k-tao+1-i-1); end kai=P*fai*inv(fai'*P*fai+mu); P=(mu^(-1))*(eye(na+nb+1)-kai*fai')*P; po_thita=po_thita+kai*(y(k)-y(k-1)-fai'*po_thita); thita=[1,po_thita']; thita1(:,k)=thita; %diophantine f=zeros(na+1,1); f(1)=1; f_f=zeros(N,na+1); E=cell(1,N); %元胞数组将不同类型的数组，这里e中元素长短不一 E{1}=1; g=zeros(1,N); G=zeros(N,Nu); f_t=zeros(N,1); g_t=zeros(N,1); w=zeros(N,1); for ite1=1:N f_e=f(1); for ite2=1:na f(ite2)=f(ite2+1)-(thita(ite2+1)-thita(ite2))*f_e; end f(na+1)=thita(na+1)*f_e; f_f(ite1,:)=f'; if ite1<N E{ite1+1}=[E{ite1},f(1)]; end end %compute g aa=zeros(1,nb+1); bb=cell(1,N); for i=1:nb+1 aa(i)=thita(na+1+i); end for i1=1:N bb{i1}=conv(aa,E{i1}); end for i=1:N g(i)=bb{N}(i); end %compute g_g，也就是G for i1=1:Nu for i2=1:N if i2<i1 G(i2,i1)=0; else G(i2,i1)=g(i2-i1+1); end end end %compute f_n*y(k) for i1=1:N for i2=1:na+1 f_t(i1)=f_t(i1)+f_f(i1,i2)*y(k+1-i2); end end %compute g_n*u(k) i2的作用考虑到两者相减剩余的项数 for i1=1:N for i2=1:nb g_t(i1)=g_t(i1)+bb{i1}(i1+i2)*(u(k-tao-i2)-u(k-tao-i2-1)); end end % set trace for i1=1:N w(i1)=alpha^(i1)*y(k)+(1-alpha^(i1))*c; end %compute u r_ft=g_t+f_t; be_u=inv(G'*G+eye(Nu))*G'; u(k-tao+1)=u(k-tao)+be_u(1,:)*(w-r_ft); end %plot y u figure plot(y); figure; plot(u);