g11=poly2tfd(12.8,[16.7 1],0,1);
g21=poly2tfd(6.6,[10.9 1],0,7);
g12=poly2tfd(-18.9,[21.0 1],0,3);
g22=poly2tfd(-19.4,[14.4 1],0,3);
gd1=poly2tfd(3.8,[14.9 1],0,8);
gd2=poly2tfd(4.9,[13.2 1],0,3);
tfinal=200; %Model horizon N
delt=1; %sampling period
ny=2; %Number of outputs
model=tfd2step(tfinal,delt,ny,g11,g21,g12,g22);
plant=model; %No plant/model mismatch
dmodel=[]; %Default disturbance model
dplant=tfd2step(tfinal,delt,ny,gd1,gd2)
P=90;M=30; %Horizons
ywt=[1 1]; uwt=[0.1 0.1]; %Q and R
tend=200; %final time for simulation
r=[-0.02 0.01]; %Set-point change inXD XB
a=zeros([1,tend]);
for i=50