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')