carima.m

基于受控自回归积分滑动平均 (CARIMA)模型控对象的预测控制程序

%-构造CARIMA模型。
A=[1 -0.7];
B=[0.9 -0.6];
a=A(2:end);
b=B;

%-控制参数。
N1=3;
Nu=3;
lanbt=0.5;

na=length(a);
nb=length(b)-1;

%-初始化向量。
aa(1)=a(1)-1;
if na>=2
   for i=2:na
     aa(i)=a(i) -a(i-1);
   end
end
aa(na+1)=-a(na);

F=zeros(N1,na+1);
F(1,:)=-aa;
e=zeros(1,N1);
e(1)=1;
f=zeros(1,na+1);
H=zeros(N1,nb);
H(1,:)=b(2:end);
g=zeros(1,N1);
g(1)=e(1)*b(1);

%-求解diophantine方程。
for j=1:N1-1

    for i=1:na
       F(j+1,i)=F(j,i+1)-aa(i)*F(j,1);
    end
    F(j+1,na+1) =-aa(na+1)*F(j,1); 
    e(j+1)=F(j,1);
    if nb>1
      for i=1:nb-1
         H(j+1,i)=e(j+1)*b(i+1)+H(j,i+1);
      end
    end
      H(j+1,nb)=e(j+1)*b(nb+1);
      g(j+1)=e(j+1)*b(1)+H(j,1);
end

%-控制量求解。
G=zeros(N1,Nu);
for i=1:Nu
   G(i:end,i)=g(1:N1-i+1) ;
end
Q=inv(G'*G+lanbt*diag(ones(Nu,1)))*G';
w=10;
Arfar=0.1;
Tstop=16;
yr=zeros(Tstop+N1,1);
yr(1)=0.11;y(2)=4.21;y(3)=6.31;
yr1=zeros(N1,1);
delta_u=zeros(1,Tstop);
u=zeros(1,Tstop);
y=zeros(1,Tstop);
y(1)=0.01;y(2)=0.12;y(3)=0.13;
for t=max(na,nb+1)+1:Tstop
    for j=1:N1
     yr(t+j)=Arfar*yr(t+j-1)+(1-Arfar)*w;
      yr1(j)=yr(t+j);
    end
    
    y(t)=0.7*y(t-1)+0.9*u(t-1)-0.6*u(t-2);
    p=Q(1,:);
    deltau=p*[yr1-H*delta_u(t-1)];
    f=p(1)*F(1,:);
    for i=2:N1
       f=f+ p(i)*F(i,:);
    end
    delta_u(t)=deltau-f(1)*y(t);
    for i=2:length(f)
      delta_u(t)=delta_u(t)-f(i)*y(t-i+1);
    end
    u(t)=u(t-1)+delta_u(t);
  
end
plot(y);
grid 
hold on
plot(yr,'r');