如何在Pyomo中针对黑箱外部仿真实现无梯度优化？

1. 排查Mindtpy的配置问题 #

Mindtpy默认可能依赖梯度信息，需强制切换到无梯度适配模式：

• 显式指定全局外逼近策略（GLOA），搭配无梯度下层求解器（如Scipy的nelder-mead或COBYLA），示例配置：

from pyomo.environ import ConcreteModel
from mindtpy import MindtPySolver

model = ConcreteModel()
# 定义变量、约束及调用黑箱仿真的目标函数
...

solver = MindtPySolver(model)
solver.solve(
    strategy='GLOA',
    nlp_solver='scipy',
    nlp_solver_args={'method': 'nelder-mead'},
    mip_solver='cbc',
    max_iter=200
)

• 检查变量边界是否贴合仿真实际可行域：不合理的上下限会导致求解器搜索方向偏离，需根据仿真的输入限制设置model.x.setlb(...)和model.x.setub(...)。

2. 切换至专用无梯度求解器 #

如果Mindtpy适配不佳，直接用Pyomo对接成熟的无梯度优化工具：

• Scipy无梯度求解器：通过SolverFactory('scipy')调用，支持多种无梯度算法，示例：

from pyomo.environ import SolverFactory

solver = SolverFactory('scipy')
result = solver.solve(
    model,
    method='COBYLA',  # 适合带约束的无梯度优化
    options={'maxiter': 1000, 'rhobeg': 0.1}
)

• NLopt全局无梯度求解器：支持DIRECT、粒子群等全局算法，适合黑箱问题的全局寻优，示例：

from pyomo.environ import SolverFactory

solver = SolverFactory('nlopt')
solver.solve(
    model,
    algorithm='GN_DIRECT_L',  # 全局DIRECT算法
    options={'maxeval': 500}
)

3. 优化黑箱仿真的调用效率 #

无梯度优化迭代次数多，需降低仿真耗时：

• 缓存已计算的输入-输出结果：用字典记录已运行过的input.json对应的目标函数值，避免重复仿真。
• 变量降维：通过敏感性分析筛选对目标函数影响显著的变量，减少优化维度。
• 代理模型预搜索：先用响应面、高斯过程等代理模型快速定位最优区域，再用真实仿真精细优化。

4. 自定义无梯度优化循环 #

若现有求解器无法满足需求，可手动实现启发式优化逻辑（如遗传算法），伪代码示例：

import random
from pyomo.environ import value

# 初始化变量种群
pop_size = 20
lb, ub = 0, 10  # 变量边界
population = [random.uniform(lb, ub) for _ in range(pop_size)]

max_iter = 50
best_obj = float('inf')
best_x = None

for _ in range(max_iter):
    fitness = []
    for x in population:
        # 生成input.json
        with open('input.json', 'w') as f:
            f.write(f'{{"variable": {x}}}')
        # 运行外部仿真（假设仿真命令为run_simulation.exe）
        import subprocess
        subprocess.run(['run_simulation.exe'], check=True)
        # 读取output.json计算目标函数
        import json
        with open('output.json', 'r') as f:
            output = json.load(f)
        obj_val = output['target_param']
        fitness.append(obj_val)
    
    # 更新最优解
    current_best_idx = fitness.index(min(fitness))
    if fitness[current_best_idx] < best_obj:
        best_obj = fitness[current_best_idx]
        best_x = population[current_best_idx]
    
    # 遗传算法选择、交叉、变异生成新种群
    # 示例：保留前50%最优个体，其余通过交叉变异生成
    sorted_pop = [x for _, x in sorted(zip(fitness, population))]
    new_pop = sorted_pop[:pop_size//2]
    for _ in range(pop_size//2):
        parent1 = random.choice(sorted_pop[:pop_size//4])
        parent2 = random.choice(sorted_pop[:pop_size//4])
        child = (parent1 + parent2) / 2 + random.uniform(-0.5, 0.5)
        child = max(lb, min(ub, child))  # 约束在边界内
        new_pop.append(child)
    population = new_pop

print(f"最优解：x={best_x}, 目标函数值={best_obj}")

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