\[\boxed{\begin{align} &3 \le 5x_1+3x_2+3x_3+5x_4 \le 5 \\ &0 \le x_2+x_3 \le 1 \\ &0 \le x_1+x_4 \le 1 \\ &1 \le 3x_3+5x_4 \le 3 \\ &0 \le 5x_1+3x_2 \le 2 \\ &0 \le x_1,x_2,x_3,x_4 \le 1 \end{align}}\]
The approach I follow is based on this post. That is, we use binary variables to explicitly encode the LP basis, and then solve a series of MIP models where we add cuts to the model to exclude the currently found basis. There are two main differences with the model shown in the original post:
- The current problem has no objective. This is not a major hurdle: we can add a dummy objective. The algorithm becomes actually slightly easier as we don't have to check if the objective starts to deteriorate.
- The current problem has bounded variables. In the original problem we assumed we have positive variables without an upper bound. That meant we could use a single binary variable to encode the basis status: a variable is basic or non-basic. In the current case we could do the same thing and implement the bounds as explicit constraints. However, instead I choose for a slightly more sophisticated approach: we allow the basis status to assume three different values: basic, non-basic at lower bound and non-basic at upper bound. This is exactly how Simplex based linear programming solvers deal with bounded variables.
Positive variables/slacks
Let’s recapitulate what we did in the case of positive variables. When we have positive variables (or slacks) \(x_k\ge 0\), we can use the following basis encoding:
\[\beta_k = \begin{cases} 1 \> \text{if \(x_k\) is basic}\\
0 \> \text{if \(x_k\) is non-basic} \end{cases}\]
If a variable is non-basic it should be equal to zero. This is implemented as
\[0 \le x_k \le M \beta_k\]
In addition we know that we have \(m\) (number of LP rows) basics:
\[\sum_k \beta_k = m\]
Bounded variables/slacks
To handle bounded variables (and slacks): \(\ell_k \le x_k \le u_k\), we introduce an index \(b=\{\text{B,NBLO,NBUP}\}\), where B means basic (variable is between its bounds), NBLO means nonbasic at lower bound (i.e. \(x_k = \ell_k\)) and NBUP means nonbasic upper (\(x_k=u_k\)). Now we can write:
\[\beta_{k,b} = \begin{cases} 1 \> \text{if \(x_k\) has basis status \(b\)}\\
0 \> \text{otherwise}\end{cases}\]
Of course each variable can assume only one (and exactly one) basis status:
\[\sum_b \beta_{k,b} = 1 \> \forall k\]
If a variable is non-basic, it should be at bound:
\[\begin{align} & \ell_k \le x_k \le \ell_k + (u_k-\ell_k)(1-\beta_{k,\text{NBLO}})\\
& u_k \ge x_k \ge u_k - (u_k-\ell_k)(1-\beta_{k,\text{NBUP}}) \end{align}\]
Again, we have \(m\) basics:
\[\sum_k \beta_{k,\text{B}} = m\]
Solution
After running the model we have a large number of feasible bases, 767 to be precise. Some partial results are listed below:
The top part are the values of the columns and rows. The bottom part shows the basis status. We see we get some duplicate values (they have a different basis, but the same values; this is called degeneracy in linear programming).
We can easily add some code to the model to find the unique values and filter out the duplicates:
There are only 20 unique corner point solutions.
We had to run 768 models for that (the last model was infeasible so we stopped), and every model became larger by one constraint. That also means models are becoming more expensive to solve. This is illustrated in this graph: both the number of Simplex Iterations needed to solve each model (orange line) and the solution times (blue line) show an increase over time.
Model
Below the complete model is listed:
| *----------------------------------- |
No comments:
Post a Comment