<figure><img src="image_1.png" /></figure>

HAL
open science

# Efficient computation of solution space and conflicts detection for linear systems

## Alexandre Borthomieu

<figure><img src="image_2.png" /></figure>

### ► To cite this version:

Alexandre Borthomieu. Efficient computation of solution space and conflicts detection for linear systems. Modeling and Simulation. 2019. hal-02433237

HAL Id: hal-02433237

https://inria.hal.science/hal-02433237

Submitted on 9 Jan 2020

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.
# Efficient computation of solution space and conflicts detection for linear systems

Alexandre BORTHOMIEU

University of Grenoble

INRIA - STEEP Team

Supervised by Jean-Yves Courtonne

alexandre.borthomieu@etu.univ-grenoble-alpes.fr

## ABSTRACT

In this paper, we analyze a system of linear inequations $l \le A.x \le u$, used for the purpose of Material Flow Analysis, with three different but complementary goals: (i) given some known variables $x_i$, efficiently compute the solution space of unknown variables, (ii) if the set of constraints is infeasible, efficiently identify the conflicts, (iii) efficiently classify variables to determine whereas they are redundant, just measured, determinable or non-determinable. In each case we compare the efficiency of different algorithms or languages.

## 1. INTRODUCTION

In a context of global environmental changes which questions our modes of production and consumption, the STEEP team of INRIA aims at developing decision-making tools to promote ecological transition at local (sub-national) levels. One of these tools, called Supply Chain Material Flow Analysis (MFA) tries to provide a consistent image of a given chain (e.g., the wood supply chain) in order to stress the potential competition between different types of uses (e.g., energy vs. construction) as well as the material dependency on imports from the rest of the world. In practice, flows are the variables of the problem (some are measured, some are unknown) and the constraints are given by the law of mass conservation and by other rules linking variables (e.g., yields of transformation). The MFA follows the following steps :

1. Reconcile the data with constraint optimization, i.e., minimize the distance between measured flows (inputs) and model results (outputs) given the set of constraints. Reconciliation is necessary because the data comes from various sources which are generally inconsistent. Constraint optimization moreover “fills the gap” for missing flows. However, it provides a unique solution for each variable even when the solution is not unique: for the so-called “free” variables, other values would also respect constraints and lead to the same value of the objective function. It is thus important to identify them.

2. For this purpose, variables are classified based on the set of constraints and on the measured variables in order to determine which ones are determinable (not measured but deductible), free (not measured and not deductible), redundant (measured and staying determinable if the measure is dropped), just measured (measured and becoming free if the measure is dropped).

3. Knowing which variables are “free” and knowing the unique solution for other variables, an interval (solution space) is computed for each free variable.

STEEP already developed a python program to conduct these 3 steps. However it proves very inefficient for large problems (about 20,000 variables). Another issue is that in practice, when one builds a large supply chain model, it often occurs that the reconciliation is infeasible because of conflicting constraints, and the manual identification of these constraints is a tedious job. The purpose of our work was therefore to make improvements on these two fronts.

In the following sections, we start by presenting the problem of interval reduction (step 3 above) before turning to conflicts identification. We finish with ongoing and future work, particularly, in relationship with variable classification.

### 1.1 Our examples

To illustrate our algorithms, we will take two examples: the first one to computing the solution space of “free” variables and the second one for the conflicts detection algorithm. Both will take the form of a simple flow-graph as a reminder of our context.

<figure><img src="image_3.png" /><figcaption>Figure 1: Without conflicts</figcaption></figure>

In this first example, we suppose that we have measured or determined some variables ($v_0, v_6, v_5$). Others remain missing. We will apply our algorithm of Intervals Reduction (IR) with the hypothesis that every variable should be between 0 and 100.

We will use this second example with our algorithm of Conflicts Detection (CD). Indeed, in this example, we have added two constraints on $v_1$ ($C_4$) and $v_6$ ($C_5$) in order to create an unfeasible system. Those constraints can not be represented by flows because they constrain only one variable. So, we can observe them by looking at the intervals of $v_1$ ([0, 30]) and $v_6$ ([50, 100]).

## 2. INTERVALS REDUCTION
<figure><img src="image_4.png" /><figcaption>Figure 2: With conflicts</figcaption></figure>

To compute new intervals, we work on a linear system. We try to determine for each variables a maximum and a minimum to surround the solution space. Thus, let say we have N constraints $c_0...c_{N-1}$ and M variables $v_0...v_{M-1}$ such as :

