function [uh, xh, yh] = PR_omezeni(Tp, R, Q, ref, N, umin, umax, Ad,Bd,Cd,Dd)

Ns = size(Ad,1);

%predikce vystupu
V = [];
for i=0:1:Tp-1
   V=[V;Cd*Ad^i];
end

% matice S

S = zeros(Tp,Tp);
for j=1:Tp
   for i=1:Tp
      if j==i
         S(i,j)=Dd;
      elseif j==i-1
         S(i,j)=Cd*Bd;
      elseif j<i-1
         S(i,j)=Cd*Ad^(i-j-1)*Bd;
      end
   end
end

       
	% matice analytickeho regulatoru
   
   H = S'*Q*S + R;
   H = (H+H')/2;
   W = zeros(Tp);
   z = zeros(Tp, 1);
   
   Umin = (ones(1,Tp)')*umin;
   Umax = (ones(1,Tp)')*umax;   
  
	% historie

	uh = zeros(1,N);
	xh = zeros(Ns,N);
	yh = zeros(1,N);

	x = zeros(Ns,1);

%simulace
	start = Ns+1;

   for t=start:N
      ref_vyber = ref(t:t+Tp-1)';

     yx = V*x;
     j = (yx - ref_vyber)'*Q*S;
            
      %      up = quadprog(H, j, W, z', Umin, Umax);
            up = qp(H, j', W, z, Umin, Umax);

      
		y = Cd*x + Dd*up(1);
	   x = Ad*x + Bd*up(1);

   
   	xh(:,t) = x;
	   yh(t) = y;
   	uh(t) = up(1);
   
   end