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

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

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?

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.

• Didn't you mean $\sum_i a_i\cdot \delta_i = x$ ? Feb 20 at 12:46
• No. Can you explain why you think that? Feb 20 at 12:55
• I've done it by introducing binary,, and done by use of element constraint (in YALMIP, but not called that, which does the linearization under the hood)), but was not aware of the term "element constraint", which I guess comes from the constraint programming world. (see cs.brown.edu/courses/csci2580/HTML/lecture2-cp.pdf ) Feb 20 at 13:10
• @ErwinKalvelagen Oh okay, I was a bit confused but you're right. I just do not yet see how this translates to the modelling. I also assume $\delta$ is then indexed on both $i$ and $N$ for my problem? Feb 20 at 13:32
• @Steven01123581321 YALMIP is an optimization modeling system under MATLAB which calls various solvers, and is an example of what Erwin was referring to. It has some (but not unlimited) capability of handling element constraints by automatically reformulating "under the hood" before calling a solver. I think the biggest benefit of the systems which can automatically reformulate is for more complicated models in which the manual reformulation would be very complicated and perhaps error-prone. Feb 20 at 13:49