$$ \forall i \in [0..N-1], \quad c_i : l_i \le \sum_{k=0}^{M-1} \alpha_{i,k} v_k \le u_i \quad (1) $$

with $l_i$ and $u_i$ the bounds of each constraint and $\alpha_{i,k}$ the coefficient of the variable k in the constraint i.

In our context, we said that some variables would be fixed as far as we know their values. Thus, let F be the number of variables known. All variables that are not "free" are known variables (they are uniquely determined). Let L be the number of unknown variables such as $M = F + L$. We can now decompose our first equation (1) :

$$ \forall i \in [0..N-1], \quad c_i : l_i \le \sum_{f=0}^{F-1} \alpha_{i,f} v_f + \sum_{l=0}^{L-1} \alpha_{i,l} v_l \le u_i \quad (2) $$

Our goal here is to compute the maximum and minimum values that each $v_l$ can achieve. To do so, a few algorithms already exist.

## 2.1 Similar work

An interesting algorithm can be found in the literature. Indeed, this problem appears similar to the optimization problem of the Simplex<fnref n="1" />. The Simplex method consists in approaching the optimal solution of an objective function in a problem under duress. Many Simplex algorithms have been developed and are extremely efficient. Thus, in their paper, Puranik and Sahinidis [2] develop the idea that if we run two Simplex programs for each variable (one to compute the maximum and one for the minimum), we are able to get the interval of each unknown variable. Therefore, we have to run 2 * L Simplex to solve our system. In term of complexity, running a such powerful tool two times on variables sub-problem could take some time. But we will detail the efficiency of that algorithm in the part 2.4.

When I joined the STEEP Team, one algorithm was already developed and was particularly long. This algorithm uses the property of each node in the flow-graph. For each constraint, if we want to isolate one unknown variable of

index k to compute his maximum or minimum, we get :

$$ \text{Let } k \in [0..L-1], \quad \forall i \in [0..N-1], $$

