在 Pyomo 的 objective 函数中使用分段函数

Question

我正在尝试在 objective 函数中使用 Pyomo 的分段线性函数。这个分段线性函数实际上是对一个名为 macc 的值数组进行插值，该数组具有 401 个值（macc[i]，i 从 0 到 400）。您可以在附图中看到macc的值

我的 objective 函数正在寻找 i 值 macc[i] 遵守约束。为此，我对数组 macc 进行插值以使其具有连续函数 f。见下文：

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)

但是当我尝试在上面的 objective 函数中调用此函数 f 时，我在 objective 函数中收到以下关于表达式 ab = f(m.x) 的消息：

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.

非常欢迎任何关于如何解决该问题的想法。

如果需要，这里是完整的代码。对于这个例子，我用一个函数创建了数组 macc，但实际上它不是来自函数而是来自内部数据。

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

Answer 1

@RonB

正如 AirSquid 评论的那样，您正在使用 kernel 和 environ 命名空间。您应该避免这种混合，因为几种方法可能不兼容。

您可以使用输入、输出参数（在环境层中称为 xvar、yvar) 在定义的变量中输出评估。

使用环境层，分段函数在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)

在这种建模方法中，您不会遇到在每个方程（约束或 objective）中计算 model.piecewise(model.x) 的问题，而是您将只使用 model.y，这相当于评估。

现在，我不知道你的问题，但我猜你的 objective 不是凸的，这可能是优化中的进一步问题。您可以使用 Gurobi 来解决此类问题，但在这种情况下，由于 model.y 取决于 model.x 并且 model.x 是有界的，因此会到达 model.x 上限绑定以使 objective 尽可能低（因为您没有在 objective 中声明任何意义，我假设您想要最小化）。我认为你应该检查你的 objective 是否代表你的想法。

你的 objective 函数正在做这样的事情

在 Pyomo 的 objective 函数中使用分段函数

Using piecewise function in objective function in Pyomo

python

pyomo