clear all;
close all;
a=cell(1,2);
b=cell(1,2);
c=cell(1,1);
d=cell(1,1);
syms k;
k=length(k);
if (0<=k<=150)
a=[1 0.9234];b=[7.2402 0.9485];c=1;d=1;
elseif (150<k<=300)
a=[1 0.8981]; b=[9.9901 0.14142]; c=1; d=1;
elseif (300<k<=450)
a=[1 0.8838]; b=[9.6041 0.34067]; c=1; d=1;
else (450<k<=600)
a=[1 0.9234]; b=[7.2402 0.9485]; c=1; d=1;
end
%na,nb is the order of
na=length(a)-1;
b=[zeros(1,d-1) b];
nb =length(b)-1;
aa=conv(a,[1 -1]);naa=na+1;
N1=d;N=15;Nu=5;% 最小输出长度、预测长度、控制长度
gamma=1*eye(Nu);alpha=0.11;% 控制加权矩阵、输出柔化系数
L=600;%控制步数
uk=zeros(d+nb,1);% 输入初值: uk(i) 表示 u(k-i)
duk=zeros(d+nb,1);% 控制增量初值
yk=zeros(naa,1);% 输出初值
w=10*[ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1)];
% 设定值
xi=sqrt(0.01)*randn(L,1);%白噪声序列
% 求解多步 Diophantine 方程并构建 F1、 F2、G
[E,F,G]=multidiophantine(aa,b,c,N);
G=G(N1: N, : );
F1=zeros(N-N1+1,Nu); F2=zeros(N-N1+1,nb);
for i=1:N-N1+1
for j=1:min(i,Nu); F1(i,j)=F(i+N1-1,i+N1-1-j+1);end
for j=1:nb;
F2(i,j)=F(i+N1-1,i+N1-1+j);end
end
for k=1:L
if (1<=k<=150)
time(k)=k;
a=[1 0.9234]; b=[7.2402 0.9485]; c=1; d=1;
y(k)=-aa(2:naa+1)*yk+b*duk(1:nb+1)+xi(k);%
Yk=[y(k);yk(1:na)];%构建向量 Y(k)
dUk=duk(1:nb);% 构建向量△ U(k-j)
elseif (150<k<=300)
time(k)=k;
a=[1 0.8981]; b=[9.9901 0.14142]; c=1; d=1;
y(k)=-aa(2:naa+1)*yk+b*duk(1:nb+1)+xi(k);%
Yk=[y(k);yk(1:na)];%采集输出数据构建向量 Y(k)
dUk=duk(1:nb);% 构建向量△ U(k-j)
elseif (300<k<=450)
time(k)=k;
a=[1 0.8838]; b=[9.6041 0.34067]; c=1; d=1;
y(k)=-aa(2:naa+1)*yk+b*duk(1:nb+1)+xi(k);%
Yk=[y(k);yk(1:na)];%采集输出数据构建向量 Y(k)
dUk=duk(1:nb);% 构建向量△ U(k-j)
else (450<k<=L)
time(k)=k;
a=[1 0.9234]; b=[7.2402 0.9485]; c=1; d=1;
y(k)=-aa(2:naa+1)*yk+b*duk(1:nb+1)+xi(k);%
Yk=[y(k);yk(1:na)];%构建向量 Y(k)
dUk=duk(1:nb);% 构建向量△ U(k-j)
end
yr(k)=y(k);
for i=1:N
yr(k+i)=alpha*yr(k+i-1)+(1-alpha)*w(k+d);
end
Yr=[yr(k+N1:k+N)]';%构建向量 Yk(k)
% 求控制量
dU=inv(F1'*F1+gamma)*F1'*(Yr-F2*dUk-G*Yk); %
du(k)=dU(1); u(k)=uk(1)+du(k);
% 更新数据
for i=1+nb:-1:2
uk(i)=uk(i-1);
duk(i)=duk(i-1);
end
uk(1)=u(k);
duk(1)=du(k);
for i=naa:-1:2
yk(i)=yk(i-1);
end
yk(1)=y(k);
end
subplot(2,1,1);
plot(time,w(1:L),'m:',time,y);
xlabel('k');ylabel('w(k),y(k)');
legend('w(k)','y(k)');
subplot(2,1,2);
plot(time,u);
xlabel('k');ylabel('u(k)');
评论0