matlab编的DMC程序
clear all;
% close all;
%系统模型建立
num=[0.8];
den=[225 1];
[a,b,c,d]=tf2ss(num,den);
% step(num,den);
Ts=30;
lambda=60;
[ad,bd,cd,dd]=c2dt(a,b,c,Ts,lambda);
[numd,dend]=ss2tf(ad,bd,cd,dd);
[a,x]=dstep(ad,bd,cd,dd);
P=10;
M=5;
N=50;
%动态矩阵A
for i=1:P
for j=1:M
if j>i
A(i,j)=0;
else
A(i,j)=a(i-j+1);
end
end
end
%A0
g(1)=a(1);
for i=2:N
g(i)=a(i)-a(i-1);
end
for i=1:P
for j=1:N-1;
if i>j
A0(i,j)=0;
elseif j==N-1
A0(i,j)=a(i+1);
else
A0(i,j)=g(i-j+N);
end
end
end
%权矩阵Q、R
Q=0.5*eye(P);
R=0.01*eye(M);
h(1)=0.5;
for i=2:P
h(i)=0.5;
end
for i=1:M
if i==1
D1(i)=1;
else
D1(i)=0;
end
end
D=D1*(inv(A'*Q*A+R)*A'*Q);
y_1=0;
u_3=0;
%初始化
for i=1:N-1
u(i)=0;
end
for i=1:M
du(i)=0;
end
for k=1:100
t(k)=k*Ts;
for i=1:P
r(i)=10;
end
y(k)=0.9149*y_1+0.0680*u_3;
%预测控制
ym=0;
for i=1:N-1
ym=ym+g(i)*u(N-i);
end
e=y(k)-ym;
dU=D*(r'-(A0*u'+h'*e));
U(k)=u(N-1)+dU;
%加扰动
if k==50
U(k)=U(k)+5;
end
for i=1:N-2
u(i)=u(i+1);
end
u(N-1)=U(k);
u_3=u(N-3);
% end
y_1=y(k);
EE(k)=r(1)-y(k);
end
figure(1)
% subplot(211)
plot(t,y,'r')
hold on;
grid on;
% subplot(212)
% plot(t,U,'g');
% hold on
% grid on;