function [Up3,up3] = problem3(N,Q,R,x_measuringP3,d0,alpha)

%% System Model
% alpha = 0.3;
% number of disturbances
nd = 2;
% d0 = [0;0];
% N = 8; 
nfull = 16;
% number of inputs
nu = 2;
[Gfull] = model_baltrunc(nfull);
% make a discrete-time approximation
Ts	= 1;
Gfull = c2d(Gfull,Ts);
[A B C Bd] = deal(Gfull.a,Gfull.b(:,1:nu),Gfull.c,Gfull.b(:,nu+1:end));
% number of states
nx = size(Gfull.a, 2);
x0 = x_measuringP3;
% weight matrices
% Q = 1*(diag(ones(16,1),0));
% R = diag(ones(2,1),0);
umax = 1;
umin = -1;
% X,Y,Z,QQ,RR
[X,Y,Z]=xyfun(alpha,N,A,B,Bd);
[QQ,RR]=QRfun(Q,R,N);

% Matrixes for quadprog
Fx = Z'*QQ*Y;
Fd = X'*QQ*Y;
H = (Y'*QQ*Y)+RR;
G = [eye(nu*N);-1*eye(nu*N)];
W = [umax*ones(nu*N,1);-umin*ones(nu*N,1)];

% Define variables 
    % variables for full-order contoller
    U 	= sdpvar(nu*N,1);   % inputs
    lambda 	= sdpvar(length(W),1);

% KKT conditions and objective value for full-order problem
slack 	= W - G*U;

KKT 	= [H*U + Fx'*x0 + Fd'*d0 + G'*lambda == 0,
           G*U <= W,
           lambda >= 0];

       
s 	= binvar(length(slack),1);
M  	= 100;	% 'big-M'
KKT 	= [KKT, -M*s <= lambda <= M*s];
KKT 	= [KKT, -M*(1-s) <= slack <= M*(1-s)];
KKT 	= [KKT, -M <= lambda <= M];
       
F = [KKT, umin <=U<= umax];

% formulation of MPC problem
% objective function of MPC problem
objective_function = U'*H*U+d0'*Fd*U+x0'*Fx*U;
u = [];
u1 = [];
u2 = [];

%% Solving the MPC problem
Options = sdpsettings('verbose', 0);
   
    sol = solvesdp(F, [], Options);
        if sol.problem == 1
            disp('MPC problem is not solvable');
        end


%   optimal value   
    Uopt = double(U(:));
    u0opt = double(U(1:2));
    u1 = [u1; u0opt(1)];
    u2 = [u2; u0opt(2)];
   
    Up3 = Uopt;
    up3 = u0opt;
    