function CON = intersect_old(CONfeas,CONadd,clean_up_flag)% Compute the intersection of two linear constraint sets%% Syntax:%   "C = intersect(a,b,clean_up_flag)"%% Description:%   "intersect(a,b,clean_up_flag)" returns a linearcon object containing%   the intersection of "a" and "b".  If no intersection exists,%   "intersect(a,b,clean_up_flag)" returns an empty linear constraint%   object.  The value of "clean_up_flag" determines whether or not%   redundant constraints are removed from the returned linear constraint%   object;  "clean_up_flag=1"--remove all redundant constraints,%   "clean_up_flag=0"--do nothing.%% Note:%   At most one equality constraint is allowed in "a", and if "b" has an%   equality constraint, it must be the same as the one present in "a".%% See Also:%   linearcon,isfeasible,andglobal GLOBAL_APPROX_PARAMglobal GLOBAL_OPTIM_PARepsilon = GLOBAL_APPROX_PARAM.poly_epsilon;hyperplane_tol = GLOBAL_APPROX_PARAM.poly_hyperplane_tol;CON = linearcon;if isempty(CONfeas) | isempty(CONadd)  returnendif (length(CONfeas.dE) > 1) | (length(CONadd.dE) > 1)  fprintf('\007intersect: Invalid constraints given, more than 1 equality found\n')  returnendif (length(CONfeas.dE) == 0) & (length(CONadd.dE) == 1)  fprintf('\007intersect: Invalid constraints given, additional equality constraints\n')  returnend% If an equality constraint is found in both sets of constraints,% check if they're the same constraint, if not return an empty% cell arrayif (length(CONfeas.dE) == 1) & (length(CONadd.dE) == 1)  MATRIX = [CONfeas.CE CONfeas.dE            CONadd.CE  CONadd.dE];  if rank(MATRIX,hyperplane_tol) > 1%    fprintf('intersect: Patches w/ different eq constraints found\n')    return  else    CONadd.CE = []; CONadd.dE = [];  endend  % Start with the feasible constraintsCE = CONfeas.CE; dE = CONfeas.dE;CI = CONfeas.CI; dI = CONfeas.dI;% Find out if each new inequality constraint is feasibleCIadd = CONadd.CI; dIadd = CONadd.dI;for k = 1:size(CIadd,1)  cIk = CIadd(k,:);  dIk = dIadd(k);  if is_new_ineq(CONfeas,cIk,dIk)    xmin = linprog(cIk',CI,dI,CE,dE,[],[],[],GLOBAL_OPTIM_PAR); %  xmin = lp(cIk',[CE; CI],[dE; dI],[],[],[],size(CE,1));    fmin = cIk*xmin;    isfeasible = (fmin < dIk-epsilon);    if isfeasible      CI = [CI; cIk];      dI = [dI; dIk];    else      return    end  endendCON = linearcon(CE,dE,CI,dI);if clean_up_flag  CON = clean_up(CON);endreturn% ----------------------------------------------------------------------------function new = is_new_ineq(CONfeas,c,d)% Check if the given constraint cTx <= d% is already in the feasible constraint CONfeasglobal GLOBAL_APPROX_PARAMhyperplane_tol = GLOBAL_APPROX_PARAM.poly_hyperplane_tol;% First search the equality constraintsCE = CONfeas.CE; dE = CONfeas.dE;new = 1;for k = 1:length(dE)  MATRIX = [CE(k,:) dE(k); c d];  if (rank(MATRIX,hyperplane_tol) < 2)    new = 0;    break;  endendif new  % Now search the inequality constraints  CI = CONfeas.CI; dI = CONfeas.dI;  for k = 1:length(dI)    MATRIX = [CI(k,:) dI(k); c d];    if (rank(MATRIX,hyperplane_tol) < 2) & (CI(k,:)*c' > 0)      new = 0;      break;    end  endendreturn