Main Content

Investigate Linear Infeasibilities

Open Live Script

This example shows how to investigate the linear constraints that cause a problem to be infeasible. For further details about these techniques, see Chinneck [1] and [2].

If linear constraints cause a problem to be infeasible, you might want to find a subset of the constraints that is infeasible, but removing any member of the subset makes the rest of the subset feasible. The name for such a subset is Irreducible Infeasible Subset of Constraints, abbreviated IIS. Conversely, you might want to find a maximum cardinality subset of constraints that is feasible. This subset is called a Maximum Feasible Subset, abbreviated MaxFS. The two concepts are related, but not identical. A problem can have many different IISs, some with different cardinality.

This example shows two ways of finding an IIS, and two ways of obtaining a feasible set of constraints. The example assumes that all given bounds are correct, meaning the lb and ub arguments have no errors.

Infeasible Example

Create a random matrix A representing linear inequalities of size 150-by-15. Set the corresponding vector b to a vector with entries of 10, and change 5% of those values to –10.

N = 15;
rng default
A = randn([10*N,N]);
b = 10*ones(size(A,1),1);
Aeq = [];
beq = [];
b(rand(size(b)) <= 0.05) = -10;
f = ones(N,1);
lb = -f;
ub = f;

Check that problem is infeasible.

[x,fval,exitflag,output,lambda] = linprog(f,A,b,Aeq,beq,lb,ub);
No feasible solution found.

Linprog stopped because no point satisfies the constraints.

Deletion Filter

To identify an IIS, perform the following steps. Given a set of linear constraints numbered 1 through N, where all problem constraints are infeasible:

For each i from 1 to N:

  • Temporarily remove constraint i from the problem.

  • Test the resulting problem for feasibility.

  • If the problem is feasible without constraint i, return constraint i to the problem.

  • If the problem is not feasible without constraint i, do not return constraint i to the problem

Continue to the next i (up to value N).

At the end of this procedure, the constraints that remain in the problem form an IIS.

For MATLAB® code that implements this procedure, see the deletionfilter helper function at the end of this example.

Note: If you use the live script file for this example, the deletionfilter function is already included at the end of the file. Otherwise, you need to create this function at the end of your .m file or add it as a file on the MATLAB path. The same is true for the other helper functions used later in this example.

See the effect of deletionfilter on the example data.

[ineqs,eqs,ncall] = deletionfilter(A,b,Aeq,beq,lb,ub);

The problem has no equality constraints. Find the indices for the inequality constraints and the value of b(iis).

iis = find(ineqs)
iis = 
114

b(iis)
ans = 
-10

Only one inequality constraint causes the problem to be infeasible, along with the bound constraints. The constraint is

A(iis,:)*x <= b(iis).

Why is this constraint infeasible together with the bounds? Find the sum of the absolute values of that row of A.

disp(sum(abs(A(iis,:))))
    8.4864

Due to the bounds, the x vector has values between –1 and 1, and so A(iis,:)*x cannot be less than b(iis) = –10.

How many linprog calls did deletionfilter perform?

disp(ncall)
   150

The problem has 150 linear constraints, so the function called linprog 150 times.

Elastic Filter

As an alternative to the deletion filter, which examines every constraint, try the elastic filter. This filter works as follows.

First, allow each constraint i to be violated by a positive amount e(i), where equality constraints have both additive and subtractive positive elastic values.

Aineqx≤bineq+eAeqx=beq+e1-e2

Next, solve the associated linear programming problem (LP)

minx,e∑ei

subject to the listed constraints and with ei≥0.

  • If the associated LP has a solution, remove all constraints that have a strictly positive associated ei, and record those constraints in a list of indices (potential IIS members). Return to the previous step to solve the new, reduced associated LP.

  • If the associated LP has no solution (is infeasible) or has no strictly positive associated ei, exit the filter.

The elastic filter can exit in many fewer iterations than the deletion filter, because it can bring many indices at once into the IIS, and can halt without going through the entire list of indices. However, the problem has more variables than the original problem, and its resulting list of indices can be larger than an IIS. To find an IIS after running an elastic filter, run the deletion filter on the result.

For MATLAB® code that implements this filter, see the elasticfilter helper function at the end of this example.

See the effect of elasticfilter on the example data.

[ineqselastic,eqselastic,ncall] = ...
    elasticfilter(A,b,Aeq,beq,lb,ub);

The problem has no equality constraints. Find the indices for the inequality constraints.

iiselastic = find(ineqselastic)
iiselastic = 5×1

     2
    60
    82
    97
   114

