function [X,u,mu]=steadyst(beta,A,B,D,E,param,overlap,crit,damp,maxit, ..
                            X0,u0,mu0);
%STEADYST computes the steady state values X, u, mu of the following problem:
%            
%       (P)      max sum{beta^t u(x(t),param) }
%                s.t.  x(t)=D*X(t)+E*u(t)       X(t): ns x 1,  u(t): nc x 1
%                      X(t+1)=A*X(t)+B*u(t)    
%                      X(1,t)=1  for all t
%                      mu(t)=lagrange mult. for X(t+1)=A*X(t)+B*u(t)
%                      NOTE: only mu elements for non-redundent or nontrivial
%                            constraints can be computed.  The others will 
%                            be set at -9999.
%                                                      
%  INPUTS:
%     1. parameters for problem (P): beta,A,B,D,E,param
%     2. overlap={ 0                        if no elements of u(t) coincide 
%                {                          with elements of X(t+1)
%                { (roverl x 2) matrix        otherwise,
%        where the roverl= # of times a redundancy of elements of X(t+1) and
%        u(t) occurs (i.e. X(j,t+1)=u(k,t) for some j,k at time t.)
%        The ith row of overlap is [j k], where j is the element of X(t+1) and
%        k is the element of u(t) such that X(j,t+1)=u(k,t), and this occur-
%        rence of redundency is the ith, i=1,..,roverl.
%        Example: X(t)=[1 k(t) n1(t-1)]', u(t)=[n1(t) n2(t)]' implies that:
%        overlap=[3 1].
%        ***CAUTION: if there are more examples of redundencies in your
%           system, try to set the parameters so that the constraints are
%           not redundent or trivial.  If this cannot be done, copy this 
%           program and edit the new version.***
%     3. crit=convergence criterion
%     4. damp=damping factor for slower (but more stable) convergence. damp=1
%        is full damping (no changes allowed), damp=0 is no damping.
%     5. maxit=maximum # of iterations
%     6. X0,u0,mu0: initial conditions for steady states.  These are optional
%        arguments.  If you do not have a guess for any subset of these 
%        vectors, put a -9999 in the argument.  If length(mu0)=number of 
%        states then the elements of mu0 corresponding to trivial or redundent
%        constraints are ignored. If length(mu0)=number of independent
%        constraints, all are used in computation.
%
%  OUTPUTS:
%     1. X,u,mu: steady state values of the problem (P)
%
[ns,col]=size(A);
[row,nc]=size(B);
[roverl,coverl]=size(overlap);
if X0==-9999;  X0=ones(ns,1); end;
if u0==-9999;  u0=ones(nc,1); end;
if mu0==-9999; mu0=ones(ns,1); end;
[rmu,cmu]=size(mu0);
if (rmu~=ns & rmu~=ns-roverl-1); 
  error('Wrong number of elements in initial multiplier vector.'); 
end;
X0=X0(:);
u0=u0(:);
mu0=mu0(:);
Xs0=X0;
A0=A(:,1);
As=A;
Bs=B;
C1=eye(ns);
F0=zeros(nc,roverl);
if overlap>0; 
  for i=roverl:-1:1;
    Xs0=reducmat(Xs0,overlap(i,1),0);
    A0=reducmat(A0,overlap(i,1),0);
    As=reducmat(As,overlap(i,1),overlap(i,1));
    Bs=reducmat(Bs,overlap(i,1),0);
    C1=reducmat(C1,0,overlap(i,1));
    C2(:,i)=[zeros(overlap(i,1)-1,1); 1; zeros(ns-overlap(i,1),1)];
    F0(overlap(i,2),i)=1;
    if rmu==ns; mu0=reducmat(mu0,overlap(i,1),0); end;
  end;
else;
  C2=zeros(ns,1);
end;
[rmu,cmu]=size(mu0);
if (overlap > 0 & rmu==ns-roverl); mu0=reducmat(mu0,1,0); end;
if (overlap ==0 & rmu==ns); mu0=reducmat(mu0,1,0); end;
Xs0=reducmat(Xs0,1,0);
A0=reducmat(A0,1,0);
As=reducmat(As,1,1);
Bs=reducmat(Bs,1,0);
C1=reducmat(C1,0,1);
[nss,col]=size(As);
tnn=2*nss+nc;
e=ones(tnn,1);
M20=9999999*ones(tnn,1);
%
%  Newton-Raphson is used to solve the first-order conditions of 
%     
%      max( sum{beta^t*ufn(D*(C0+C1*Xs(t)+C2*Xo(t))+E*u(t)) } )
%      st.  Xs(t+1)=As*Xs(t)+Bs*u(t)  (mult: mu),  Xo(t+1)=F0'*u(t) (mult: lam)
%
%  These conditions (in addition to the constraints) are:
%      
%  Xs:   beta*C1'*D'*uj(t) +mu(t) -beta*As'*mu(t+1) =0
%  u:    E'*uj -Bs'*mu(t) -F0*lam(t) =0
%  Xo:   beta*C2'*D'+lam(t) =0
%
for i=1:maxit;
  [us,uj,uh]=ufn(D*X0+E*u0,param); 
  M1=[eye(nss)-As           -Bs           zeros(nss);
   beta*C1'*D'*uh*D*C1   beta*C1'*D'*uh*(E+beta*D*C2*F0') eye(nss)-beta*As';
   -(E+beta*D*C2*F0')'*uh*D*C1  -(E+beta*D*C2*F0')'*uh*(E+beta*D*C2*F0') Bs'];
  M2=[A0-(eye(nss)-As)*Xs0+Bs*u0;
      -(eye(nss)-beta*As')*mu0-beta*C1'*D'*uj;
       -Bs'*mu0+(E+beta*D*C2*F0')'*uj];
  %  Are we diverging?
  if (abs(M2) > abs(M20)); 
    disp('The procedure is diverging.  Try new starting values.')
    return
  end;
  MIM=(1-damp)*(M1\M2);
  Xs1=Xs0+MIM(1:nss,1);
  u1=u0+MIM(nss+1:nss+nc,1);
  mu1=mu0+MIM(nss+nc+1:tnn,1);
  %  check for convergence:
  xum0=[Xs0;u0;mu0];
  if max(abs( ([Xs1;u1;mu1]-xum0)./max(abs(xum0),e) )) < crit; break; end;
  X0=Xs1;
  mu=mu1;
  if overlap>0;
    for j=1:roverl;
      X0=insrows(X0,overlap(j,1)-1,u1(overlap(j,2),1));
      mu=insrows(mu,overlap(j,1)-1,-9999);
    end;
  end;
  X0=insrows(X0,1,1);
  mu=insrows(mu,1,-9999);
  Xs0=Xs1;
  u0=u1;
  mu0=mu1;
  M20=M2;
end;
if i>=maxit; 
  disp('Iteration limit reached in STEADYST.'); 
else;
  sprintf('Convergence achieved in STEADYST at iteration: %g',i)
end;
X=X0;
u=u1;
