Fitting 16 variables in 17 inequalities in Python

Question

I have 17 inequalities:

1:  (x[0] + x[1] + x[2] + x[4] + x[6] + x[7])/(x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[7]) >= 0.4
2:  (x[3] + x[5])/(x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[7]) <= 0.6
3:  x[0]/(x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[7]) <= 0.1
4:  (x[0] + x[2] + x[4] + x[6])/(x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[7]) <= 0.1
5:  (1.2*x[11]*x[3] + 1.2*x[13]*x[5])/(x[3] + x[5]) <= 520
6:  (1.2*x[11]*x[3] + 1.2*x[13]*x[5])/(x[3] + x[5]) >= 470
7:  (x[0]*x[8] + x[1]*x[9] + x[10]*x[2] + x[12]*x[4] + x[14]*x[6] + x[15]*x[7])/(x[0] + x[1] + x[2] + x[4] + x[6] + x[7]) <= 420
8:  (x[3] + x[7])/x[5] >= 0.05
9:  (x[3] + x[7])/x[5] <= 0.2
10:  x[2]/(x[3] + x[7]) >= 0.05
11:  x[2]/(x[3] + x[7]) <= 0.15
12:  x[5]/(x[3] + x[5]) >= 0.95
13:  0.833333333333333/x[11] >= 376
14:  0.833333333333333/x[11] <= 424 
15:  x[13]/x[11] >= 0.7
16:  x[13]/x[11] <= 0.82
17:  1.2*x[11]*x[3] + 1.2*x[13]*x[5] <= 317300.0

And 16 variables.

I need to fit values of variables to these constraints.

Also: x[0-7] values must be between 20 and 600 x[8-15] values must be between 200 and 600

I tried scipy.optimize.minimize() optimizer with SLSQP method with using simple function as sum of all variables. But I guess it is not guarantied that minimize() satisfy all constrains. Some of my constraints are ignored.

I don't even need to minimize all values but I can't find any other solution.

P.S. I am far from math and programming, so sorry for simple mistakes.

Answer 1

Erwin is correct: no solver will be able to solve this completely because it's infeasible. Picking and choosing your constraints, relaxing one side each of the double-sided constraints that are non-satisfiable (11, 13), a partially valid solution is

import numpy as np
import scipy.optimize
from scipy.optimize import Bounds, LinearConstraint, NonlinearConstraint


def nop(_) -> int:
    return 0


def nonlin_constraint(x: np.ndarray) -> tuple[float, ...]:
    c56 = (1.2*x[11]*x[3] + 1.2*x[13]*x[5])/(x[3] + x[5])
    c7 = (
         x[0]*x[8] + x[1]*x[9] + x[10]*x[2] + x[12]*x[4] + x[14]*x[6] + x[15]*x[7]
     ) / (x[0] + x[1] + x[2] + x[4] + x[6] + x[7])
    c17 = 1.2*x[11]*x[3] + 1.2*x[13]*x[5]
    return c56, c7, c17


result = scipy.optimize.minimize(
    fun=nop,
    x0=np.full(shape=16, fill_value=200),
    bounds=Bounds(
        lb=[20]*8 + [200]*8,
        ub=[600]*16,
    ),
    constraints=[
        LinearConstraint(
            A=np.array((
                #    0      1      2      3      4      5      6      7   8-15
                (1-0.4, 1-0.4, 1-0.4,  -0.4, 1-0.4,  -0.4, 1-0.4, 1-0.4,  0,0,0,0,0,0,0,0),  # 1
                ( -0.6,  -0.6,  -0.6, 1-0.6,  -0.6, 1-0.6,  -0.6,  -0.6,  0,0,0,0,0,0,0,0),  # 2
                (1-0.1,  -0.1,  -0.1,  -0.1,  -0.1,  -0.1,  -0.1,  -0.1,  0,0,0,0,0,0,0,0),  # 3
                (1-0.1,  -0.1, 1-0.1,  -0.1, 1-0.1,  -0.1, 1-0.1,  -0.1,  0,0,0,0,0,0,0,0),  # 4
                (    0,     0,     0,     1,     0, -0.05,     0,     1,  0,0,0,0,0,0,0,0),  # 8
                (    0,     0,     0,     1,     0,  -0.2,     0,     1,  0,0,0,0,0,0,0,0),  # 9
                (    0,     0,     1, -0.05,     0,     0,     0, -0.05,  0,0,0,0,0,0,0,0),  # 10
                (    0,     0,     1, -0.15,     0,     0,     0, -0.15,  0,0,0,0,0,0,0,0),  # 11
                (    0,     0,     0, -0.95,     0,1-0.95,     0,     0,  0,0,0,0,0,0,0,0),  # 12
                (    0,     0,     0,     0,     0,     0,     0,     0,  0,0,0, -376/(5/6), 0,0,0,0),  # 13
                (    0,     0,     0,     0,     0,     0,     0,     0,  0,0,0, -424/(5/6), 0,0,0,0),  # 14
                (0,0,0,0, 0,0,0,0,                                        0,0,0, -0.7, 0, 1, 0,0),  # 15
                (0,0,0,0, 0,0,0,0,                                        0,0,0,-0.82, 0, 1, 0,0),  # 16
            )),
            lb=np.array((
                0,        # 1
                -np.inf,  # 2
                -np.inf,  # 3
                -np.inf,  # 4
                0,        # 8
                -np.inf,  # 9
                0,        # 10
                -np.inf,  # 11
                0,        # 12
                -1e6,     # 13: unsatisfiable; this should be '0'
                -np.inf,  # 14
                0,        # 15
                -np.inf,  # 16
            )),
            ub=np.array((
                np.inf,  # 1
                0,       # 2
                0,       # 3
                0,       # 4
                np.inf,  # 8
                0,       # 9
                np.inf,  # 10
                10,      # 11: unsatisfiable; this should be '0'
                np.inf,  # 12
                np.inf,  # 13
                0,       # 14
                np.inf,  # 15
                0,       # 16
            ))
        ),
        NonlinearConstraint(
            fun=nonlin_constraint,
            lb=(
                470,      # 5, 6
                -np.inf,  # 7
                -np.inf,  # 17
            ),
            ub=(
                520,      # 5, 6
                420,      # 7
                317300,   # 17
            ),
        ),
    ],
)

