Pyomo抽象模型结合树基回归模型时约束定义报错求助

你遇到的错误核心是：Pyomo的决策变量（如r[1]）是符号变量，并非普通Python浮点数。线性回归、Ridge这类模型的predict方法本质是线性组合（矩阵乘法+截距），Pyomo符号变量支持重载这类运算符，因此能正常生成约束表达式；但KNN、随机森林这类树基模型的predict逻辑依赖样本距离计算、树分支判断，内部会强制将输入转换为浮点数，而Pyomo默认禁止这种隐式转换，因此触发报错。

方案1：用Pyomo黑箱非线性接口包装树模型 #

适合搭配非线性求解器（如Ipopt），将树模型的预测逻辑包装为Pyomo可识别的外部函数。示例代码如下：

from pyomo.environ import *
from pyomo.opt import SolverFactory
import pandas as pd
import joblib

# 加载树模型
filename = 'model_DecisionTree.pkl'
estimator = joblib.load(filename)

# 定义黑箱预测函数：输入为浮点数列表，返回模型预测结果
def tree_predict(input_list):
    return estimator.predict([input_list])[0]

x_inp = [0.8101265822784809,
        0.802431611347517,
        1.145142857599998,
        0.5823611752600839,
        0.7556756748240172,
        0.37747819179144226,
        0.6501687291338578,
        0.3826800178066474,
         0.7857142833333342,
         0.30769230821005866,
         0.4729729730277574,
         0.1079607425922493,
         0.2260967374953724,
         -0.08108108128560954,
         0.5316094312055879,
         0.12634964888607925,
         -0.3808357060152776]

out_range = pd.read_csv('Output_range.csv', header=0)
inp_range = pd.read_csv('Input_ratio_range.csv', header=0)
n_p = 6

model = AbstractModel()
model.m = Param(within=NonNegativeIntegers)
model.n = Param(within=NonNegativeIntegers)
model.I = RangeSet(1, model.m)
model.J = RangeSet(1, model.n)
model.p = Param(model.I)
model.r = Var(model.I, domain=Reals)

# 目标函数
model.OBJ = Objective(expr=summation(model.p, model.r))

# 定义外部函数：输入维度为用户输入长度+决策变量数量，输出维度为模型输出数量
model.predict_fn = ExternalFunction(
    input_dimension=len(x_inp)+n_p,
    output_dimension=model.n,
    function=tree_predict
)

# 输出下界约束
def low_constraint_rule(m,j):
    input_vec = x_inp + [m.r[i] for i in range(1, n_p+1)]
    return m.predict_fn(*input_vec)[j-1] >= out_range.iat[0,j-1]

# 输出上界约束
def high_constraint_rule(m,j):
    input_vec = x_inp + [m.r[i] for i in range(1, n_p+1)]
    return m.predict_fn(*input_vec)[j-1] <= out_range.iat[1,j-1]

# 其他约束保持不变
model.AxbConstraint_l = Constraint(model.J, rule=low_constraint_rule)
model.AxbConstraint_h = Constraint(model.J, rule=high_constraint_rule)
model.AxbConstraint_dv = Constraint(expr=summation(model.r) == 1)
model.AxbConstraint_dv_l = Constraint(model.I, rule=lambda m,i: m.r[i] >= inp_range.iat[0,i-1])
model.AxbConstraint_dv_h = Constraint(model.I, rule=lambda m,i: m.r[i] <= inp_range.iat[1,i-1])

# 加载数据并使用非线性求解器Ipopt求解
data = DataPortal()
data.load(filename='values.dat', model=model)
instance = model.create_instance(data)
solver = SolverFactory('ipopt')  # 替代GLPK，GLPK仅支持线性问题
results = solver.solve(instance)
instance.display()
print(results)

方案2：将树模型转换为MILP约束 #

若必须使用线性求解器（如GLPK），可借助skompiler库将树模型编译为混合整数线性规划（MILP）约束，示例如下：

from skompiler import skompile
import pyomo.environ as pyo

# 加载树模型
estimator = joblib.load('model_DecisionTree.pkl')

# 将树模型编译为Pyomo兼容的表达式
predict_expr = skompile(estimator.predict, backend='pyomo')

# 修改约束规则，使用编译后的表达式
def low_constraint_rule(m,j):
    input_vec = x_inp + [m.r[i] for i in range(1, n_p+1)]
    return predict_expr(*input_vec)[j-1] >= out_range.iat[0,j-1]

def high_constraint_rule(m,j):
    input_vec = x_inp + [m.r[i] for i in range(1, n_p+1)]
    return predict_expr(*input_vec)[j-1] <= out_range.iat[1,j-1]

注意：复杂树模型（如随机森林）生成的MILP约束数量会非常多，可能导致求解速度变慢。

方案3：迭代优化框架 #

若上述方法不适用，可采用「预测-优化」循环迭代：

1. 初始化决策变量初始值（如均匀分布）
2. 用当前决策变量值输入树模型，得到预测输出，生成线性约束
3. 用GLPK求解线性优化问题，得到新的决策变量值
4. 重复步骤2-3，直到决策变量收敛（变化量小于设定阈值）

核心代码片段：

initial_r = [1/n_p for _ in range(n_p)]
tolerance = 1e-5
max_iter = 50

for iter in range(max_iter):
    # 用当前决策变量值预测输出
    pred_out = estimator.predict([x_inp + initial_r])[0]
    
    # 创建临时线性模型，基于当前预测生成近似约束（示例需根据实际情况调整）
    model = AbstractModel()
    model.m = Param(within=NonNegativeIntegers)
    model.n = Param(within=NonNegativeIntegers)
    model.I = RangeSet(1, model.m)
    model.J = RangeSet(1, model.n)
    model.p = Param(model.I)
    model.r = Var(model.I, domain=Reals)
    
    model.OBJ = Objective(expr=summation(model.p, model.r))
    # 此处需根据模型梯度或预测值构建线性约束，示例仅作参考
    model.AxbConstraint_l = Constraint(model.J, rule=lambda m,j: m.r[j] >= out_range.iat[0,j-1] - pred_out[j-1])
    model.AxbConstraint_h = Constraint(model.J, rule=lambda m,j: m.r[j] <= out_range.iat[1,j-1] - pred_out[j-1])
    model.AxbConstraint_dv = Constraint(expr=summation(model.r) == 1)
    model.AxbConstraint_dv_l = Constraint(model.I, rule=lambda m,i: m.r[i] >= inp_range.iat[0,i-1])
    model.AxbConstraint_dv_h = Constraint(model.I, rule=lambda m,i: m.r[i] <= inp_range.iat[1,i-1])
    
    # 求解得到新的决策变量值
    data = DataPortal()
    data.load(filename='values.dat', model=model)
    instance = model.create_instance(data)
    solver = SolverFactory('glpk')
    solver.solve(instance)
    new_r = [instance.r[i].value for i in range(1, n_p+1)]
    
    # 检查收敛
    if sum(abs(new_r[i]-initial_r[i]) for i in range(n_p)) < tolerance:
        break
    initial_r = new_r

• 线性求解器（GLPK、CBC）仅支持线性约束，树模型本质是非线性/非凸的，因此必须使用非线性求解器或转换模型为线性形式。
• Pyomo的ExternalFunction需要求解器支持黑箱函数，Ipopt是开源常用选择，部分商业求解器（如Gurobi）也支持该特性。

内容的提问来源于stack exchange，提问作者priya kedia