solving a system of inequalities where the solution is infeasible

I am trying to solve the attached a non-linear optimization problem. I would like to give mystic a try, since SciPy.optimize does not work as expected (further details below).

Because c_1 = 3, the problem is infeasible. This is because (p_1 - 3 - 0.22) / p_1 < 0.05 means that p_1 would have to be larger than 3.22 which would conflict with p_1 / 2.2 <= 0.65.

There exists an unanswered question addressing the fundamental problem of SciPy successfully terminating for infeasible problems here. Unfortunately, I need this specific problem solved with python, which is why I am trying my luck here again.

When I used SciPy.optimize, the program also terminated successfully, breaking the constraints without raising an issue (even when setting keep_feasible=True). My problem is larger and other constraints could not be met, either. Accordingly, SciPy appears to be the wrong tool for the job.

My first question: am I doing anything wrong? If not: Are there any alternatives to SciPy.optimize? I have also been looking into mystic, but not been able to get it to work at all.

I know this is a rather specific problem, so my thanks to everyone who wants to chip in. Finally, my apologies if I have offended any mathematicians with the probably incorrect formulation of the problem.

This is a bit of a difficult problems, since it requires solving a system of symbolic inequality equations with denominators. A symbolic solve with inequalities with denominators is tricky because multiplying by the denominator can flip the sign of the inequality, depending on the value of the yet-unknown variable. Solving one of these alone is difficult... but then adding in that the solution is infeasible makes it even worse. So, bear in mind that the following is a bit shaky... but it's probably the close to the best you can do to approach this problem directly.

>>> def objective(x):
...   x0,x1,x2 = x
...   return (x0 - 3 - .23)/x0 + (x1 - 1.17 - .23)/x1 + (x2 - 0.71 - .23)/x2
...
>>> def cost(x): # maximize
...   return -objective(x)
...
>>> bounds = [(None,.65*2.2),(None,.65*2.9),(None,.65*1.91)]
>>>
>>> equations = '''
... (x0 - 3.00 - .23)/(38000 * x0) >= .05
... (x1 - 1.17 - .23)/(33000 * x1) >= .05
... (x2 - 0.71 - .23)/(29000 * x2) >= .05
... (38000 * x0 + 33000 * x1 + 29000 * x2) / (38000 * (2.2 - x0) + 33000 * (2.9 - x1) + 29000 * (1.91 - x2)) <= -0.5
... (38000 * (2.2 - x0) + 33000 * (2.9 - x1) + 29000 * (1.91 - x2)) / (38000 * x0 + 33000 * x1 + 29000 * x2) >= 0.18
... '''
>>> from mystic.solvers import diffev2
>>> from mystic.symbolic import generate_constraint, generate_solvers, simplify

That should set up the problem... now to the tricky part. Simplify the equations, and then solve.

>>> eqn = simplify(equations)
>>> cf = generate_constraint(generate_solvers(eqn))
>>> result = diffev2(cost, x0=bounds, bounds=bounds, constraints=cf, npop=40, gtol=50, disp=True, full_output=True)
Warning: Maximum number of iterations has been exceeded
>>> print(result[1])
inf

We see that we get inf, which is an indication that mystic is running into an infeasible solution. But, I'm glossing over some stuff above, and was a bit lucky, actually. Let's look at the simplified equation.

>>> print(eqn)
x0 > 0
x0 > -0.868421052631579*x1 - 0.763157894736842*x2 + 6.17605263157895
x2 >= -0.000648723257418910
x1 >= -0.000848999393571862
x0 >= -0.868421052631579*x1 - 0.763157894736842*x2 - 6.17605263157895
x2 < 0
x1 < 0
x0 < -33*x1/38 - 29*x2/38
x0 <= -0.00170089520800421
x0 >= -0.868421052631579*x1 - 0.763157894736842*x2 + 5.23394290811775

That's only one possible of the several simplified equations, due to the denominators and inequalities, as mentioned above... it just happened to be the right one. We can get all of the candidate simplified equations, as follows:

>>> all_eqn = simplify(equations, all=True)
>>> len(all_eqn)
32

Note we have 32 possibilities -- some of which are "good", and some aren't -- depending on the values of the variables. You can plug all_eqn into cf = generate_constraint(generate_solvers(all_eqn)), and then mystic will do it's best to find the solution that minimizes all the possible simplified equations. This typically works ok... and I'm sure may fail. The situations gets worse with there being a maximization, and more importantly, the solution you are looking for is infeasible.

I'm going to say that this is an area of active development, and mystic could use some improvement to handle cases like this better.

EDIT: In the solution above, I forgot use the keyword join. The default is join=None which applies each equation sequentially. While this is faster, it's not what you want if there's conflicting equations. I should probably have used:

>>> from mystic.constraints import and_
>>> cf = generate_constraint(generate_solvers(eqn), join=and_)

which should better ensure that all constraints are respected.