### Highlights from bnb

• ... NLCONST Helper function to find the constrained minimum of a function
• BNB20(fun,x0,xstat,xl,xu,...BNB20 Finds the constrained minimum of a function of several possibly integer variables.
• BNBGUI(file)
• BNBGUICB(action,file); BNBGUICB Callback function for BNBGUI 2.0.
• [X,lambda,exitflag,output...QP Quadratic programming subproblem. Handles qp and constrained
• guierr() This is the machine-generated representation of a Handle Graphics object
• guifun() This is the machine-generated representation of a Handle Graphics object
• guimain() This is the machine-generated representation of a Handle Graphics object
• guiset() This is the machine-generated representation of a Handle Graphics object
• guiupd() This is the machine-generated representation of a Handle Graphics object
• View all files

# bnb

### Koert Kuipers (view profile)

14 Feb 2000 (Updated )

BNB20 solves mixed integer nonlinear optimization problems

[X,lambda,exitflag,output,how]=qpsub(H,f,A,B,lb,ub,X,neqcstr,verbosity,caller,ncstr,numberOfVariables,options)
```function [X,lambda,exitflag,output,how]=qpsub(H,f,A,B,lb,ub,X,neqcstr,verbosity,caller,ncstr,numberOfVariables,options)
%QP Quadratic programming subproblem. Handles qp and constrained
%   linear least-squares as well as subproblems generated from NLCONST.
%
%   X=QP(H,f,A,b) solves the quadratic programming problem:
%
%            min 0.5*x'Hx + f'x   subject to:  Ax <= b
%             x
%

%   Copyright (c) 1990-98 by The MathWorks, Inc.
%   \$Revision: 1.21 \$  \$Date: 1998/09/01 21:37:56 \$
%   Andy Grace 7-9-90. Mary Ann Branch 9-30-96.

% Define constant strings
global maxSQPiter;
NewtonStep = 'Newton';
SteepDescent = 'steepest descent';
Conls = 'lsqlin';
Lp = 'linprog';
Qpsub = 'qpsub';
how = 'ok';
exitflag = 1;
output = [];
iterations = 0;
if nargin < 13
options = [];
end

lb=lb(:); ub = ub(:);

msg = nargchk(12,12,nargin);
if isempty(verbosity), verbosity = 1; end
if isempty(neqcstr), neqcstr = 0; end

LLS = 0;
if strcmp(caller, Conls)
LLS = 1;
[rowH,colH]=size(H);
numberOfVariables = colH;
end
if strcmp(caller, Qpsub)
normalize = -1;
else
normalize = 1;
end

simplex_iter = 0;
if  norm(H,'inf')==0 | isempty(H), is_qp=0; else, is_qp=1; end

if LLS==1
is_qp=0;
end

normf = 1;
if normalize > 0
% Check for lp
if ~is_qp & ~LLS
normf = norm(f);
if normf > 0
f = f./normf;
end
end
end

% Handle bounds as linear constraints
arglb = ~eq(lb,-inf);
lenlb=length(lb); % maybe less than numberOfVariables due to old code
if nnz(arglb) > 0
lbmatrix = -eye(lenlb,numberOfVariables);

A=[A; lbmatrix(arglb,1:numberOfVariables)]; % select non-Inf bounds
B=[B;-lb(arglb)];
end

argub = ~eq(ub,inf);
lenub=length(ub);
if nnz(argub) > 0
ubmatrix = eye(lenub,numberOfVariables);
A=[A; ubmatrix(argub,1:numberOfVariables)];
B=[B; ub(argub)];
end
ncstr=ncstr + nnz(arglb) + nnz(argub);

% Figure out max iteration count
% For now, don't limit the iterations when qpsub is called from nlconst
if exist('maxSQPiter','var'), maxSQPiters = maxSQPiter; else maxSQPiters=inf; end;
maxiter = optimget(options,'MaxIter',maxSQPiters);

% Used for determining threshold for whether a direction will violate
% a constraint.
normA = ones(ncstr,1);
if normalize > 0
for i=1:ncstr
n = norm(A(i,:));
if (n ~= 0)
A(i,:) = A(i,:)/n;
B(i) = B(i)/n;
normA(i,1) = n;
end
end
else
normA = ones(ncstr,1);
end
errnorm = 0.01*sqrt(eps);

tolDep = 100*numberOfVariables*eps;
lambda=zeros(ncstr,1);
aix=lambda;
ACTCNT=0;
ACTSET=[];
ACTIND=0;
CIND=1;
eqix = 1:neqcstr;

%------------EQUALITY CONSTRAINTS---------------------------
Q = zeros(numberOfVariables,numberOfVariables);
R = [];
indepInd = 1:ncstr;

if neqcstr>0
% call equality constraint solver
[Q,R,A,B,CIND,X,Z,actlambda,how,...
ACTSET,ACTIND,ACTCNT,aix,eqix,neqcstr,ncstr,remove,exitflag]= ...
eqnsolv(A,B,eqix,neqcstr,ncstr,numberOfVariables,LLS,H,X,f,normf,normA,verbosity, ...
aix,how,exitflag);

if ~isempty(remove)
indepInd(remove)=[];
normA = normA(indepInd);
end

if ACTCNT >= numberOfVariables - 1
simplex_iter = 1;
end
[m,n]=size(ACTSET);

if strcmp(how,'infeasible')
% Equalities are inconsistent, so X and lambda have no valid values
% Return original X and zeros for lambda.
output.iterations = iterations;
return
end

err = 0;
if neqcstr > numberOfVariables
err = max(abs(A(eqix,:)*X-B(eqix)));
if (err > 1e-8)  % Equalities not met
how='infeasible';
% was exitflag = 7;
exitflag = -1;
if verbosity > 0
disp('Exiting: The equality constraints are overly stringent;')
disp('         there is no feasible solution.')
end
% Equalities are inconsistent, X and lambda have no valid values
% Return original X and zeros for lambda.
output.iterations = iterations;
return
else % Check inequalities
if (max(A*X-B) > 1e-8)
how = 'infeasible';
% was exitflag = 8;
exitflag = -1;
if verbosity > 0
disp('Exiting: The constraints or bounds are overly stringent;')
disp('         there is no feasible solution.')
disp('         Equality constraints have been met.')
end
end
end
if is_qp
actlambda = -R\(Q'*(H*X+f));
elseif LLS
actlambda = -R\(Q'*(H'*(H*X-f)));
else
actlambda = -R\(Q'*f);
end
lambda(indepInd(eqix)) = normf * (actlambda ./normA(eqix));
output.iterations = iterations;
return
end
if isempty(Z)
if is_qp
actlambda = -R\(Q'*(H*X+f));
elseif LLS
actlambda = -R\(Q'*(H'*(H*X-f)));
else
actlambda = -R\(Q'*f);
end
lambda(indepInd(eqix)) = normf * (actlambda./normA(eqix));
if (max(A*X-B) > 1e-8)
how = 'infeasible';
% was exitflag = 8;
exitflag = -1;
if verbosity > 0
disp('Exiting: The constraints or bounds are overly stringent;')
disp('         there is no feasible solution.')
disp('         Equality constraints have been met.')
end
end
output.iterations = iterations;
return
end

% Check whether in Phase 1 of feasibility point finding.
if (verbosity == -2)
cstr = A*X-B;
mc=max(cstr(neqcstr+1:ncstr));
if (mc > 0)
X(numberOfVariables) = mc + 1;
end
end
else
Z=1;
end

% Find Initial Feasible Solution
cstr = A*X-B;
mc=max(cstr(neqcstr+1:ncstr));
if mc>eps
A2=[[A;zeros(1,numberOfVariables)],[zeros(neqcstr,1);-ones(ncstr+1-neqcstr,1)]];
quiet = -2;
[XS,lambdaS,exitflagS,outputS] = qpsub([],[zeros(numberOfVariables,1);1],A2,[B;1e-5], ...
[],[],[X;mc+1],neqcstr,quiet,Qpsub,size(A2,1),numberOfVariables+1);
X=XS(1:numberOfVariables);
cstr=A*X-B;
if XS(numberOfVariables+1)>eps
if XS(numberOfVariables+1)>1e-8
how='infeasible';
% was exitflag = 4;
exitflag = -1;
if verbosity > 0
disp('Exiting: The constraints are overly stringent;')
disp('         no feasible starting point found.')
end
else
how = 'overly constrained';
% was exitflag = 3;
exitflag = -1;
if verbosity > 0
disp('Exiting: The constraints are overly stringent;')
disp(' initial feasible point found violates constraints ')
disp(' by more than eps.');
end
end
lambda(indepInd) = normf * (lambdaS((1:ncstr)')./normA);
output.iterations = iterations;
return
end
end

if (is_qp)
gf=H*X+f;
%  SD=-Z*((Z'*H*Z)\(Z'*gf));
[SD, dirType] = compdir(Z,H,gf,numberOfVariables,f);

% Check for -ve definite problems:
%  if SD'*gf>0, is_qp = 0; SD=-SD; end
elseif (LLS)
HXf=H*X-f;
gf=H'*(HXf);
HZ= H*Z;
[mm,nn]=size(HZ);
if mm >= nn
%   SD =-Z*((HZ'*HZ)\(Z'*gf));
[QHZ, RHZ] =  qr(HZ);
Pd = QHZ'*HXf;
% Now need to check which is dependent
if min(size(RHZ))==1 % Make sure RHZ isn't a vector
depInd = find( abs(RHZ(1,1)) < tolDep);
else
depInd = find( abs(diag(RHZ)) < tolDep );
end
end
if mm >= nn & isempty(depInd) % Newton step
SD = - Z*(RHZ(1:nn, 1:nn) \ Pd(1:nn,:));
dirType = NewtonStep;
else % steepest descent direction
SD = -Z*(Z'*gf);
dirType = SteepDescent;
end
else % lp
gf = f;
SD=-Z*Z'*gf;
dirType = SteepDescent;
if norm(SD) < 1e-10 & neqcstr
% This happens when equality constraint is perpendicular
% to objective function f.x.
actlambda = -R\(Q'*(gf));
lambda(indepInd(eqix)) = normf * (actlambda ./ normA(eqix));
output.iterations = iterations;
return;
end
end

oldind = 0;

% The maximum number of iterations for a simplex type method is when ncstr >=n:
% maxiters = prod(1:ncstr)/(prod(1:numberOfVariables)*prod(1:max(1,ncstr-numberOfVariables)));

%--------------Main Routine-------------------
while 1 & iterations <= maxiter
iterations = iterations + 1;
if isinf(verbosity)
curr_out = sprintf('Iter: %5.0f, Active: %5.0f, step: %s, proc: %s',iterations,ACTCNT,dirType,how);
disp(curr_out);
end

% Find distance we can move in search direction SD before a
% constraint is violated.
% Gradient with respect to search direction.
GSD=A*SD;

% Note: we consider only constraints whose gradients are greater
% than some threshold. If we considered all gradients greater than
% zero then it might be possible to add a constraint which would lead to
% a singular (rank deficient) working set. The gradient (GSD) of such
% a constraint in the direction of search would be very close to zero.
indf = find((GSD > errnorm * norm(SD))  &  ~aix);

if isempty(indf) % No constraints to hit
STEPMIN=1e16;
dist=[]; ind2=[]; ind=[];
else % Find distance to the nearest constraint
dist = abs(cstr(indf)./GSD(indf));
[STEPMIN,ind2] =  min(dist);
ind2 = find(dist == STEPMIN);
% Bland's rule for anti-cycling: if there is more than one
% blocking constraint then add the one with the smallest index.
ind=indf(min(ind2));
% Non-cycling rule:
% ind = indf(ind2(1));
end
%-----Update X-------------

% Assume we do not delete a constraint
delete_constr = 0;

if ~isempty(indf)& isfinite(STEPMIN) % Hit a constraint
if strcmp(dirType, NewtonStep)
% Newton step and hit a constraint: LLS or is_qp
if STEPMIN > 1  % Overstepped minimum; reset STEPMIN
STEPMIN = 1;
delete_constr = 1;
end
X = X+STEPMIN*SD;
else
% Not a Newton step and hit a constraint: is_qp or LLS or maybe lp
X = X+STEPMIN*SD;
end
else %  isempty(indf) | ~isfinite(STEPMIN)
% did not hit a constraint
if strcmp(dirType, NewtonStep)
% Newton step and no constraint hit: LLS or maybe is_qp
STEPMIN = 1;   % Exact distance to the solution. Now delete constr.
X = X + SD;
delete_constr = 1;
else % Not a Newton step: is_qp or lp or LLS
if is_qp
% Is it semi-def, neg-def or indef?
eigoptions.disp = 0;
ZHZ = Z'*H*Z;
if numberOfVariables < 400 % only use EIGS on large problems
[VV,DD] = eig(ZHZ);
[smallRealEig, eigind] = min(diag(DD));
ev = VV(:,eigind(1));
else
[ev,smallRealEig,flag] = eigs(ZHZ,1,'sr',eigoptions);
if flag  % Call to eigs failed
[VV,DD] = eig(ZHZ);
[smallRealEig, eigind] = min(diag(DD));
ev = VV(:,eigind(1));
end
end

else % define smallRealEig for LLS
smallRealEig=0;
end

if (~is_qp & ~LLS) | (smallRealEig < -100*eps) % LP or neg def: not LLS
% neg def -- unbounded
if norm(SD) > errnorm
if normalize < 0
STEPMIN=abs((X(numberOfVariables)+1e-5)/(SD(numberOfVariables)+eps));
else
STEPMIN = 1e16;
end
X=X+STEPMIN*SD;
how='unbounded';
% was exitflag = 5;
exitflag = -1;
else % norm(SD) <= errnorm
how = 'ill posed';
% was exitflag = 6;
exitflag = -1;

end
if verbosity > 0
if norm(SD) > errnorm
disp('Exiting: The solution is unbounded and at infinity;')
disp('         the constraints are not restrictive enough.')
else
disp('Exiting: The search direction is close to zero; ')
disp('      the problem is ill-posed.')
disp('      The gradient of the objective function may be zero')
disp('         or the problem may be badly conditioned.')
end
end % if verbosity > 0
output.iterations = iterations;
return
else % singular: solve compatible system for a solution: is_qp or LLS
if is_qp
projH = Z'*H*Z;
Zgf = Z'*gf;
projSD = pinv(projH)*(-Zgf);
else % LLS
projH = HZ'*HZ;
Zgf = Z'*gf;
projSD = pinv(projH)*(-Zgf);
end

% Check if compatible
if norm(projH*projSD+Zgf) > 10*eps*(norm(projH) + norm(Zgf))
% system is incompatible --> it's a "chute": use SD from compdir
% unbounded in SD direction
if norm(SD) > errnorm
if normalize < 0
STEPMIN=abs((X(numberOfVariables)+1e-5)/(SD(numberOfVariables)+eps));
else
STEPMIN = 1e16;
end
X=X+STEPMIN*SD;
how='unbounded';
% was exitflag = 5;
exitflag = -1;
else % norm(SD) <= errnorm
how = 'ill posed';
%was exitflag = 6;
exitflag = -1;
end
if verbosity > 0
if norm(SD) > errnorm
disp('Exiting: The solution is unbounded and at infinity;')
disp('         the constraints are not restrictive enough.')
else
disp('Exiting: The search direction is close to zero; ')
disp('      the problem is ill-posed.')
disp('      The gradient of the objective function may be zero')
disp('         or the problem may be badly conditioned.')
end
end % if verbosity > 0
output.iterations = iterations;
return
else % Convex -- move to the minimum (compatible system)
SD = Z*projSD;
dirType = 'singular';
% First check if constraint is violated.
GSD=A*SD;
indf = find((GSD > errnorm * norm(SD))  &  ~aix);
if isempty(indf) % No constraints to hit
STEPMIN=1;
delete_constr = 1;
dist=[]; ind2=[]; ind=[];
else % Find distance to the nearest constraint
dist = abs(cstr(indf)./GSD(indf));
[STEPMIN,ind2] =  min(dist);
ind2 = find(dist == STEPMIN);
% Bland's rule for anti-cycling: if there is more than one
% blocking constraint then add the one with the smallest index.
ind=indf(min(ind2));
end
if STEPMIN > 1  % Overstepped minimum; reset STEPMIN
STEPMIN = 1;
delete_constr = 1;
end
X = X + STEPMIN*SD;
end
end % if ~is_qp | smallRealEig < -eps
end % if strcmp(dirType, NewtonStep)
end % if ~isempty(indf)& isfinite(STEPMIN) % Hit a constraint

if delete_constr
% Note: only reach here if a minimum in the current subspace found
if ACTCNT>0
if ACTCNT>=numberOfVariables-1,
% Avoid case when CIND is greater than ACTCNT
if CIND <= ACTCNT
ACTSET(CIND,:)=[];
ACTIND(CIND)=[];
end
end
if is_qp
rlambda = -R\(Q'*(H*X+f));
elseif LLS
rlambda = -R\(Q'*(H'*(H*X-f)));
% else: lp does not reach this point
end
actlambda = rlambda;
actlambda(eqix) = abs(rlambda(eqix));
indlam = find(actlambda < 0);
if (~length(indlam))
lambda(indepInd(ACTIND)) = normf * (rlambda./normA(ACTIND));
output.iterations = iterations;
return
end
% Remove constraint
lind = find(ACTIND == min(ACTIND(indlam)));
lind=lind(1);
ACTSET(lind,:) = [];
aix(ACTIND(lind)) = 0;
[Q,R]=qrdelete(Q,R,lind);
ACTIND(lind) = [];
ACTCNT = ACTCNT - 2;
simplex_iter = 0;
ind = 0;
else % ACTCNT == 0
output.iterations = iterations;
return
end
delete_constr = 0;
end

% Calculate gradient w.r.t objective at this point
if is_qp
gf=H*X+f;
elseif LLS % LLS
gf=H'*(H*X-f);
% else gf=f still true.
end

% Update X and calculate constraints
cstr = A*X-B;
cstr(eqix) = abs(cstr(eqix));
% Check no constraint is violated
if normalize < 0
if X(numberOfVariables,1) < eps
output.iterations = iterations;
return;
end
end

if max(cstr) > 1e5 * errnorm
if max(cstr) > norm(X) * errnorm
if ( verbosity > 0 ) & ( exitflag == 1 )
disp('Note: The problem is badly conditioned;')
disp('         the solution may not be reliable')
% verbosity = 0;
end
how='unreliable';
% exitflag = 2;
exitflag = -1;
if 0
X=X-STEPMIN*SD;
output.iterations = iterations;
return
end
end
end

if ind % Hit a constraint
aix(ind)=1;
ACTSET(CIND,:)=A(ind,:);
ACTIND(CIND)=ind;
[m,n]=size(ACTSET);
[Q,R] = qrinsert(Q,R,CIND,A(ind,:)');
end
if oldind
aix(oldind) = 0;
end
if ~simplex_iter
% Z = null(ACTSET);
[m,n]=size(ACTSET);
Z = Q(:,m+1:n);
ACTCNT=ACTCNT+1;
if ACTCNT == numberOfVariables - 1, simplex_iter = 1; end
CIND=ACTCNT+1;
oldind = 0;
else
rlambda = -R\(Q'*gf);

if isinf(rlambda(1)) & rlambda(1) < 0
fprintf('         Working set is singular; results may still be reliable.\n');
[m,n] = size(ACTSET);
rlambda = -(ACTSET + sqrt(eps)*randn(m,n))'\gf;
end
actlambda = rlambda;
actlambda(eqix)=abs(actlambda(eqix));
indlam = find(actlambda<0);
if length(indlam)
if STEPMIN > errnorm
% If there is no chance of cycling then pick the constraint
% which causes the biggest reduction in the cost function.
% i.e the constraint with the most negative Lagrangian
% multiplier. Since the constraints are normalized this may
% result in less iterations.
[minl,CIND] = min(actlambda);
else
% Bland's rule for anti-cycling: if there is more than one
% negative Lagrangian multiplier then delete the constraint
% with the smallest index in the active set.
CIND = find(ACTIND == min(ACTIND(indlam)));
end

[Q,R]=qrdelete(Q,R,CIND);
Z = Q(:,numberOfVariables);
oldind = ACTIND(CIND);
else
lambda(indepInd(ACTIND))= normf * (rlambda./normA(ACTIND));
output.iterations = iterations;
return
end
end %if ACTCNT<numberOfVariables

if (is_qp)
Zgf = Z'*gf;
if ~isempty(Zgf) & (norm(Zgf) < 1e-15)
SD = zeros(numberOfVariables,1);
else
[SD, dirType] = compdir(Z,H,gf,numberOfVariables,f);
end
elseif (LLS)
Zgf = Z'*gf;
HZ = H*Z;
if (norm(Zgf) < 1e-15)
SD = zeros(numberOfVariables,1);
else
HXf=H*X-f;
gf=H'*(HXf);
[mm,nn]=size(HZ);
if mm >= nn
[QHZ, RHZ] =  qr(HZ);
Pd = QHZ'*HXf;
% SD = - Z*(RHZ(1:nn, 1:nn) \ Pd(1:nn,:));
% Now need to check which is dependent
if min(size(RHZ))==1 % Make sure RHZ isn't a vector
depInd = find( abs(RHZ(1,1)) < tolDep);
else
depInd = find( abs(diag(RHZ)) < tolDep );
end
end
if mm >= nn & isempty(depInd) % Newton step
SD = - Z*(RHZ(1:nn, 1:nn) \ Pd(1:nn,:));
dirType = NewtonStep;
else % steepest descent direction
SD = -Z*(Z'*gf);
dirType = SteepDescent;
end
end
else % LP
if ~simplex_iter
SD = -Z*(Z'*gf);
else
SD = -Z;
else
SD = Z;
end
end
if abs(gradsd) < 1e-10  % Search direction null
% Check whether any constraints can be deleted from active set.
% rlambda = -ACTSET'\gf;
if ~oldind
rlambda = -R\(Q'*gf);
end
actlambda = rlambda;
actlambda(1:neqcstr) = abs(actlambda(1:neqcstr));
indlam = find(actlambda < errnorm);
lambda(indepInd(ACTIND)) = normf * (rlambda./normA(ACTIND));
if ~length(indlam)
output.iterations = iterations;
return
end
cindmax = length(indlam);
cindcnt = 0;
newactcnt = 0;
while (abs(gradsd) < 1e-10) & (cindcnt < cindmax)
cindcnt = cindcnt + 1;
if oldind
% Put back constraint which we deleted
[Q,R] = qrinsert(Q,R,CIND,A(oldind,:)');
else
simplex_iter = 0;
if ~newactcnt
newactcnt = ACTCNT - 1;
end
end
CIND = indlam(cindcnt);
oldind = ACTIND(CIND);

[Q,R]=qrdelete(Q,R,CIND);
[m,n]=size(ACTSET);
Z = Q(:,m:n);

if m ~= numberOfVariables
SD = -Z*Z'*gf;
else
SD = -Z;
else
SD = Z;
end
end
end
if abs(gradsd) < 1e-10  % Search direction still null
output.iterations = iterations;
return;
end
lambda = zeros(ncstr,1);
if newactcnt
ACTCNT = newactcnt;
end
end
end

if simplex_iter & oldind
% Avoid case when CIND is greater than ACTCNT
if CIND <= ACTCNT
ACTIND(CIND)=[];
ACTSET(CIND,:)=[];
CIND = numberOfVariables;
end
end
end % while 1
if iterations > maxiter
exitflag = 0;
how = 'ill-conditioned';
end

output.iterations = iterations;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function[Q,R,A,B,CIND,X,Z,actlambda,how,...
ACTSET,ACTIND,ACTCNT,aix,eqix,neqcstr,ncstr,remove,exitflag]= ...
eqnsolv(A,B,eqix,neqcstr,ncstr,numberOfVariables,LLS,H,X,f,normf,normA,verbosity, ...
aix,how,exitflag)
% EQNSOLV Helper function for QPSUB.
%    Finds a feasible point with respect to the equality constraints.
%    If the equalities are dependent but not consistent, warning
%    messages are given. If the equalities are dependent but consistent,
%    the redundant constraints are removed and the corresponding variables

% set tolerances
tolDep = 100*numberOfVariables*eps;
tolCons = 1e-10;

actlambda = [];
aix(eqix)=ones(neqcstr,1);
ACTSET=A(eqix,:);
ACTIND=eqix;
ACTCNT=neqcstr;
CIND=neqcstr+1;
Z=[]; Anew=[]; Bnew=[]; remove =[];

% See if the equalities form a consistent system:
%   QR factorization of A
[Qa,Ra,Ea]=qr(A(eqix,:));
% Now need to check which is dependent
if min(size(Ra))==1 % Make sure Ra isn't a vector
depInd = find( abs(Ra(1,1)) < tolDep);
else
depInd = find( abs(diag(Ra)) < tolDep );
end
if neqcstr > numberOfVariables
depInd = [depInd; ((numberOfVariables+1):neqcstr)'];
end

if ~isempty(depInd)
if verbosity > 0
disp('The equality constraints are dependent.')
end
how='dependent';
exitflag = 1;
bdepInd =  abs(Qa(:,depInd)'*B(eqix)) >= tolDep ;

if any( bdepInd ) % Not consistent
how='infeasible';
exitflag = 9;exitflag = -1;
if verbosity > 0
disp('The system of equality constraints is not consistent.');
if ncstr > neqcstr
disp('The inequality constraints may or may not be satisfied.');
end
disp('  There is no feasible solution.')
end
else % the equality constraints are consistent
numDepend = nnz(depInd);
% delete the redundant constraints:
% By QR factoring the transpose, we see which columns of A'
%   (rows of A) move to the end
[Qat,Rat,Eat]=qr(ACTSET');
[i,j] = find(Eat); % Eat permutes the columns of A' (rows of A)
remove = i(depInd);
if verbosity > 0
disp('The system of equality constraints is consistent. Removing');
disp('the following dependent constraints before continuing:');
disp(remove)
end
A(eqix(remove),:)=[];
B(eqix(remove))=[];
neqcstr = neqcstr - nnz(remove);
ncstr = ncstr - nnz(remove);
eqix = 1:neqcstr;
aix=[ones(neqcstr,1); zeros(ncstr-neqcstr,1)];
ACTIND = eqix;
ACTSET=A(eqix,:);

CIND = neqcstr+1;
ACTCNT = neqcstr;
end % consistency check
end % dependency check

if ~strcmp(how,'infeasible')
% Find a feasible point
if max(abs(A(eqix,:)*X-B(eqix))) > tolCons
X = A(eqix,:)\B(eqix);
end
end

[Q,R]=qr(ACTSET');
Z = Q(:,neqcstr+1:numberOfVariables);

% End of eqnsolv.m

```