Search code examples
pythonpyomo

Using piecewise function in objective function in Pyomo

I am trying to use a piecewise linear function from Pyomo in my objective function. This piecewise linear function is actually interpolating an array of values called macc, having 401 values (macc[i], i from 0 to 400). You can see the values of macc in the picture attached

enter image description here

My objective function is looking for the value i for which macc[i] respect the constraints. To do that I interpolate the array macc to have a continuous function f. See below:

c = np.arange(401)
f = pyopiecewise.piecewise(c,macc,validate=False)
model = pyo.ConcreteModel()

#Declare variable
model.x = pyo.Var(domain=pyo.NonNegativeReals, bounds=(5,395), initialize = cp0)

#Declare parameters
model.s = pyo.Param(domain=pyo.NonNegativeReals,initialize=s0)
model.b = pyo.Param(domain=pyo.NonNegativeReals,initialize=b0)
model.tnac = pyo.Param(domain=pyo.NonNegativeReals,initialize=tnac0)

#Objective function
def objective_(m):
    ab = f(m.x)

    e = m.b - ab

    return (e * m.x)

#Constraints
def constraint1(m):

    ab = f(m.x)

    e = m.b - ab

    return e <= (m.tnac + m.s)

But when I try to call this function f in my objective function above, I get the message below for the expression ab = f(m.x) in the objective function:

ERROR: Rule failed when generating expression for Objective Obj with index
    None: PyomoException: Cannot convert non-constant expression to bool. This
    error is usually caused by using an expression in a boolean context such
    as an if statement. For example,
        m.x = Var() if m.x <= 0:
        ...
would cause this exception.

ERROR: Constructing component 'Obj' from data=None failed: PyomoException:
    Cannot convert non-constant expression to bool. This error is usually
caused by using an expression in a boolean context such as an if
statement. For example,
        m.x = Var() if m.x <= 0:
        ...
would cause this exception.

Any idea on how to solve that would be very welcome.

Here is the complete code if needed. For this example I created the array macc with a function, but in reality it does not come from a function but from internal data.

import numpy as np
import pyomo.environ as pyo
import pyomo.core.kernel.piecewise_library.transforms as pyopiecewise

#Create macc
# logistic sigmoid function
def logistic(x, L=1, x_0=0, k=1):
    return L / (1 + np.exp(-k * (x - x_0)))


c = np.arange(401)
macc = 2000*logistic(c,L=0.5,x_0 = 60,k=0.02)
macc = macc -macc[0]

f = pyopiecewise.piecewise(c,macc,validate=False)

s0 = 800
b0 = 1000
tnac0 = 100

cp0 = 10
ab0 = 100

model = pyo.ConcreteModel()

#Declare variable
model.x = pyo.Var(domain=pyo.NonNegativeReals, bounds=(5,395), initialize = cp0)

#Declare parameters
model.s = pyo.Param(domain=pyo.NonNegativeReals,initialize=s0)
model.b = pyo.Param(domain=pyo.NonNegativeReals,initialize=b0)
model.tnac = pyo.Param(domain=pyo.NonNegativeReals,initialize=tnac0)

#Objective function
def objective_(m):
    ab = f(m.x)

    e = m.b - ab

    return (e * m.x)

model.Obj = pyo.Objective(rule=objective_)

#Constraints
def constraint1(m):

    ab = f(m.x)

    e = m.b - ab

    return e <= (m.tnac + m.s)

def constraint2(m): 

    ab = f(m.x)

    e = m.b - ab

    return e >= 1

def constraint3(m):

    ab = f(m.x)

    return ab >= 0


model.con1 = pyo.Constraint(rule = constraint1)
model.con2 = pyo.Constraint(rule = constraint2)
model.con3 = pyo.Constraint(rule = constraint3)

And here is my objective function

Solution

  • @RonB

    As AirSquid commented, you're using the kernel and environ namespace. You should avoid this mixing, since several approach may not be compatible.

    Instead of explicit evaluating the piecewised funtion using the __call__() (f(model.x)) method you can use the input, output args (in environ layer are called xvar, yvar) to output the evaluation in a defined variable.

    Using the environ layer, piecewise function are available in pyo.Piecewise

    import numpy as np
import pyomo.environ as pyo

#Create macc
# logistic sigmoid function
def logistic(x, L=1, x_0=0, k=1):
    return L / (1 + np.exp(-k * (x - x_0)))

c = np.linspace(0,400,400)
macc = 2000*logistic(c,L=0.5,x_0 = 60,k=0.02)
macc = macc -macc[0]

s0 = 800
b0 = 1000
tnac0 = 100

cp0 = 10
ab0 = 100

model = pyo.ConcreteModel()

#Declare variable
model.x = pyo.Var(domain=pyo.NonNegativeReals, bounds=(5,395), initialize = cp0)
model.y = pyo.Var()
model.piecewise = pyo.Piecewise(model.y, model.x, pw_pts=list(c), f_rule=list(macc), pw_constr_type='EQ', pw_repn='DCC')

#Declare parameters
model.s = pyo.Param(domain=pyo.NonNegativeReals,initialize=s0)
model.b = pyo.Param(domain=pyo.NonNegativeReals,initialize=b0)
model.tnac = pyo.Param(domain=pyo.NonNegativeReals,initialize=tnac0)


model.Obj = pyo.Objective(expr= model.b*model.x - model.y*model.x, sense=pyo.minimize)

model.con1 = pyo.Constraint(expr=model.b - model.y <= model.tnac + model.s)
model.con2 = pyo.Constraint(expr=model.b - model.y >= 1)
model.con3 = pyo.Constraint(expr= model.y >= 0)
solver = pyo.SolverFactory('ipopt')
solver.solve(model, tee=True)

    In this modeling approach you don't have the problem of evaluate model.piecewise(model.x) in each equation (constraint or objective), instead you will just use model.y which is equivalent to evaluate.

    Now, I don't know your problem, but that I guess your objective is not convex, which may be a further problem in the optimization. you can use Gurobi to solve such problems, but in this case, as model.y depends upon model.x and model.x is bounded, is going to the model.x upper bound in order to make objective as low as possible (since you don't declare any sense in the objective, I assume that you want to minimize). I think you should check your objective if it represents what you think.

    Your objective function is doing something like this objective-function-w.r.t-xvar