How to generate multiple solutions to linear program in Python?

Question

I recently posted a question here where I had written some Pyomo code to solve a simple linear program, however I feel the Operations Research platform is perhaps a better fit for this question. My problem is that every time I run this, I get just one optimal solution. I am wondering whether there is a way to produce a different solution every time I run the code (it doesn't have to always be a unique solution)?

The problem:

$$ 3a+2b+2c=15$$

where $$-20\leqslant a \leqslant 0$$

$$ 0 \leqslant b \leqslant 20$$

$$ 20 \leqslant c \leqslant 50$$

I am also open to solutions in Pulp as well of Pyomo.

Answer 1

When I understand your question correctly you are not interested in the full solution space. In that case you can generate alternative optimal solutions by first solving your problem to optimality. Then you can add the objective function with the optimal value as a constraint and resolve with a random objective vector. For example for the problem.

$ \begin{align} \min && c^Tx\\ \text{s.t.}&& Ax\leq b\\ && x\geq 0 \end{align} $

Let $C^*$ be the optimal objective value you get from solving the initial problem and let $c'$ be an other random objective vector. Then solving

$ \begin{align} \min && c'^Tx\\ \text{s.t.}&& Ax&\leq b\\ && c^Tx &= C^*\\ && x&\geq 0 \end{align} $

will give you a solution that is optimal with respect to the original objective $\min c^Tx$.

If you are interested in getting all optimal solutions you will need to find all optimal basis solutions as pointed you by @Erwin Kalvelagen here.

There are multiple ways to do so. Either you solve multiple MIPs where in each iteration you add cuts that prevent you from finding previously found basis solutions again. It is described in this blog post. As pointed out by Erwin Kalvelagen some solvers like gurobi provide features that allow to extract multiple optimal solutions. If these are not available multiple MIPs might to be solved, which could result in quite big runtimes.

The more direct way would be to modify the Simplex Algorithm to not exchange for better basis solutions, but for equal ones. Starting with an optimal solution (which you can get from any LP solver you want) you choose entering variables with entry in the objective row $=0$. For leaving variables you proceed as usual. To ensure you find all solutions you can use classic graph traversal techniques.

Note that any optimal solution can then be expressed as linear combinations of the optimal basis solutions.

Answer 2

Your problem does not have any objective function, so you do not have to use software to solve linear programs here. It is just a system of linear inequalities. (Of course, you still can use LP tools by just minimizing a constant objective function.)

May I suggest a solution that neither requires Pyomo nor Pulp? You can solve those systems by using the Fourier-Motzkin Method. This would give you either no solution (if solution space is empty), one solution (if unique), or infinitely many (if solid). In the latter case, you should be able to get parametrized bounds for each variable depending on which order you eliminate the variables.