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;