解决MindtPy中Pyomo非凸INLP的KeyError及优化策略

pip install --upgrade pyomo

from pyomo.environ import SolverFactory

opt = SolverFactory('mindtpy')
results = opt.solve(
    M,
    mip_solver='cplex',  # 或 'gurobi', 'glpk' 等
    nlp_solver='ipopt',  # 或 'bonmin', 'baron' 等
    strategy='GOA',      # 指定使用全局外逼近策略
    tee=True
)

import math # 导入math模块以使用log函数，虽然Pyomo通常会自动处理

def obj(M):
  # 确保对数函数的参数始终大于一个非常小的正数
  # 使用一个小的epsilon值，例如1e-9，以避免log(0)
  epsilon = 1e-9
  return sum(-log(max(epsilon, 1 - prod((1-pyp[i,j,k])**(M.n[i,j,k]) for j in M.e for k in M.s))) for i in M.m)

M.obj = Objective(rule=obj, sense=minimize)

def obj(M):
  # 在Pyomo中直接使用 log(expr + epsilon) 是更常见的做法
  # 假设 1 - prod(...) 的理论最小值是0，添加epsilon确保参数 > 0
  epsilon = 1e-9
  return sum(-log(1 - prod((1-pyp[i,j,k])**(M.n[i,j,k]) for j in M.e for k in M.s) + epsilon) for i in M.m)
M.obj = Objective(rule=obj, sense=minimize)

import pandas as pd
import random as r
import numpy as np
# glpk 通常用于MILP，这里MindtPy会调用其MIP求解器
from pyomo.environ import *
# amplpy.AMPL 如果不直接使用AMPL求解器，可能不是必需的

def pyblock(pyp, pytau, pyr, pys):
  M = ConcreteModel()

  M.m = Set(initialize = list(range(int(len(pyp)))))
  M.e = Set(initialize = list(range(int(len(pyr)))))
  M.s = Set(initialize = list(range(int(pys))))

  M.r = Param(M.e, initialize = pyr)
  M.tau = Param(M.m, initialize = pytau)
  # p 参数的初始化方式可能需要根据实际数据调整，确保与M.m, M.e, M.s的索引匹配
  # 如果pyp已经是正确格式的字典或列表，直接传递可能更合适
  # 示例中pyp[i,j,k]的访问方式暗示pyp是一个多维数组或字典
  M.p = Param(M.m, M.e, M.s, initialize = 0) # 假设pyp会在目标函数中直接使用

  M.n = Var(M.m, M.e, M.s, domain=NonNegativeIntegers, initialize=0)

  def obj_rule(M):
    epsilon = 1e-9 # 添加一个小的正数以避免log(0)
    # 确保pyp的索引方式与M.m, M.e, M.s一致
    # 这里的pyp[i,j,k]假设是一个外部传入的字典或列表，与M.p不同
    return sum(-log(1 - product((1-pyp[i,j,k])**(M.n[i,j,k]) for j in M.e for k in M.s) + epsilon) for i in M.m)
  M.obj = Objective(rule=obj_rule, sense=minimize)

  def fire_rate_rule(M, j, k):
    return sum(M.n[i,j,k] for i in M.m) <= M.r[j]
  M.fire_rate = Constraint(M.e, M.s, rule = fire_rate_rule)

  opt = SolverFactory('mindtpy')
  results = opt.solve(
    M,
    mip_solver = 'cplex', # 确保您的系统已安装并配置了cplex
    nlp_solver = 'ipopt', # 确保您的系统已安装并配置了ipopt
    strategy = 'GOA',     # 针对非凸INLP问题使用GOA策略
    tee=True
  )

  # 检查求解状态
  if (results.solver.status == SolverStatus.ok) and \
     (results.solver.termination_condition == TerminationCondition.optimal):
      print("MindtPy 求解成功，找到最优解。")
  elif results.solver.termination_condition == TerminationCondition.infeasible:
      print("模型不可行。")
  else:
      print(f"MindtPy 求解器状态：{results.solver.status}, 终止条件：{results.solver.termination_condition}")

  return M.n.extract_values()

# 示例数据（请根据实际情况提供）
# 假设pyp是一个字典，键为(i,j,k)
_pyp_data = {
    (0,0,0): 0.1, (0,0,1): 0.2, (0,0,2): 0.15,
    (0,1,0): 0.05, (0,1,1): 0.1, (0,1,2): 0.08,
    (1,0,0): 0.2, (1,0,1): 0.1, (1,0,2): 0.25,
    (1,1,0): 0.15, (1,1,1): 0.2, (1,1,2): 0.1
}
# 假设pyp是一个函数或lambda表达式，根据(i,j,k)返回对应值
_pyp_func = lambda i,j,k: _pyp_data.get((i,j,k), 0) # 确保有默认值

_pytau_data = {0: 10, 1: 15}
_pyr_data = {0: 5, 1: 8}
_pys_val = 3

# 运行优化
# optimal_n_values = pyblock(_pyp_func, _pytau_data, _pyr_data, _pys_val)
# print(optimal_n_values)