The elastic IIS lists five constraints, whereas the deletion filter found only one. Run the deletion filter on the returned set to find a genuine IIS.

A_1 = A(ineqselastic > 0,:);
b_1 = b(ineqselastic > 0);
[dineq_iis,deq_iis,ncall2] = deletionfilter(A_1,b_1,Aeq,beq,lb,ub);
iiselasticdeletion = find(dineq_iis)
iiselasticdeletion = 
5

The fifth constraint in the elastic filter result, inequality 114, is the IIS. This result agrees with the answer from the deletion filter. The difference between the approaches is that the combined elastic and deletion filter approach uses many fewer linprog calls. Display the total number of linprog calls used by the elastic filter followed by the deletion filter.

disp(ncall + ncall2)
     7

Remove IIS in a Loop

Generally, obtaining a single IIS does not enable you to find all the reasons that your optimization problem fails. To correct an infeasible problem, you can repeatedly find an IIS and remove it from the problem until the problem becomes feasible.

The following code shows how to remove one IIS at a time from a problem until the problem becomes feasible. The code uses an indexing technique to keep track of constraints in terms of their positions in the original problem, before the algorithm removes any constraints.

The code keeps track of the original variables in the problem by using a Boolean vector activeA to represent the current constraints (rows) of the A matrix, and a Boolean vector activeAeq to represent the current constraints of the Aeq matrix. When adding or removing constraints, the code indexes into A or Aeq so that the row numbers do not change, even though the number of constraints changes.

Running this code returns idx2, a vector of the indices of the nonzero elements in activeA:

idx2 = find(activeA)

Suppose that var is a Boolean vector that has the same length as idx2. Then

idx2(find(var))

expresses var as indices into the original problem variables. In this way, the indexing can take a subset of a subset of constraints, work with only the smaller subset, and still unambiguously refer to the original problem variables.

