# Implementing a MINLP problem in Pyomo: giving an infeasible solution

I'm trying to implement a MINLP problem which is described in a previous post here: How do we formulate a problem where the decision variable has an index that is also a decision variable? The only change is that the objective function is different and is to be minimized instead of to be maximized (the function has changed but is given below for completeness).

The problem is about minimizing a (nonlinear) objective function w.r.t. $$x_i$$:

$$\min \sum_{i=1}^N f(x_i, s_i, p_i),$$

where $$s_i$$ and $$p_i$$ are constants for each $$i=1, 2, \dots, N$$ and $$x_i$$ is a decision variable, where $$0 \le x_i < 1$$. The function $$f$$ is

$$f(x_i, s_i, p_i) = p_is_i\frac{x_i^{0.135}-(1-x_i)^{0.135}}{0.1975}.$$

We also have a binary variable $$y_{ij}$$ and a continuous variable $$a_j$$ where $$0 \le a_j < 1$$ and $$j = 1, 2, 3$$.

The constraints are the following:

\begin{align} \beta= \frac{\sum_{i=1}^N x_id_i}{\sum_{i=1}^N d_i},\tag1\label1 \\ \end{align} where $$d_i$$ is again a constant for each $$i=1, 2, \dots, N$$ and $$\beta$$ is a constant where $$0 \le \beta < 1$$.

\begin{align} \sum_j y_{ij} &= 1 &&\text{for all i}, \tag2\label2 \\ -(1 - y_{ij}) \le x_i - a_j &\le 1 - y_{ij} &&\text{for all i and j}. \tag3\label3 \end{align}

I'm implementing it in Pyomo as interface and I use mindtpy as a solver. Here is the code I have implemented (note that I split up constraint $$3$$).

f = lambda x, s, p: ((x**0.135-(1-x)**0.135)/0.1975)*p*s
Beta = 0.96
model = ConcreteModel(name="ORExchangeExample")

model.x = Var(N, bounds=(0.00, 0.9999), initialize = 0.50)
model.a = Var(J, bounds=(0.00, 0.9999), initialize = 0.50)
model.y = Var(N, J, domain=Binary, initialize=0)

def obj_rule(model):
return sum(f(model.x[n], s[n], p[n]) for n in N)
model.obj = Objective(rule=obj_rule, sense=minimize)

def constraint_1(model):
return sum(d[n]*model.x[n] for n in N) == Beta*(sum(d[n] for n in N))
model.constraint_1= Constraint(rule=constraint_1)

def constraint_2(model, i):
return sum(model.y[i, j] for j in J) == 1
model.constraint_2= Constraint(N, rule=constraint_2)

def constraint_3(model, i, j):
return -(1-model.y[i, j]) <= model.x[i]-model.a[j]
model.constraint_3= Constraint(N, J, rule=constraint_3)

def constraint_4(model, i, j):
return model.x[i]-model.a[j] <= 1-model.y[i, j]
model.constraint_4= Constraint(N, J, rule=constraint_4)

solver.solve(model, tee=True,
mip_solver='glpk',
nlp_solver='ipopt',)
solver.solve(model)


The variable $$N$$ is a list containing $$1, 2 , \dots, N$$ and $$J$$ is also a list $$= [1, 2, 3]$$. The constant variables $$p, s$$ and $$d$$ are dictionaries with the elements of $$N$$ as keys and the actual values as values, for example:

p = {1: 19.06, 2: 4.0607, 3: 9.75, 4: 0.2716, ...}


When I run this, it runs without errors but gives me back as a solution the initialized values for $$x$$ and $$a$$. So it doesn't solve it, because constraint (1) isn't met and the objective function is definitely not minimized. The result (in this code set up) will give [0.5, 0.5, ...] for all values of $$x_i$$.

Am I missing something?

EDIT: I used the Ipopt as the nonlinear solver but still the same result. However, if I do it for a small sample ($$N = 10$$), then it all works fine...

This is what the solver has to say, but I don't really know how to interpret this. In my data, $$N = 4831$$.

