Pyomo calling functions when creating constraints

Question

I'm trying to use a loop that creates constraints, where within each constraint, 2 functions are called that perform simple operations on some pyomo variables at a particular index, and then return a vector. An example code is as follows:

import pyomo.environ as pyo

def f(x, y):
    x0 = 2.0*x[1]
    x1 = 0.5*y
    return x0, x1

def g(x):
    x0 = x[0] - 5.0
    x1 = x[1] + 10.0
    return x0, x1

model = pyo.ConcreteModel()

N = 100
num_rows = range(2)
num_cols = range(N)

# Creating decision variables
model.x = pyo.Var(num_rows, num_cols, domain=pyo.Reals)
model.y = pyo.Var(num_cols, domain=pyo.Reals)

# Creating constraints in a loop
model.constraints = pyo.ConstraintList()
for k in range(N-1):
    model.constraints.add(expr=model.x[:,k+1] == model.x[:,k] + 0.5*f(model.x[:,k+1], model.y[k+1]) + 2.0*g(model.x[:,k]))

In the above example, the pyomo variables x and y at index k+1 are passed into the function f. This function then assigns x[0,k+1] = 2.0*x[1,k+1], and x[1,k+1] = 0.5*y[k+1]. The vector x[:,k+1] of dimension 2x1 is then returned by f.

Similarly, for the function g, the pyomo variable x at index k is passed into it, where we then assign x[0,k] = x[0,k] - 5.0 and x[1,k] = x[1,k] + 10.0. The vector x[:,k] of dimension 2x1 is then returned by g.

Currently the above code does not work, as I get the following error:

    x0 = 2.0*x[1]
TypeError: unsupported operand type(s) for *: 'float' and 'IndexedComponent_slice'

Any help on how to use Python functions that manipulate and return pyomo variables would be appreciated. Thanks very much.

Answer 1

When making constraints, I don't believe you can use slices or "vectorize" like you are trying to do. But your construct is crashing before the constraint is made...

I also don't think it is possible to index into a slice of a variable like you are trying to do in your functions... You are thinking "too numpy" here.

Stick with simple expressions for constraints during model construction. So, just break this out into two equality constraints within the loop, and you can probably avoid the complication of the helper function along the way...

# Creating constraints in a loop
model.constraints = pyo.ConstraintList()
for k in range(N-1):
    model.constraints.add(expr=model.x[0, k+1] == ...
    model.constraints.add(expr=model.x[1, k+1] == ...