x = result.x
c56, c7, c17 = nonlin_constraint(x)
print(result.message)
print(x)

# All inequalities evaluated for validation
print('1:', (x[0] + x[1] + x[2] + x[4] + x[6] + x[7])/(x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[7]), '>= 0.4')
print('2:', (x[3] + x[5])/(x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[7]), '<= 0.6')
print('3:', x[0]/(x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[7]), '<= 0.1')
print('4:', (x[0] + x[2] + x[4] + x[6])/(x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[7]), '<= 0.1')
print('5:', c56, '<= 520')
print('6:', c56, '>= 470')
print('7:', c7, '<= 420')
print('8:', (x[3] + x[7])/x[5], '>= 0.05')
print('9:', (x[3] + x[7])/x[5], '<= 0.2')
print('10:', x[2]/(x[3] + x[7]), '>= 0.05')
print('11:', x[2]/(x[3] + x[7]), '<= 0.15 (unsatisfiable)'),
print('12:', x[5]/(x[3] + x[5]), '>= 0.95')
print('13:', 0.833333333333333/x[11], '>= 376 (unsatisfiable)')
print('14:', 0.833333333333333/x[11], '<= 424')
print('15:', x[13]/x[11], '>= 0.7')
print('16:', x[13]/x[11], '<= 0.82')
print('17:', c17, '<= 317300')

Optimization terminated successfully
[ 20.         267.44163498  20.          20.48131036  20.
 389.14489686  20.          48.2563213  200.         200.
 200.         472.45677523 200.         387.41455569 200.
 200.        ]
1: 0.49135239474239856 >= 0.4
2: 0.5086476052576016 <= 0.6
3: 0.024834719863914737 <= 0.1
4: 0.09933887945565895 <= 0.1
5: 470.000000000485 <= 520
6: 470.000000000485 >= 470
7: 200.0 <= 420
8: 0.17663762833141627 >= 0.05
9: 0.17663762833141627 <= 0.2
10: 0.2909614357900084 >= 0.05
11: 0.2909614357900084 <= 0.15 (unsatisfiable)
12: 0.95 >= 0.95
13: 0.0017638297872322214 >= 376 (unsatisfiable)
14: 0.0017638297872322214 <= 424
15: 0.8199999999999998 >= 0.7
16: 0.8199999999999998 <= 0.82
17: 192524.31739218766 <= 317300

Removing the relaxation constraint values from constraints 11 or 13 will instead produce this approximated result violating constraints 3, 4, 6, 9, 11, 12, 13, and 16:

Positive directional derivative for linesearch
[200. 200. 200. 200. 200. 200. 200. 200. 200. 200. 200. 200. 200. 200.
 200. 200.]
1: 0.75 >= 0.4
2: 0.25 <= 0.6
3: 0.125 <= 0.1
4: 0.5 <= 0.1
5: 240.0 <= 520
6: 240.0 >= 470
7: 200.0 <= 420
8: 2.0 >= 0.05
9: 2.0 <= 0.2
10: 0.5 >= 0.05
11: 0.5 <= 0.15
12: 0.5 >= 0.95
13: 0.004166666666666665 >= 376
14: 0.004166666666666665 <= 424
15: 1.0 >= 0.7
16: 1.0 <= 0.82
17: 96000.0 <= 317300

Even constraint 13 on its own doesn't make any sense. Given the bounds of variable 11, expression 13 can only vary between

5/6/200 = 0.004166666666666667

and

5/6/600 = 0.001388888888888889

so there's absolutely no way it's going above 376.