2
$\begingroup$

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$. Output Solver

$\endgroup$
12
  • $\begingroup$ Is there any other constraint involving the vector $a$? $\endgroup$ Feb 17, 2023 at 14:37
  • $\begingroup$ If you can show the log, that might provide some clues. $\endgroup$
    – RobPratt
    Feb 17, 2023 at 14:45
  • $\begingroup$ Can you add a small working example with all the values for s, p, N, ... ? $\endgroup$ Feb 17, 2023 at 14:53
  • $\begingroup$ 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, ... $\endgroup$ Feb 17, 2023 at 15:02
  • $\begingroup$ 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 $\endgroup$ Feb 17, 2023 at 15:09

1 Answer 1

2
$\begingroup$

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}

$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$ Feb 17, 2023 at 16:22
  • $\begingroup$ 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... $\endgroup$ Feb 17, 2023 at 17:01
  • $\begingroup$ Perhaps it is the time limit, as Rob Pratt suggested. Could that be the case? That it just takes a lot of time? $\endgroup$ Feb 17, 2023 at 17:03
  • $\begingroup$ I tried N = 20 and it has been solved immediately $\endgroup$ Feb 17, 2023 at 17:10
  • 1
    $\begingroup$ 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. $\endgroup$ Feb 17, 2023 at 17:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.