Python -- Optimize system of inequalities

Question

I am working on a program in Python in which a small part involves optimizing a system of equations / inequalities. Ideally, I would have wanted to do as can be done in Modelica, write out the equations and let the solver take care of it.

The operation of solvers and linear programming is a bit out of my comfort zone, but I decided to try anyway. The problem is that the general design of the program is object-oriented, and there are many different possibilities of combinations to form up the equations, as well as some non-linearities, so I have not been able to translate this into a linear programming problem (but I might be wrong) .

After some research I found that the Z3 solver seemed to do what I wanted. I came up with this (this looks like a typical case of what I would like to optimize):

from z3 import *

a = Real('a')
b = Real('b')
c = Real('c')
d = Real('d')
e = Real('e')
g = Real('g')
f = Real('f')
cost = Real('cost')

opt = Optimize()
opt.add(a + b - 350 == 0)
opt.add(a - g == 0)
opt.add(c - 400 == 0)
opt.add(b - d * 0.45 == 0)
opt.add(c - f - e - d == 0)
opt.add(d <= 250)
opt.add(e <= 250)

opt.add(cost == If(f > 0, f * 50, f * 0.4) + e * 40 + d * 20 + 
  If(g > 0, g * 50, g * 0.54))

h = opt.minimize(cost)
opt.check()
opt.lower(h)
opt.model()

Now this works, and gives me the result I want, despite it not being extremely fast (I need to solve such systems several thousands of times). But I am not sure I am using the right tool for the job (Z3 is a "theorem prover").

The API is basically exactly what I need, but I would be curious if other packages allow a similar syntax. Or should I try to formulate the problem in a different way to allow a standard LP approach? (although I have no idea how)

Microsecond

An LP solver should be able to solve this in about 0 seconds, as this is really tiny. – Integrand

Yes, but how can I manage a number of "if" conditions in the function? The solutions I've seen seem kinda hackish, and would not really work in my case. – Microsecond

I think this can be implemented using variable splitting. Ie. introduce non-negative variables fplus,fmin. Add the constraint f=fplus-fmin. The first if then becomes: 50 fplus-0.4 fmin. – Integrand

Answer 1

Z3 is the most powerful solver I have found for such flexible systems of equations. Z3 is an excellent choice now that it is released under the MIT license.

There are a lot of different types of tools with overlapping use cases. You mentioned linear programming -- there are also theorem provers, SMT solvers, and many other types of tools. Despite marketing itself as a theorem prover, Z3 is often marketed as an SMT solver. At the moment, SMT solvers are leading the pack for the flexible and automated solution of coupled algebraic equations and inequalities over the booleans, reals, and integers, and in the world of SMT solvers, Z3 is king. Take a look at the results of the last SMT comp if you want evidence of this. That being said, if your equations are all linear, then you might also find better performance with CVC4. It doesn't hurt to shop around.

If your equations have a very controlled form (for example, minimize some function subject to some constraints) then you might be able to get better performance using a numerical library such as GSL or NAG. However, if you really need the flexibility, then I doubt you are going to find a better tool than Z3.

Putnem

Answer 2

The best solution will probably be to use an ILP solver. Your problem can be formulated as an integer linear programming (ILP) instance. There are many ILP solvers, and some might perform better than Z3. For only 7 variables, any decent ILP solver should find a solution very rapidly.

The only tricky bit are the conditional expressions (If(...)). However, as @Erwin Kalvelagen suggests, the conditionals can be handled using variable splitting. For instance, introduce variables fplus and fminus, with the constraints f = fplus - fminus and fplus >= 0 and fminus >= 0. Now you can replace If(f > 0, f * 50, f * 0.4) with 50 * fplus - 0.4 * fminus. In this case, that will be equivalent.

Variable splitting doesn't always work. You have to think about whether it might introduce spurious solutions (where both fplus > 0 and fminus > 0). In this case, though, spurious solutions will never be optimal -- one can show that the optimal solution will never be optimal. Consequently, variable splitting works fine here.

If you have a situation where you do have conditional statements but variable splitting doesn't work, you can often use the techniques at https://cs.stackexchange.com/q/12102/755 to formulate the problem as an instance of ILP.

Glaikit