• Is there any other constraint involving the vector $a$? Commented Feb 17, 2023 at 14:37
• If you can show the log, that might provide some clues. Commented Feb 17, 2023 at 14:45
• Can you add a small working example with all the values for s, p, N, ... ? Commented Feb 17, 2023 at 14:53
• There is no constraint involving the vector $a$. I will edit the question with the log and add a small working example for s, p, N, ... Commented Feb 17, 2023 at 15:02
• ok, because without any other constraint involving a or if $a$ isn't part of the objective, it acts as bounds for x. So x is constrained to be equal to a=0.5 Commented Feb 17, 2023 at 15:09

Try ipopt as the nonlinear solver. With ipopt I got the following solution from your code (with the following assumptions):

12 Set Declarations
a_index : Size=1, Index=None, Ordered=Insertion
Key  : Dimen : Domain : Size : Members
None :     1 :    Any :    3 : {1, 2, 3}
constraint_2_index : Size=1, Index=None, Ordered=Insertion
Key  : Dimen : Domain : Size : Members
None :     1 :    Any :    2 : {1, 2}
constraint_3_index : Size=1, Index=None, Ordered=True
Key  : Dimen : Domain                                    : Size :
Members
None :     2 : constraint_3_index_0*constraint_3_index_1 :    6 : {(1,
1), (1, 2), (1, 3), (2, 1), (2, 2), (2, 3)}
constraint_3_index_0 : Size=1, Index=None, Ordered=Insertion
Key  : Dimen : Domain : Size : Members
None :     1 :    Any :    2 : {1, 2}
constraint_3_index_1 : Size=1, Index=None, Ordered=Insertion
Key  : Dimen : Domain : Size : Members
None :     1 :    Any :    3 : {1, 2, 3}
constraint_4_index : Size=1, Index=None, Ordered=True
Key  : Dimen : Domain                                    : Size :
Members
None :     2 : constraint_4_index_0*constraint_4_index_1 :    6 : {(1,
1), (1, 2), (1, 3), (2, 1), (2, 2), (2, 3)}
constraint_4_index_0 : Size=1, Index=None, Ordered=Insertion
Key  : Dimen : Domain : Size : Members
None :     1 :    Any :    2 : {1, 2}
constraint_4_index_1 : Size=1, Index=None, Ordered=Insertion
Key  : Dimen : Domain : Size : Members
None :     1 :    Any :    3 : {1, 2, 3}
x_index : Size=1, Index=None, Ordered=Insertion
Key  : Dimen : Domain : Size : Members
None :     1 :    Any :    2 : {1, 2}
y_index : Size=1, Index=None, Ordered=True
Key  : Dimen : Domain              : Size : Members
None :     2 : y_index_0*y_index_1 :    6 : {(1, 1), (1, 2), (1, 3),
(2, 1), (2, 2), (2, 3)}
y_index_0 : Size=1, Index=None, Ordered=Insertion
Key  : Dimen : Domain : Size : Members
None :     1 :    Any :    2 : {1, 2}
y_index_1 : Size=1, Index=None, Ordered=Insertion
Key  : Dimen : Domain : Size : Members
None :     1 :    Any :    3 : {1, 2, 3}

3 Var Declarations
a : Size=3, Index=a_index
Key : Lower : Value              : Upper  : Fixed : Stale : Domain
1 :   0.0 : 0.7327896229031757 : 0.9999 : False : False :  Reals
2 :   0.0 : 0.7327896229031757 : 0.9999 : False : False :  Reals
3 :   0.0 : 0.7327896229031757 : 0.9999 : False : False :  Reals
x : Size=2, Index=x_index
Key : Lower : Value              : Upper  : Fixed : Stale : Domain
1 :   0.0 : 0.9788067730557521 : 0.9999 : False : False :  Reals
2 :   0.0 :  0.950596613472124 : 0.9999 : False : False :  Reals
y : Size=6, Index=y_index
Key    : Lower : Value               : Upper : Fixed : Stale : Domain
(1, 1) :     0 : 0.33333333333333337 :     1 : False : False : Binary
(1, 2) :     0 : 0.33333333333333337 :     1 : False : False : Binary
(1, 3) :     0 : 0.33333333333333337 :     1 : False : False : Binary
(2, 1) :     0 :  0.3333333333333333 :     1 : False : False : Binary
(2, 2) :     0 :  0.3333333333333333 :     1 : False : False : Binary
(2, 3) :     0 :  0.3333333333333333 :     1 : False : False : Binary