opts = optimoptions('linprog','Display',"none");
activeA = true(size(b));
activeAeq = true(size(beq));
[~,~,exitflag] = linprog(f,A,b,Aeq,beq,lb,ub,opts);
ncl = 1;
while exitflag < 0
    [ineqselastic,eqselastic,ncall] = ...
        elasticfilter(A(activeA,:),b(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub);
    ncl = ncl + ncall;
    idxaa = find(activeA);
    idxae = find(activeAeq);
    tmpa = idxaa(find(ineqselastic));
    tmpae = idxae(find(eqselastic));
    AA = A(tmpa,:);
    bb = b(tmpa);
    AE = Aeq(tmpae,:);
    be = beq(tmpae);
    [ineqs,eqs,ncall] = ...
        deletionfilter(AA,bb,AE,be,lb,ub);
    ncl = ncl + ncall;
    activeA(tmpa(ineqs)) = false;
    activeAeq(tmpae(eqs)) = false;
    disp(['Removed inequalities ',int2str((tmpa(ineqs))'),' and equalities ',int2str((tmpae(eqs))')])
    [~,~,exitflag] = ...
        linprog(f,A(activeA,:),b(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub,opts);
    ncl = ncl + 1;
end
Removed inequalities 114 and equalities 
Removed inequalities 97 and equalities 
Removed inequalities 64  82 and equalities 
Removed inequalities 60 and equalities 

fprintf('Number of linprog calls: %d\n',ncl)
Number of linprog calls: 28

Notice that the loop removes inequalities 64 and 82 simultaneously, which indicates that these two constraints form an IIS.

Find MaxFS

Another approach for obtaining a feasible set of constraints is to find a MaxFS directly. As Chinneck [1] explains, finding a MaxFS is an NP-complete problem, meaning the problem does not necessarily have efficient algorithms for finding a MaxFS. However, Chinneck proposes some algorithms that can work efficiently.

Use Chinneck's Algorithm 7.3 to find a cover set of constraints that, when removed, gives a feasible set. The algorithm is implemented in the generatecover helper function at the end of this example.

[coversetineq,coverseteq,nlp] = generatecover(A,b,Aeq,beq,lb,ub)
coversetineq = 5×1

   114
    97
    60
    82
     2


coverseteq =

     []

nlp = 
40

Remove these constraints and solve the LP.

usemeineq = true(size(b));
usemeineq(coversetineq) = false; % Remove inequality constraints
usemeeq = true(size(beq));
usemeeq(coverseteq) = false; % Remove equality constraints
[xs,fvals,exitflags] = ...
    linprog(f,A(usemeineq,:),b(usemeineq),Aeq(usemeeq),beq(usemeeq),lb,ub);
Optimal solution found.

Notice that the cover set is exactly the same as the iiselastic set from Elastic Filter. In general, the elastic filter finds too large a cover set. Chinneck's Algorithm 7.3 starts with the elastic filter result and then retains only the constraints that are necessary.

Chinneck's Algorithm 7.3 takes 40 calls to linprog to complete the calculation of a MaxFS. This number is a bit more than 28 calls used earlier in the process of deleting IIS in a loop.

Also, notice that the inequalities removed in the loop are not exactly the same as the inequalities removed by Algorithm 7.3. The loop removes inequalities 114, 97, 82, 60, and 64, while Algorithm 7.3 removes inequalities 114, 97, 82, 60, and 2. Check that inequalities 82 and 64 form an IIS (as indicated in Remove IIS in a Loop), and that inequalities 82 and 2 also form an IIS.

usemeineq = false(size(b));
usemeineq([82,64]) = true;
ineqs = deletionfilter(A(usemeineq,:),b(usemeineq),Aeq,beq,lb,ub);
disp(ineqs)
   1
   1

usemeineq = false(size(b));
usemeineq([82,2]) = true;
ineqs = deletionfilter(A(usemeineq,:),b(usemeineq),Aeq,beq,lb,ub);
disp(ineqs)
   1
   1

References

[1] Chinneck, J. W. Feasibility and Infeasibility in Optimization: Algorithms and Computational Methods. Springer, 2008.

[2] Chinneck, J. W. "Feasibility and Infeasibility in Optimization." Tutorial for CP-AI-OR-07, Brussels, Belgium. Available at https://www.sce.carleton.ca/faculty/chinneck/docs/CPAIOR07InfeasibilityTutorial.pdf.

Helper Functions

This code creates the deletionfilter helper function.

function [ineq_iis,eq_iis,ncalls] = deletionfilter(Aineq,bineq,Aeq,beq,lb,ub)
ncalls = 0;
[mi,n] = size(Aineq); % Number of variables and linear inequality constraints
f = zeros(1,n);
me = size(Aeq,1); % Number of linear equality constraints
opts = optimoptions("linprog","Algorithm","dual-simplex","Display","none");

ineq_iis = true(mi,1); % Start with all inequalities in the problem
eq_iis = true(me,1); % Start with all equalities in the problem

for i=1:mi
    ineq_iis(i) = 0; % Remove inequality i
    [~,~,exitflag] = linprog(f,Aineq(ineq_iis,:),bineq(ineq_iis),...
                                Aeq,beq,lb,ub,opts);
    ncalls = ncalls + 1;
    if exitflag == 1 % If now feasible
        ineq_iis(i) = 1; % Return i to the problem
    end
end
for i=1:me
    eq_iis(i) = 0; % Remove equality i
    [~,~,exitflag] = linprog(f,Aineq,bineq,...
                                Aeq(eq_iis,:),beq(eq_iis),lb,ub,opts);
    ncalls = ncalls + 1;
    if exitflag == 1 % If now feasible
        eq_iis(i) = 1; % Return i to the problem
    end
end
end

This code creates the elasticfilter helper function.

function [ineq_iis,eq_iis,ncalls,fval0] = elasticfilter(Aineq,bineq,Aeq,beq,lb,ub)
ncalls = 0;
[mi,n] = size(Aineq); % Number of variables and linear inequality constraints
me = size(Aeq,1);
Aineq_r = [Aineq -1.0*eye(mi) zeros(mi,2*me)];
Aeq_r = [Aeq zeros(me,mi) eye(me) -1.0*eye(me)]; % Two slacks for each equality constraint
lb_r = [lb(:); zeros(mi+2*me,1)];
ub_r = [ub(:); inf(mi+2*me,1)];
ineq_slack_offset = n;
eq_pos_slack_offset = n + mi;
eq_neg_slack_offset = n + mi + me;
f = [zeros(1,n) ones(1,mi+2*me)];
opts = optimoptions("linprog","Algorithm","dual-simplex","Display","none");
tol = 1e-10;

ineq_iis = false(mi,1);
eq_iis = false(me,1);
[x,fval,exitflag] = linprog(f,Aineq_r,bineq,Aeq_r,beq,lb_r,ub_r,opts);
fval0 = fval;
ncalls = ncalls + 1;
while exitflag == 1 && fval > tol % Feasible and some slacks are nonzero
    c = 0;
    for i = 1:mi        
        j = ineq_slack_offset+i;
        if x(j) > tol
            ub_r(j) = 0.0;
            ineq_iis(i) = true;
            c = c+1;
        end
    end
    for i = 1:me
        j = eq_pos_slack_offset+i;
        if x(j) > tol            
            ub_r(j) = 0.0;
            eq_iis(i) = true;
            c = c+1;
        end
    end
    for i = 1:me
        j = eq_neg_slack_offset+i;
        if x(j) > tol            
            ub_r(j) = 0.0;
            eq_iis(i) = true;
            c = c+1;
        end
    end
    [x,fval,exitflag] = linprog(f,Aineq_r,bineq,Aeq_r,beq,lb_r,ub_r,opts);
    if fval > 0
        fval0 = fval;
    end
    ncalls = ncalls + 1;
end 
end

This code creates the generatecover helper function. The code uses the same indexing technique for keeping track of constraints as the Remove IIS in a Loop code.

function [coversetineq,coverseteq,nlp] = generatecover(Aineq,bineq,Aeq,beq,lb,ub)
% Returns the cover set of linear inequalities, the cover set of linear
% equalities, and the total number of calls to linprog.
% Adapted from Chinneck [1] Algorithm 7.3. Step numbers are from this book.
coversetineq = [];
coverseteq = [];
activeA = true(size(bineq));
activeAeq = true(size(beq));
% Step 1 of Algorithm 7.3
[ineq_iis,eq_iis,ncalls] = elasticfilter(Aineq,bineq,Aeq,beq,lb,ub);
nlp = ncalls;
ninf = sum(ineq_iis(:)) + sum(eq_iis(:));
if ninf == 1
    coversetineq = ineq_iis;
    coverseteq = eq_iis;
    return
end
holdsetineq = find(ineq_iis);
holdseteq = find(eq_iis);
candidateineq = holdsetineq;
candidateeq = holdseteq;
% Step 2 of Algorithm 7.3
while sum(candidateineq(:)) + sum(candidateeq(:)) > 0
    minsinf = inf;
    ineqflag = 0;
    for i = 1:length(candidateineq(:))
        activeA(candidateineq(i)) = false;
        idx2 = find(activeA);
        idx2eq = find(activeAeq);
        [ineq_iis,eq_iis,ncalls,fval] = elasticfilter(Aineq(activeA,:),bineq(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub);
        nlp = nlp + ncalls;
        ineq_iis = idx2(find(ineq_iis));
        eq_iis = idx2eq(find(eq_iis));
        if fval == 0
            coversetineq = [coversetineq;candidateineq(i)];
            return
        end
        if fval < minsinf
            ineqflag = 1;
            winner = candidateineq(i);
            minsinf = fval;
            holdsetineq = ineq_iis;
            if numel(ineq_iis(:)) + numel(eq_iis(:)) == 1
                nextwinner = ineq_iis;
                nextwinner2 = eq_iis;
                nextwinner = [nextwinnner,nextwinner2];
            else
                nextwinner = [];
            end
        end
        activeA(candidateineq(i)) = true;
    end
    for i = 1:length(candidateeq(:))
        activeAeq(candidateeq(i)) = false;
        idx2 = find(activeA);
        idx2eq = find(activeAeq);
        [ineq_iis,eq_iis,ncalls,fval] = elasticfilter(Aineq(activeA,:),bineq(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub);
        nlp = nlp + ncalls;
        ineq_iis = idx2(find(ineq_iis));
        eq_iis = idx2eq(find(eq_iis));
        if fval == 0
            coverseteq = [coverseteq;candidateeq(i)];
            return
        end
        if fval < minsinf
            ineqflag = -1;
            winner = candidateeq(i);
            minsinf = fval;
            holdseteq = eq_iis;
            if numel(ineq_iis(:)) + numel(eq_iis(:)) == 1
                nextwinner = ineq_iis;
                nextwinner2 = eq_iis;
                nextwinner = [nextwinnner,nextwinner2];
            else
                nextwinner = [];
            end
        end
        activeAeq(candidateeq(i)) = true;
    end
% Step 3 of Algorithm 7.3
    if ineqflag == 1
        coversetineq = [coversetineq;winner];
        activeA(winner) = false;
        if nextwinner
            coversetineq = [coversetineq;nextwinner];
            return
        end
    end
    if ineqflag == -1
        coverseteq = [coverseteq;winner];
        activeAeq(winner) = false;
        if nextwinner
            coverseteq = [coverseteq;nextwinner];
            return
        end
    end
    candidateineq = holdsetineq;
    candidateeq = holdseteq;
end
end

See Also

Related Topics