# Investigate Linear Infeasibilities

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.

`$\begin{array}{l}{A}_{ineq}x\le {b}_{ineq}+e\\ {A}_{eq}x={b}_{eq}+{e}_{1}-{e}_{2}\end{array}$`

Next, solve the associated linear programming problem (LP)

`$\underset{x,e}{\mathrm{min}}\sum {e}_{i}$`

subject to the listed constraints and with ${e}_{i}\ge 0$.

• If the associated LP has a solution, remove all constraints that have a strictly positive associated ${e}_{i}$, 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 ${e}_{i}$, 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```