1 Objective Declarations
obj : Size=1, Index=None, Active=True
Key  : Active : Sense    : Expression
None :   True : minimize : (x[1]**0.135 - (1 - x[1])**0.135)/0.1975 +
(x[2]**0.135 - (1 - x[2])**0.135)/0.1975*2*2

4 Constraint Declarations
constraint_1 : Size=1, Index=None, Active=True
Key  : Lower : Body          : Upper : Active
None :  2.88 : x[1] + 2*x[2] :  2.88 :   True
constraint_2 : Size=2, Index=constraint_2_index, Active=True
Key : Lower : Body                     : Upper : Active
1 :   1.0 : y[1,1] + y[1,2] + y[1,3] :   1.0 :   True
2 :   1.0 : y[2,1] + y[2,2] + y[2,3] :   1.0 :   True
constraint_3 : Size=6, Index=constraint_3_index, Active=True
Key    : Lower : Body                           : Upper : Active
(1, 1) :  -Inf : - (1 - y[1,1]) - (x[1] - a[1]) :   0.0 :   True
(1, 2) :  -Inf : - (1 - y[1,2]) - (x[1] - a[2]) :   0.0 :   True
(1, 3) :  -Inf : - (1 - y[1,3]) - (x[1] - a[3]) :   0.0 :   True
(2, 1) :  -Inf : - (1 - y[2,1]) - (x[2] - a[1]) :   0.0 :   True
(2, 2) :  -Inf : - (1 - y[2,2]) - (x[2] - a[2]) :   0.0 :   True
(2, 3) :  -Inf : - (1 - y[2,3]) - (x[2] - a[3]) :   0.0 :   True
constraint_4 : Size=6, Index=constraint_4_index, Active=True
Key    : Lower : Body                       : Upper : Active
(1, 1) :  -Inf : x[1] - a[1] - (1 - y[1,1]) :   0.0 :   True
(1, 2) :  -Inf : x[1] - a[2] - (1 - y[1,2]) :   0.0 :   True
(1, 3) :  -Inf : x[1] - a[3] - (1 - y[1,3]) :   0.0 :   True
(2, 1) :  -Inf : x[2] - a[1] - (1 - y[2,1]) :   0.0 :   True
(2, 2) :  -Inf : x[2] - a[2] - (1 - y[2,2]) :   0.0 :   True
(2, 3) :  -Inf : x[2] - a[3] - (1 - y[2,3]) :   0.0 :   True

20 Declarations: x_index x a_index a y_index_0 y_index_1 y_index y obj
constraint_1 constraint_2_index constraint_2 constraint_3_index_0
constraint_3_index_1 constraint_3_index constraint_3 constraint_4_index_0
constraint_4_index_1 constraint_4_index constraint_4
Objective value:  8.660160748261013
x[1] value is 0.978807
x[2] value is 0.950597

Process finished with exit code 0


Assumptions: n_ = 3 N = list(range(1, n_)) J = [1, 2, 3] s = {1:1, 2:2, 3:3} p = {1:1, 2:2, 3:3} d = {1:1, 2:2, 3:3}

• Thanks! I used the Ipopt as the nonlinear solver, but still the same results. Do you have any other pointers I can take a look at? If I restrict $N=10$, it works fine. So it seems to have something to do with the magnitude of the problem. Commented Feb 17, 2023 at 16:22
• I don't think changing N from 10 to 11 makes the problem infeasible unless there is some other hidden errors in the indexing of the variables... Commented Feb 17, 2023 at 17:01
• Perhaps it is the time limit, as Rob Pratt suggested. Could that be the case? That it just takes a lot of time? Commented Feb 17, 2023 at 17:03
• I tried N = 20 and it has been solved immediately Commented Feb 17, 2023 at 17:10
• Yes, it seems to be hard when $N$ grows larger and it depends a lot on the values in it I see. I'll accept the answer, because it the implementation seems correct as you showed. Commented Feb 17, 2023 at 17:44