How to generate multiple solutions to linear program in Python?

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.

• If you're going to down-vote, please offer a reason for this so I can adjust the question accordingly. It is not helpful otherwise Jan 13 at 5:08
• Please, be aware that you posted your question in another community yesterday and duplicated here. Noted that, you should take time to the community folk to answer your question. :) Jan 13 at 9:05
• Cross-site dupe (SO): How to solve for multiple solutions to linear program in Python? Jan 13 at 11:32
• Thanks @A.Omidi, I am not familiar with etiquette here but I will note this for next time. I only duplicated because I realized later that the OR group would probably be a better fit for this type of question. Thanks again Jan 14 at 4:44

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.

• I'd like to add that, instead of solving multiple MIPs (and adding cuts) it is also possible to use the Solution Pool in solvers like Cplex and Gurobi. This is usually way more efficient. Anyway, I think using binary variables to encode the basis is a neat trick. Jan 15 at 15:24
• Thanks for your recommendations @SimonT and your assistance Erwin, I think this method of cutting and excluding past solutions could be useful and relatively straightforward to implement. Ideally, each new solution shouldn't be too similar to the previous solution (I'd like some degree of variability each time it is run), so I am hoping this method can achieve that. I had been solving this previously by using a random sample within each variable's bounds, and back-solving for the third variable (only breaking once the final back-solved variable is in the correct range). Jan 17 at 13:16

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.

• Thanks @ttnick, I will need to look more into how this would work for my given case, but I do find this approach interesting if I can turn it into an algorithm for these types of problems. Jan 17 at 13:41