资料来自于网络,如有侵权请联系删除!
控制对象为:7.5e^(150s)/(120s+1); 其余参数见仿真程序;
clc
clear
G=tf(1.2,[208 1],'inputdelay',80); %input('输入传递函数G=')
%设置DMC参数
Ts=5; %input('采样周期Ts='); %采样时间
P=50; %input('预测周期P='); %预测步长
M=2; %input('控制时域M='); %控制步长
umax=500;%input('控制量最大值umax='); %控制量最大值
%delta_umax=input('控制量变化量最大值delata_max=')%空置量变化最大值
N=200; %截断步长
Yrk1=zeros(1,P);
c=3; %%%设定参考值;阶跃值
yr0=1000;
startvalue=0; %系统初始输出值
x1=startvalue;
x2=0;
for n=1:P
x2=(0.9^n)*x1+(1-0.9^n)*c; %轨迹参考参数a=0.8
x1=x2;
Yrk1(n)=x2;
end
%建立系统的阶跃响应模型
[y0,x0]=step(G,0:5:5*N);
%初始化
A=zeros(P,M);
for i=1:N
a(i)=y0(i);
end
for i=1:P
for j=1:M
if i-j+1>0
A(i,j)=a(i-j+1); %构造矩阵A
end
end
end
%初始化向量ys,y,u,e和矩阵A0
ys=ones(N,1);
y=zeros(N,1);
u=zeros(N,1);
e=zeros(N,1);
A0=zeros(P,N-1);
for i=1:P
for j=N-2:-1:1
if N-j+1+i-1<=N
A0(i,j)=a(N-j+1+i-1)-a(N-j+i-1); %构造矩阵A0
else
A0(i,j)=0;
end
end
A0(i,N-1)=a(i+1);
end
%DMC程序
for k=2:N
Uk_1=zeros(N-1,1);
for i=1:N-1
if k-N+i<=0
Uk_1(i)=0;
else
Uk_1(i)=u(k-N+i);
end
end
%Y0为Ym
Y0=A0*Uk_1;
e(k)=y(k-1)-Y0(1);
if k==50
Y0=Y0+2;
end
Yr=zeros(P,1);
%Yr为期望输出
for i=1:P
Yr(i)=Yrk1(i); %yr为期望的输出
end
Ek=zeros(P:1);
for i=1:P
Ek(i)=e(k);
end
%inv----------求逆矩阵
A2=zeros(P:1);
for i=1:P
A2(i)=1-0.9^i;
end
% A3=zeros(P:1);
A3=A2*(yr0-y(k-1));
delta_u=inv(A'*A+eye(M))*A'*(Yr-Y0); %控制增量计算
% for i=1:M
% if delta_u(i)>delta_umax
% delta_u(i)=delta_umax;
% end
% end
for i=1:M
if k+i-1<=N
u(k+i-1)=u(k+i-1-1)+delta_u(i); %控制率的计算
if u(k+i-1)>umax
u(k+i-1)=umax;
end
end
end
temp=0;%设置在k-j+1时刻以点的控制率
for j=1:N-1
if k-j<=0
temp;
else
if k-j-1<=0
temp=temp+a(j)*u(k-j);
else
temp=temp+a(j)*(u(k-j)-u(k-j-1));
end
end
end
if k-N<=0
y(k)=temp+e(N);
else
y(k)=temp+a(N)*u(k-N)+e(N);
end
end
%画图
t=5*(1:N);
subplot(211);
plot(t,y);
title('DMC控制输出曲线');
xlabel('t');
ylabel('y');
grid on
subplot(212);
plot(t,u,'r');
title('控制作用');
xlabel('t');
ylabel('u');
grid on