function [F,P] = ddlqr(beta,Ay,Az,By,Q,R,W,Theta,Psi)
%DLQRD   Linear quadratic regulator design for discrete-time systems
%        with side constraints.
%        [F,P] = ddlqr(beta,Ay,Az,By,Q,R,W,Theta,Psi)
%	 calculates the optimal feedback gain matrix F such that 
%        the feedback law:
%
%        (*)    u = -F y  
%
%        minimizes the cost function:
%
%		J = Sum beta^n{x[n]'Qx[n] + u[n]'Ru[n] + 2x[n]'Wu[n]}
%
%        subject to the constraint equations: 
%		
%	        y[n+1] = Ay*y[n] + Az*z[n] + By*u[n]
%
%        with x=[y;z] and {z[n]} given.  The solution (*) 
%        is found by imposing the following side constraints (e.g.
%        market clearing conditions) on first-order conditions
%
%               z[n] = Theta*y[n] + Psi*u[n].
%
%	 Also returned is P, the steady-state solution to the associated 
%	 discrete matrix Riccati equation.

%	 Ellen McGrattan 1-22-89
%	 Revised 6-5-92, ERM

%  Convert to problem without cross products or discounting.

[ny,nz]  = size(Az);
n  = ny+nz;
Wy = W(1:ny,:);
Wz = W(ny+1:n,:);
ay = sqrt(beta)*(Ay-By/R*Wy');
az = sqrt(beta)*(Az-By/R*Wz');
by = sqrt(beta)*By;
q  = Q - W/R*W';
qy = q(1:ny,1:ny);
qz = q(1:ny,ny+1:n);
p  = inv(eye(nz)+Psi/R*Wz');
t  = p*(Theta-Psi/R*Wy');
p  = p*Psi;

%  Compute matrices needed for Hamiltonian

ahat = ay + az*t;
bhat = by + az*p;
qhat = qy + qz*t;
atil = ay - by/R*p'*qz';
brb  = bhat/R*by';

%  Construct Hamiltonian matrix, compute its eigenvalues and
%  eigenvectors and select vectors outside the unit circle. 

if rank(atil) == ny;
  itil  = inv(atil');
  h     = [ ahat+brb*itil*qhat, -brb*itil;
           -itil*qhat, itil];
  [v,d] = eig(h);
else;   % compute generalized eigenvalues if ahat singular
  h1    = [eye(ny),brb;zeros(ny),atil'];
  h2    = [ahat,zeros(ny);-qhat,eye(ny)];
  [v,d] = eig(h2,h1);   
end;

d      = diag(d);
[e,in] = sort(abs(d));
if (~( (e(ny) < 1) & (e(ny+1) > 1) ) )
  error('Can''t order eigenvalues');
end

W11    = v(1:ny,in(1:ny));
W21    = v(ny+1:2*ny,in(1:ny));
P      = real(W21/W11);

F      = (R+by'*P*bhat)\(by'*P*ahat);

%  The decision function for the discounted problem with cross-products:

F      = (R+Wz'*Psi)\(R*F+Wy'+Wz'*Theta);
