在PYOMO约束中使用外部功能

Question

问题摘要

是PYOMO的新手，我正在尝试解决一个非线性优化问题，其中约束是任意的Python函数，其值必须在Intervel [1/8，1]中。就我而言，这些功能是scipy.interpaly.pchipinterpolator s，并且由于问题给了它们，所以我不能太多更改它们。

我知道 s，但是我必须做错了什么，因为我得到了keyError：'pyomo_socket_server'。据我所知：

• f可以是任何python callable。
• g是梯度向量。就我而言，所有功能都是r-＆gt; r，这是f的第一个导数的长度1列表。
• H是Hessian矩阵，但我们只能返回上部三角形部分。就我而言，所有功能都是r-＆gt; r，这是f的第二个导数的长度1列表。

我正在使用ipopt，它可以正常工作，因为我尝试了使用非线性约束的其他问题（但没有外部功能 s）。

您发现任何明显的错误吗？我该如何完成这项工作？

最小（非）工作示例

import numpy as np
from pyomo.environ import (
    Any,
    ConcreteModel,
    ExternalFunction,
    NonNegativeIntegers,
    NonNegativeReals,
    Param,
    Set,
    SolverFactory,
    Var,
    maximize,
    summation,
)
from scipy.interpolate import PchipInterpolator

# Generate example data
bounds = {
    "g": (800, 2150),
    "v": (1200, 2150),
    "em": (1300, 2350),
    "e": (1400, 2550),
    "ep": (1500, 2450),
    "d": (1800, 2650),
}

y_bounds = {
    "g": (0, 2),
    "v": (0, 4),
    "em": (0, 3.5),
    "e": (0, 3.5),
    "ep": (0, 1.75),
    "d": (0, 1.5),
}

curves = {}
for (product_, (lower_x, upper_x)), (lower_y, upper_y) in zip(bounds.items(), y_bounds.values()):
    x = np.arange(lower_x, upper_x, 1)
    y = (1 - 1 / (1 + np.exp(-0.01 * (x - x.mean())))) * (upper_y - lower_y) + lower_y
    curves[product_] = PchipInterpolator(x, y)

curves_derivatives = {product_: curve.derivative(nu=1) for product_, curve in curves.items()}
curves_derivatives_2 = {product_: curve.derivative(nu=2) for product_, curve in curves.items()}

inventory = dict(zip(bounds, (1, 4, 6, 10, 4, 1)))
initial_guesses = {k: lower_x for k, (lower_x, _) in bounds.items()}

# Create Pyomo model
opt_model = ConcreteModel(name="optimizer")

opt_model.item_ids = Set(initialize=bounds.keys())
opt_model.inventory = Param(
    opt_model.item_ids, initialize=inventory, domain=NonNegativeIntegers, mutable=False
)

opt_model.curves = Param(opt_model.item_ids, initialize=curves, mutable=False, domain=Any)

opt_model.curves_first_derivative = Param(
    opt_model.item_ids, initialize=curves_derivatives, domain=Any, mutable=False
)

opt_model.curves_second_derivative = Param(
    opt_model.item_ids, initialize=curves_derivatives_2, domain=Any, mutable=False
)


opt_model.x = Var(
    opt_model.item_ids, bounds=bounds, initialize=initial_guesses, domain=NonNegativeReals,
)


@opt_model.Objective(sense=maximize)
def objective_rule(model):
    return summation(model.x)


# Difficult constraint.
@opt_model.Constraint(opt_model.item_ids)
def my_rule(model, item_id):
    def f(x):
        return model.curves[item_id](x) / model.inventory[item_id]

    def g(args, fixed=None):
        return [model.curves_first_derivative[item_id](args[0]) / model.inventory[item_id]]

    def h(args, fixed=None):
        return [model.curves_second_derivative[item_id](args[0]) / model.inventory[item_id]]

    external_function = ExternalFunction(function=f, gradient=g, hessian=h)
    return (1 / 8, external_function(model.x[item_id]), 1)


SolverFactory("ipopt").solve(opt_model)