Modelling a decision variable as an index of a (fixed) set

Question

I'm trying to model the following MINLP problem in Pyomo.

We are trying to minimize a nonlinear objective function $f$ in $x_i \in \lbrace{0, 1, 2\rbrace}$ for $i= 1, 2, \dots, N$ and subject to a constraint involving $x_i$. Then, $a_j \in \lbrace{a_0, a_1, a_2\rbrace}$ is a fixed set where $0 \le a_0, a_1, a_2 < 1$ and the problem can be described as

$$\min f(a_{x_i}, s_i, p_i),$$

s.t.

$$\beta = \frac{\sum_{i=1}^N d_ia_{x_i}}{\sum_{i=1}^N d_i},$$ where $d_i, p_i$ and $s_i$ are constants for each $i=1, 2, \dots, N$ and $\beta$ is a fixed number where $0 \le \beta < 1$.

The function $f(a_{x_i}, s_i, p_i)$ is

$$p_is_i\frac{a_{x_i}^{0.135}-(1-a_{x_i})^{0.135}}{0.1975}.$$

This is an example of the code I had at first:

f = lambda x, s, p: ((x**0.135-(1-x)**0.135)/0.1975)*p*s
Beta = 0.96
N = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
a = [0.98, 0.95, 0.90]
p = {1: 2.1, 2: 3.4, 3: 1.7, 4: 2.4, 5: 5.3, 6: 0.8, 7: 8.8, 8: 4.9, 9: 10.6, 10: 3.2}
d = {1: 1000, 2: 870, 3: 500, 4: 480, 5: 370, 6: 325, 7: 280, 8: 190, 9: 145, 10: 115}
s = {1: 25, 2: 36, 3: 17, 4: 14, 5: 46, 6: 14, 7: 56, 8: 11, 9: 8, 10: 13}
model = ConcreteModel(name="SLDIFF")
model.x = Var(N, bounds=(0, 2), domain=NonNegativeIntegers, initialize=0)
def obj_rule(model):
    return sum(f(a[model.x[n]], s[n], p[n]) for n in N)
model.obj = Objective(rule=obj_rule, sense=minimize)
def constraint(model):
    return sum(d[n]*a[model.x[n]] for n in N) == Beta*(sum(d[n] for n in N))
model.constraint = Constraint(rule=constraint)
solver = SolverFactory('mindtpy')
solver.solve(model, mip_solver = 'glpk', nlp_solver='ipopt')

However, Pyomo doesn't allow this. It gives the TypeError: list indices must be integers or slices, not _GeneralVarData.

My first thought was to adjust 'model.x[n]' to 'value(model.x[n])':

f = lambda x, s, p: ((x**0.135-(1-x)**0.135)/0.1975)*p*s
Beta = 0.96
N = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
a = [0.98, 0.95, 0.90]
p = {1: 2.1, 2: 3.4, 3: 1.7, 4: 2.4, 5: 5.3, 6: 0.8, 7: 8.8, 8: 4.9, 9: 10.6, 10: 3.2}
d = {1: 1000, 2: 870, 3: 500, 4: 480, 5: 370, 6: 325, 7: 280, 8: 190, 9: 145, 10: 115}
s = {1: 25, 2: 36, 3: 17, 4: 14, 5: 46, 6: 14, 7: 56, 8: 11, 9: 8, 10: 13}
model = ConcreteModel(name="SLDIFF")
model.x = Var(N, bounds=(0, 2), domain=NonNegativeIntegers, initialize=0)
def obj_rule(model):
    return sum(f(a[value(model.x[n])], s[n], p[n]) for n in N)
model.obj = Objective(rule=obj_rule, sense=minimize)
def constraint(model):
    return sum(d[n]*a[value(model.x[n])] for n in N) == Beta*(sum(d[n] for n in N))
model.constraint = Constraint(rule=constraint)
solver = SolverFactory('mindtpy')
solver.solve(model, mip_solver = 'glpk', nlp_solver='ipopt')
print([value(model.x[key]) for key in model.x])
print(value(model.obj))

But then it gives an error when trying to build the constraint:

ValueError: Invalid constraint expression. The constraint expression resolved to a trivial 
Boolean (False) instead of a Pyomo object. Please modify your rule to return 
Constraint.Infeasible instead of False.

Error thrown for Constraint 'constraint'

But I don't see how this expression is invalid? Anyone experience with this error and how to solve it?

Answer 1

You can't use a decision variable x[n] as an index in a[x[n]].

Say you want

y = a[x] 

Instead of that, introduce binary variables $\delta_i$ and do: $$\begin{align} &y = \sum_i a_i \cdot\delta_i \\ &x = \sum_i i\cdot \delta_i \\ &\sum_i \delta_i = 1 \\ &\delta_i \in \{0,1\} \end{align} $$ This is now completely linear.

Or use a solver or tool that can handle the so-called element constraint.