$$ c_i \left\{ \begin{array}{l} l_i - \sum_{f=0}^{F-1} \alpha_{i,f} v_f - \sum_{l=0, l \neq k}^{L-1} \alpha_{i,l} v_l \le \alpha_{i,k} v_k \\ \alpha_{i,k} v_k \le u_i - \sum_{f=0}^{F-1} \alpha_{i,f} v_f - \sum_{l=0, l \neq k}^{L-1} \alpha_{i,l} v_l \end{array} \right. \quad (3) $$

For each $v_l$, we can only get an interval. In order to calculate the value of $v_k$, we need to use the maximum or the minimum of each variable $v_l$ which depends on the coefficient $\alpha_{i,k}$. Thus, we generalize this with four cases as follow. We place ourselves in a context of one constraint with two unknown variables to simplify the inequations. With the graph, on a node, we observe in-flow and out-flow. These flows are represented by the sign of the coefficient $\alpha_X$. Furthermore, when a flow $v_l$ (unknown variable) has the same sign than $v_k$ (unknown variable for which we wish to tighten his interval), we use $v_{lmin}$ (resp. $v_{lmax}$) to compute $v_{kmax}$ (resp. $v_{kmin}$). On the other hand, if coefficients have a different sign, we use $v_{lmin}$ (resp. $v_{lmax}$) to compute $v_{kmin}$ (resp. $v_{kmax}$). In this formalization, the coefficients $\alpha_X$ are positive.

* Case 1: $l \le \alpha_k v_k + \alpha_l v_l \le u$

$$ \left\{ \begin{array}{l} v_{kmax} \le \frac{1}{\alpha_k} (u - \alpha_l v_{lmin}) \\ v_{kmin} \ge \frac{1}{\alpha_k} (l - \alpha_l v_{lmax}) \end{array} \right. \quad (4) $$

* Case 2: $l \le \alpha_k v_k - \alpha_l v_l \le u$

$$ \left\{ \begin{array}{l} v_{kmax} \le \frac{1}{\alpha_k} (u + \alpha_l v_{lmax}) \\ v_{kmin} \ge \frac{1}{\alpha_k} (l + \alpha_l v_{lmin}) \end{array} \right. \quad (5) $$

* Case 3: $l \le -\alpha_k v_k - \alpha_l v_l \le u$

$$ \left\{ \begin{array}{l} v_{kmax} \le \frac{1}{-\alpha_k} (l + \alpha_l v_{lmin}) \\ v_{kmin} \ge \frac{1}{-\alpha_k} (u + \alpha_l v_{lmax}) \end{array} \right. \quad (6) $$

* Case 4: $l \le -\alpha_k v_k + \alpha_l v_l \le u$

$$ \left\{ \begin{array}{l} v_{kmax} \le \frac{1}{-\alpha_k} (l - \alpha_l v_{lmin}) \\ v_{kmin} \ge \frac{1}{-\alpha_k} (u - \alpha_l v_{lmax}) \end{array} \right. \quad (7) $$

In this algorithm, we loop over the unknown variables to compute a maximum and a minimum for each constraint with those properties. We repeat this operation while the intervals are being modified.

## 2.2 Our Algorithm

In order to find a new algorithm, we manually tried to find the intervals for our small example, without looking at existing algorithms. Then, we tried to automate our method in an algorithm. Finally, we compared other algorithms with our own method. We found out that our algorithm was really close to the STEEP Team's one. Actually, the differences came from the order of the for-loop. We loop first over the constraints and then over unknown variables to compute the intervals (instead of doing the opposite). Also, a list of constraints to check was added to increase efficiency. When a variable is modified, we put in this list all constraints in which our variable occurs. We pop out those constraints once they are checked. Thanks to that, we do not need to loop over every constraints each time a variable is modified.

<footnote n="1">¹https://en.wikipedia.org/wiki/Simplex_algorithm</footnote>
**Algorithm 1: Intervals Reduction**

1 **intervals reduction** (*M*, *I*, *V*, *u*, *l*) : 
   **Input** : *M* a matrix of coefficient, *I* a list of the intervals, *V* a list of the unknown variables, *u* and *l* the lists of the boundaries of each constraint

**Output:** a list of the reduced intervals

2 $C \leftarrow$ all the constraints

3 **while** *I* is modified **do**

4 **for** each constraint $c_i$ in *C* **do**
5 &nbsp;&nbsp;&nbsp;&nbsp;Drop $c_i$
6 &nbsp;&nbsp;&nbsp;&nbsp;**for** each unknown variable $v_k$ in $c_i$ **do**
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute $v_{kmax}$ and $v_{kmin}$ (With the four cases (4) - (7))
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $v_{kmax} \le$ **upper-bound** $v_k$ **then**
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Modify the upper-bound $v_k$
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $v_{kmin} \ge$ **lower-bound** $v_k$ **then**
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Modify the lower-bound $v_k$
13 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
14 &nbsp;&nbsp;&nbsp;&nbsp;**if** $v_k$'s interval is modified **then**
15 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Add to *C* every $c_j$ in which $v_k$ occurs
16 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
17 &nbsp;&nbsp;&nbsp;&nbsp;**end**
18 &nbsp;&nbsp;&nbsp;&nbsp;**end**
19 &nbsp;&nbsp;&nbsp;&nbsp;**end**
20 **return** *I*

Here, before this algorithm we use a pre-processing to compute the new boundaries (i.e. fixed values are being subtracted from the $l_i$, $u_i$ bounds) for each constraint as follow :

$$ \forall i \in [|0..N - 1|], \left\{ \begin{array}{l} l_i \leftarrow l_i - \sum_{f=0}^{F-1} \alpha_{i,f} v_f \\ u_i \leftarrow u_i - \sum_{f=0}^{F-1} \alpha_{i,f} v_f \end{array} \right. \qquad (8) $$

Furthermore, our algorithm stops because either we reach a max iterations that we fixed or the intervals are unmodified which end the while loop.

With this method, we approach the optimal intervals with more efficiency. Indeed, by computing in each constrain all the unknown variables, we obtain a better restriction on the other variables. That is because the intervals of more variables are reduced in one loop. And the more there are restrictions, the more the intervals are quickly tighten to their optimal values.

Moreover, we can apply a few optimization methods easily (such as parallelism) as the form of the two for-loop are simple (thus easily parallelisable). This tool can drastically increase the speed of our algorithm. To confirm such predicate, we have to create tests that will give us the performance data for our analysis. But first, let us illustrate our algorithm with our example.

## 2.3 Application to our example

With our example (Figure 3), the system will be stored in

a matrix as follow :

$$ \begin{array}{c} v_0 \quad v_1 \quad v_2 \quad v_3 \quad v_4 \quad v_5 \quad v_6 \\ C_0 \quad \begin{pmatrix} 1 & -1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & -1 & 0 & -1 \\ 0 & 0 & 1 & 0 & 1 & -1 & 0 \end{pmatrix} \\ C_1 \\ C_2 \\ C_3 \end{array} $$

$v_1, v_2, v_3$ and $v_4$ are our unknown variables here. Once this matrix is constructed and the initialization of our intervals is done, we start the loop over the constraints. $C_0$ will give us new boundaries on $v_1$ and $v_2$ such as :

$$ \left\{ \begin{array}{l} v_{1max} \le \frac{1}{1}(-98 + 0) \\ v_{1min} \ge \frac{1}{1}(-98 + 100) \end{array} \right. \Rightarrow \begin{array}{l} v_{1max} \le 98 \\ v_{1min} \ge -2 \end{array} \qquad (9) $$

We keep the upper-bound only because we want our interval tightened as possible. Thus, $v_1$ is in $[0, 98]$.

$$ \left\{ \begin{array}{l} v_{2max} \le \frac{1}{1}(-98 + 0) \\ v_{2min} \ge \frac{1}{1}(-98 + 98) \end{array} \right. \Rightarrow \begin{array}{l} v_{2max} \le 98 \\ v_{2min} \ge 0 \end{array} \qquad (10) $$

We continue the loop over the constraints while there are still some modifications. When the algorithm ends, we get the following result :

<figure><img src="image_5.png" /><figcaption>Figure 3: Result Intervals Reduction</figcaption></figure>

## 2.4 Efficiency

Our tests work as follow : for each algorithm we run four examples of different sizes and store the results. First, we want to warm-up our CPU. Thus, ten call to our algorithm are launched without any measurements of the execution time. Then, we run one hundred instances on our example and we compute the average of those launches. At last, we do this operation ten times and as a result, we get the average of these ten averages<fnref n="2" />. As for now, the final execution time is quite representative of the efficiency of this algorithm.

In this table, each average (computed like mentioned above) is stored depending on the size of the example and the algorithm employed. Each number is in milliseconds (except with the huge example).

<table><thead><tr><th></th><th>Small<sup>1</sup></th><th>Medium<sup>2</sup></th><th>Medium<sup>3</sup></th><th>Huge<sup>4</sup></th></tr></thead><tbody><tr><td>Simplex</td><td>113</td><td>/</td><td>/</td><td>/</td></tr><tr><td>STEEP</td><td>0.303</td><td>105.30</td><td>438</td><td>~1h</td></tr><tr><td>STEEP (numba)</td><td>0.04</td><td>9.04</td><td>77</td><td>~11min</td></tr><tr><td>Our Algo</td><td>0.02</td><td>2.00</td><td>6</td><td>~1s</td></tr><tr><td>Our Algo (numba)</td><td>0.03</td><td>0.75</td><td>2</td><td>76ms</td></tr><tr><td>Our Algo (pybind)</td><td>0.004</td><td>0.12</td><td>0.4</td><td>31ms</td></tr></tbody></table>

We developed in python language all those algorithms. However, to use parallelism, we use two kind of tools. The

<footnote n="1">1 around 10 constraints and 4 unknown variables</footnote>

<footnote n="2">2 around 300 constraints and 150 unknown variables</footnote>

<footnote n="3">3 around 800 constraints and 400 unknown variables</footnote>

<footnote n="4">4 around 18 000 constraints and 16 500 unknown variables</footnote>

<footnote n="2">²Theorem Central Limit</footnote>
first one we thought of was the `numba` tool. It is a python library which permits to "translate python function to optimized machine code at runtime"<fnref n="3" />. Thanks to that, we can approach the speeds of a C program. There are many constraints on types of all the arguments, pre-processing was necessary to redefine our parameters. The second tool is called `pybind`. It consists in writing a C++ code of our algorithm with some links to python types. Then compiling it in a python library with many optimization processes. Finally, we can use this library in a python program to compute our intervals faster thanks to parallelism and the C like code.

Concerning the `Simplex` program, we did not succeed to run it on bigger examples and we are trying to understand why. However, this program takes a long time on a very small example compare to the other ones.

## 3. CONFLICTS DETECTION

The conflict detection is used to identify conflicting constraints in order to modify them and make the problem feasible. A variable shouldn't be getting two disjointed intervals with two constraints. In this case, our previous algorithm will notify us. But once our algorithm detect a conflict, we need to know which constraints are involved in this conflict. Indeed, the source of a conflict can be many loop before it is notified. A variable can be modify several times before encounters a constraint in which it will get a disjointed interval. Moreover, the conflict can come from another variable that affects this one. So, how can we determine which constraints should be verified?

### 3.1 Our Algorithm

Our first idea was to store each modification of a variable when IR is running and try to reverse the path that led to a conflict. However, it was extremely difficult to analyze each possibilities to compute the correct one and the complexity of that idea in term of time and space wasn't satisfactory. Then, we discover the work of J.W. Chinneck [1] that Y. Puranik and N.V. Sahinidis [2] took over. In his article, J.W. Chinneck presents two algorithms to find a subset of constraint which is infeasible in a linear program. The first one called the *Deletion Filter*, consists in deleting constraints one by one until the linear program becomes feasible. Once we reach this point, we put back the last constraint dropped in the set and we keep deleting the other constraints. Here is their algorithm (Algorithm 2).

The second algorithm called *Additive Filter* works on the same idea. We add constraints to a set $T$ until this set becomes infeasible. Then we add the last constraint to the $I_S$ infeasible set and start over with $T = I_S$. Here is their algorithm (Algorithm 3).

To detect if a set is infeasible, we use our *Intervals Reduction* algorithm. When a variable gets two disjointed intervals, we have an infeasible set :

$$ 
\begin{align} 
v_{k \max} &\le \text{lower-bound } v_k \text{ or } v_{k \min} \ge \text{upper-bound } v_k \tag{11} \\ 
&\Rightarrow \text{infeasible set} 
\end{align} 
$$

Both algorithms end with one irreducible subset of infeasible constraints. Sometimes, our linear system could have

<caption>**Algorithm 2: Deletion Filter**</caption>

1 **deletion filter (I):**
   Input : *I* a set of infeasible constraints
   Output: an irreducible subset of infeasible constraints
2 **for each constraint cᵢ in I do**
3 &nbsp;&nbsp;&nbsp;&nbsp;Temporarily drop *cᵢ*
4 &nbsp;&nbsp;&nbsp;&nbsp;**if** the set become feasible **then**
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Return *cᵢ* to the set
6 &nbsp;&nbsp;&nbsp;&nbsp;**end**
7 &nbsp;&nbsp;&nbsp;&nbsp;**else**
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Drop the constraint permanently
9 &nbsp;&nbsp;&nbsp;&nbsp;**end**
10 **end**
11 **return** *I*

<caption>**Algorithm 3: Additive Filter**</caption>

1 **additive filter (I):**
   Input : *I* a set of infeasible constraints
   Output: an irreducible subset of infeasible constraints
2 $T = I_S = \emptyset$
3 **while** $I_S$ is feasible **do**
4 &nbsp;&nbsp;&nbsp;&nbsp;$T \leftarrow I_S$
5 &nbsp;&nbsp;&nbsp;&nbsp;**for each constraint cᵢ in I do**
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$T \cup \{cᵢ\}$
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** T is infeasible **then**
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$I_S \leftarrow I_S \cup \{cᵢ\}$
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;GOTO while
10 &nbsp;&nbsp;&nbsp;&nbsp;**end**
11 &nbsp;&nbsp;&nbsp;&nbsp;**end**
12 **end**
13 **return** $I_S$

<footnote n="3">³according to https://numba.pydata.org/</footnote>
more than one conflict. To include this possibility, we add a while loop around both algorithms which loop if the set *I* is still infeasible without the irreducible subset just found. With this modification we are able to collect all the irreducible subsets of our linear system.

Now that we have two algorithms that compute irreducible infeasible subsets of constraints, we want to imagine a new one that could give us the same result faster. Indeed, in both algorithms, constraints are added or deleted one by one and in a huge linear system that takes time. Y.Puranik and N.V. Sahinidis [2] mentioned in their paper that those algorithms can be modify by adding or deleting constraints by packs.

With this idea, we thought about a new algorithm. This algorithm consists in apply the *Additive Filter* by adding constraints by packs. Then we have $I_S$ which is not an irreducible infeasible subset but a smaller infeasible subset. We apply the *Deletion Filter* algorithm to this subset. We finally get one irreducible subset. We repeat this until *I* becomes feasible. To sum up, we add constraints by pack and we drop them one by one. We called it the *Hybrid Filter*.

## 3.2 Our example

We will process our second example 2 with the algorithm *Deletion Filter*.

First, $C_0$ is dropped from *I*. The set remains infeasible so this constraint is dropped permanently. Now, $C_1$ is processed, when this constraint is dropped, the subset becomes feasible. That means this $C_1$ is a part of the problem. So, this constraint will remain in the subset. Next, $C_2$ will be equivalent to $C_1$ and $C_3$ will be permanently removed because the subset remains infeasible without them. This algorithm will now process $C_4$ and $C_5$ and find out that they are also a source to our problem. Eventually, this algorithm will give us the following subset: $[C_1, C_2, C_4, C_5]$.

## 3.3 Efficiency

The efficiency (measured with the same process details in the section 2.4) of those three algorithms is stored in this table. Like said in the section 3.1 with the expression (11), we use our IR algorithm. In order to reduce the execution time, we choose the IR algorithm with the *numba* structure which is easier to call.

<table><thead><tr><th></th><th>Small<sup>1</sup></th><th>Medium<sup>2</sup></th><th>Medium<sup>3</sup></th><th>Huge<sup>4</sup></th></tr></thead><tbody><tr><td>Deletion Filter</td><td>0.89ms</td><td>1.2s</td><td>10s</td><td>+1h</td></tr><tr><td>Additive Filter</td><td>0.93ms</td><td>1.2s</td><td>9.5s</td><td>+1h</td></tr><tr><td>Hybrid Filter</td><td>1.60ms</td><td>2.3s</td><td>18s</td><td>+1h</td></tr></tbody></table>

For those three algorithms, a huge example takes too much time. We stop the execution once it exceeds one hour. We notice that our *Hybrid Filter* algorithm is not faster like expected. But, it was predictable. Indeed, in this algorithm, to add a pack of constraints, we are forced by our structure to add them one by one. Thus, the execution time is far more high than expected. The *Deletion Filter* and the *Additive Filter* algorithms show that on particularly large linear system, an algorithm more efficient is needed to find infeasible subsets.

## 4. ONGOING AND FUTURE WORK

We recently realize that our IR algorithm is not totally accurate. Indeed, in a particularly case, the intervals computed are incorrect.

When there are dependencies between unknown variables and constraints, some intervals can not be computed correctly. Here is an example, let say we have a constraint such as $C_0: a = b - c$ and a second one such as $C_1: x + b = y + c$ where a, b and c are unknown and x and y known. With some inequalities, we have:

$$ \begin{cases} 0 \le a - b + c \le 0 \\ y - x \le b - c \le y - x \end{cases} \tag{12} $$

If some constraints look like those in our linear system, some intervals could be computed with infinite value, which is incorrect. Indeed, in this example, our algorithm will compute a with infinite bounds instead of $y - x$. One idea is to modify our matrix to a special form called *RREF* (Reduced Row Echelon Form) which is a form that we use in our variable's classification algorithm [3]. So, this tool can be easily used for our matrix here. However, this tool also takes time on large matrix. Our current and future work is to transform our current code of *RREF* into a *pybind* structure and then apply it to our matrix. Then, we will verify our assumption with tests and demonstrations. If this idea does not work, we will have to find other methods to have an IR algorithm applicable to any example.

Secondly, like we said in our 3.3 section, the *Hybrid Filter* algorithm add constraints one by one because of our data's type. We could think of a new data structure which will allow us to add constraints by pack and truly exploit our idea. Potentially, we could think of new algorithm with a kind of dichotomy to approach infeasible subsets faster.

## 5. CONCLUSION

Both efficiency tables showed in this paper help us to understand where our theoretical ideas were inaccurate and where we could still gain execution time. This paper gives a new tool to the tightening intervals problem with improvements still possible. This algorithm can now be added to the STEEP project in order to increase their time efficiency on their MFA tool. Indeed, we divide here our execution time by thirty six hundred without any parallelism. This work also gives new ideas to develop a high efficiency on conflicts detection with a large number of constraints. Despite our wrong expectations on our experimental measurements with the *Hybrid Filter*, we know how to transform our data structure to increase the performance on this topic. Moreover, the classification method should now be easily improved by the ongoing work on the *RREF* algorithm.

## 6. REFERENCES

[1] J. W. Chinneck. Finding a useful subset of constraints for analysis in an infeasible linear program. *INFORMS Journal on Computing*, 9(2):164-174, 1997.

[2] Y. Puranik and N. V. Sahinidis. Deletion presolve for accelerating infeasibility diagnosis in optimization models. *INFORMS Journal on Computing*, 29(4):754-766, 2017.

[3] V. Ververka and F. Madron. Material and energy balancing in the process industries, 1997.

<footnote n="1">¹ around 10 constraints and 4 unknown variables</footnote>

<footnote n="2">² around 300 constraints and 150 unknown variables</footnote>

<footnote n="3">³ around 800 constraints and 400 unknown variables</footnote>

<footnote n="4">⁴ around 18 000 constraints and 16 500 unknown variables</footnote>