www.pudn.com > guangyiyucekongzhi.rar > test.m, change:2009-02-08,size:2794b

```clear all;close all;%程序开始，清空工作间关闭窗口

s=tf('s');
Ts=12;%Ts为采样周期(步长）
T_final = 125*Ts; %T_final 仿真时间
g11=0.162*(275*s+1)/((168*s+1)^2*(11.5*s+1));
figure(1);
step(g11);
%g12=-0.081*(0.01+0.99/(97*s+1));
%g21=2.116*(457*s+1)/((221*s+1)^2*(21.8*s+1));
%g22=1.483*s*(150*s+1)/((632*s^2+40*s+1)*(2.7*s+1));
%G=[g11,g12;g21,g22];
%step(G);
G11=c2d(g11,Ts);
[num,den]=tfdata(G11,'v');

%initialization
na = length(den)-1;
nb = length(num)-1;
k = 1;%从第k 步开始预测
st = 1;%系统设定值setpoint在初始时读入
p = 20; %根据需要选取预测步长P
M = 8; %控制步长为M
aifa = 0.05 ; %输入柔化因子aifa
lambda = 10; %输入增量参数lambda
b = 1;%wug's的阶梯因子
h=[];

%the first line of F and G
F(1,1:na) = den(1,1:na)-den(1,2:na+1);
F(1,na+1) = den(1,na+1);%F（1，：）(q^-1)=(1-a1)+(a1-a2)*q^-1+...+ana*q^-na
E(1,1)=1;
G_ef(1,:) = num;

%solve the diophantus equation
for j= 2:k+p-1
%get G
G_ef(j,1:j-1) = G_ef(j-1,1:j-1);
for i= j:nb+j-1
G_ef(j,i) = G_ef(j-1,i) + F(j-1,1)*num(i-j+1);
end
G_ef(j,nb+j) = F(j-1,1) * num(nb+1);
%get F
for i= 1:na
F(j,i) = F(j-1,i+1) + F(j-1,1)*F(1,i);
end
F(j,na+1) = F(j-1,1) * den(na+1);
end
%*-------------finished solving the diophantus equation-----------*

%calculate wug's G2
for i = 1:p%构造讲义中的 G
for j = 1:i
G(i,j) = G_ef(i,i-j+1);
end
end

G1 = G(1:p,1:M);%根据预测控制时域得到的G1

for i=1:M
%h=[h;b^(i-1)];
h=[h;b^(i-1)];
end
G2=G1*h;%引入阶梯因子后的G
% end calculate G2

for i= 1:p  %构造计算y^1（即G_poly*H+F*yk )时所需的矩阵G_poly
G_poly(i,1:nb+k-1) = G_ef(k-1+i,i+1:nb+k+i-1);
end

ypast(:,1) = zeros(na+1,1); %初始值ypast(na+1项)，从y(t-1)开始至y(t-na-1).
Upast(:,1) = zeros(nb+k+1,1); %预测前已知的控制量u(nb+k)，从u(t-1)开始至u(t-nb-k),deltu(t-nb-k+1)时用到u(t-nb-k)
outU = [];
outY = [];
for I = 0:Ts:T_final %开始循环计算
delt_ypast = ypast(1:na,1)-ypast(2:na+1,1);
deltU = Upast(1:nb+k,1) - Upast(2:nb+k+1,1);%产生已知输入控制信号deltu(t)
yt = ypast(1,1) - den(1,2:na+1) * delt_ypast(1:na,1) + num(1,:) * deltU(k:nb+k,1);
outY = [outY;yt];
%采样，得到仿真对象的当前输出值yt
ypast = [yt;ypast];
ypast = ypast(1:na+1,1);
for i = 1:p
ypre1(i,1) = F(k+i-1,:) * ypast(:,1) + G_poly(i,:) * deltU(1:nb+k-1,1);
end
%以上为求y^1,ypre1为

%Tender the setpoint track vector
W(1,1) = yt;
for i = 2:p+1
W(i,1)=aifa*W(i-1,1)+(1-aifa)*st;
end
deltU1 = inv(G2'*G2)*G2'*(W(2:p+1,1)-ypre1(:,1)) ;%当前控制律deltU1
ut = Upast(1,1)+deltU1(1,1);
Upast = [ut;Upast];
Upast = Upast(1:nb+k+1,1);
outU = [outU;ut];
end
t=[0:Ts:T_final];
figure(2);
plot(t,outY),hold on,
figure(3);
plot(t,outU,'r')```