function out = project(in,vars)% Compute the projection of the input linearcon object (a polyhedron)% onto the specified set of variables %% Syntax:%   "C = project(con,v)"%% Description:%   "project(con,v)" returns a linear constraint object representing%   the projection of "con" onto the variables given in the indices%   array "v".%% Examples:%   Suppose the input linear constraint object "in" has 5 variables.%%%%%   "out = project(in,[1 4 3])"%%%%   returns "out" a linear constraint object representing a projection of%   the polyhedron onto the variables [x1 x4 x3]' (in the order specified%   in the indices array). %% Note:%   The input linear constraint object must be bounded for this function%   to work properly.%% See Also:%   linearcon[CE,dE,CI,dI] = linearcon_data(in);if length(unique(vars)) ~= length(vars)  error('Variable indices must be unique')endN = size(CI,2);if any(vars > N) | any(vars < 1)  error('Invalid variable indices given')endif size(CE,1) > 1  error('Invalid constraints, more than one equality found')end% Reorder the variables so that the variables specified in vars appear% first.othervars = setdiff([1:N],vars);if ~isempty(CE)  CE = [CE(:,vars) CE(:,othervars)];endCI = [CI(:,vars) CI(:,othervars)];% Keep eliminating the last variable until we have only the variables in% vars in the constraints.while size(CI,2) > length(vars)  if isempty(CE)    [CI,dI] = elim_no_eq(CI,dI);  else    [CE,dE,CI,dI] = elim_one_eq(CE,dE,CI,dI);  endendout = linearcon(CE,dE,CI,dI);return% -----------------------------------------------------------------------------function [Celim,delim] = elim_no_eq(C,d)iplus = [];iminus = [];izero = [];N = size(C,2);% Sort indicies into zero, positive, or negative coefficient of% variable to be eliminatedfor k = 1:size(C,1)  if (C(k,N) == 0)    izero = [izero k];  else     if (C(k,N) > 0)       iplus = [iplus k];    else      iminus = [iminus k];    end  endendCelim = []; delim = [];% If coefficient of last variable to be eliminated is zero,% then remove the last column and retain the row as isfor k = izero  Celim = [Celim; C(k,1:N-1)];  delim = [delim; d(k)];end% Do <name/type of> elimination on rows with nonzero coefficientsfor k = iplus  cplus = C(k,1:N-1);  cNplus = C(k,N);   dplus = d(k);  for l = iminus    cminus = C(l,1:N-1);    cNminus = C(l,N);     dminus = d(l);    ckl = cplus/cNplus - cminus/cNminus;    dkl = dplus/cNplus - dminus/cNminus;    magnitude = sqrt(ckl*ckl');    if (magnitude == 0) % What happens if magnitude is zero and dkl is positive???      if (dkl < 0)        error('Infeasible constraint found! Something is wrong!')      end    else      ckl = ckl/magnitude;      dkl = dkl/magnitude;      Celim = [Celim; ckl];      delim = [delim; dkl];    end  endendreturn% -----------------------------------------------------------------------------function [Ceelim,deelim,Cielim,dielim] = elim_one_eq(Ce,de,Ci,di)N = size(Ce,2);Ceelim = []; deelim = [];Cielim = []; dielim = [];if (Ce(1,N) == 0)  % If coefficient of variable to be eliminated is zero, then remove last column  % and retain the equality as is.  Perform elimination as if with no equality.  Ceelim = Ce(:,1:N-1);  deelim = de;  [Cielim,dielim] = elim_no_eq(Ci,di);  returnelse  % Otherwise, solve equality for variable to be eliminated in terms of the other  % variables and substitute into the inequalities  Cehat = Ce(1,1:N-1);  CeN = Ce(1,N);  for k = 1:size(Ci,1)    cik = Ci(k,1:N-1) - Ci(k,N)*Cehat/CeN;    dik = di(k) - Ci(k,N)*de/CeN;    magnitude = sqrt(cik*cik');    if magnitude == 0      if dik < 0        error('Infeasible constraint found! Something is wrong!')      end    else      cik = cik/magnitude;      dik = dik/magnitude;      Cielim = [Cielim; cik];      dielim = [dielim; dik];    end  endendreturn