www.pudn.com > dmc.rar > dmc.m, change:2013-03-21,size:1268b


ntot=250; dt=0.5; t=dt*[0:ntot-1]';% time vector 
r = [1,3,-1,-2.5,3.5,5.4,-4,-6,5,2.2];% random square wave setpoint 
Yspt = ones(size(t(1:fix(length(t)/length(r))+1)))*r(:)'; 
Yspt = Yspt(1:length(t))'; % Cut off excess data points 
Kp=1.5;taup=2.5; tau = 2.0; zeta = 0.45; % process model 
Gpc = tf(Kp*[-taup 1],[tau^2 2*tau*zeta 1]); 
Gp = c2d(Gpc,dt); % discrete plant 
[Phi,Del,C,D] = ssdata(Gp); % state space for x0 
Hc = 8; %Control horizon,Hc 
Hp = 25;%Prediction horizon,Hp> Hc 
[g,ts] = step(Gp,dt*Hp); % step response coefficients 
g(1) = [];  
G = flipud(hankel(flipud(g))); G(:,Hc+1:end) = []; 
weight=1; %control move weightingw°›0 
GG = G'*G;  
invG = (GG + weight*eye(size(GG)))\G'; 
x = [0,0]'; y=C*x; u=0; p=0; % cheat here but Ok 
Uopt=zeros(size(t)); Y=zeros(size(t)); % initialise u(t) 
for i=1:ntot-Hp-1; 
du = invG*(Yspt(i:i+Hp-1)-p); % analytical optimisation 
u = du(1)+u ; % only take first change 
x = Phi*x + Del*u; y = C*x; % apply 1 control move 
% Calculate future free response (assume all future du=0) 
[p,px] = dlsim(Phi,Del,C,D,u*ones(Hp+1,1),x); 
p(1)=[]; % kill repeat of x0 
Y(i,:) = y; Uopt(i,:) = u; 
end % for i 
Yspt(1) = []; Yspt=[Yspt;1]; 
subplot(2,1,1);plot(t,Y,t,Yspt,':') % plot results 
subplot(2,1,2);stairs(t,Uopt);