GEKKO求解INLP遇Solution Not Found错误，求可行重构方案

问题根源分析 #

1. 非光滑函数导致求解困难：原模型使用m.max2()/m.min2()这类非光滑函数，与整数变量结合时，容易造成求解器在搜索可行解时陷入局部最优或无法收敛。
2. 设备运行逻辑不完整：原代码仅考虑了启动后1小时的设备负荷，未覆盖设备运行时长l的全部时段，导致约束逻辑错误。
3. 变量定义与目标函数矛盾：原E_gain被定义为负数（m.min2(0, ...)的结果），但售电量应为非负值，导致目标函数计算逻辑错误，进一步引发约束冲突。

具体修复步骤 #

1. 用线性约束替换非光滑max/min函数 #

将购售电的非线性约束转化为等价线性约束，消除非光滑性：

• 购电量E_spent[i]（非负）与售电量E_gain[i]（非负）满足：power_balance[i] + device_contribution[i] = E_spent[i] - E_gain[i]
• 该约束自动实现：当总功率（家庭平衡+设备贡献）为正时，E_spent[i]等于总功率，E_gain[i]为0；当总功率为负时，E_gain[i]等于总功率的绝对值，E_spent[i]为0。

2. 修正设备运行时长逻辑 #

通过辅助变量is_running[i]建模设备连续运行状态：

• is_running[i]为二进制变量，表示第i小时设备是否在运行。
• 通过线性组合关联启动变量is_start：is_running[i]等于窗口[max(0, i-l+1), i]内所有is_start的和（因仅能启动一次，结果为0或1）。
• 设备功率贡献device_contribution[i]通过启动时刻与设备功率曲线的线性组合计算，确保启动后连续l小时的负荷都被正确计入。

3. 修正变量定义与目标函数 #

• 将E_gain重新定义为非负的售电量，符合实际经济意义。
• 目标函数调整为购电成本减去售电收益，正确反映总电费的最小化需求。

4. 优化求解器参数（可选） #

适当放宽整数容忍度、增加迭代次数，帮助求解器更高效地搜索可行解。

修正后的完整代码 #

import pandas as pd
import numpy as np
from gekko import GEKKO

# 读取数据
df = pd.read_csv('hourly.csv')
power_balance = df["house_connection"].to_numpy()  # 家庭电网连接功率（可正可负）

df = pd.read_csv('device_hourly.csv')
device_profile = df["device"].to_numpy()  # 设备功率曲线

# 参数定义
N = power_balance.size 
l = device_profile.size  # 设备运行时长，满足 l < N
priceBuy = 0.23  # 购电价（欧元/kWh）
priceSell = 0.063  # 售电价（欧元/kWh）

# 初始化GEKKO模型
m = GEKKO(remote=False) 
m.options.SOLVER = 1  # 使用APOPT求解器处理整数规划

# APOPT求解器配置
m.solver_options = ['minlp_maximum_iterations 2000',
                   'minlp_max_iter_with_int_sol 1000',
                   'minlp_as_nlp 0',
                   'nlp_maximum_iterations 1500',
                   'minlp_integer_tol 0.01',
                   'minlp_gap_tol 0.05']

# 变量定义
is_start = m.Array(m.Var, N, integer=True, lb=0, ub=1)  # 二进制启动变量（1表示该时刻启动）
is_running = m.Array(m.Var, N, integer=True, lb=0, ub=1)  # 设备运行状态变量
device_contribution = m.Array(m.Var, N)  # 设备在第i小时的功率贡献
E_spent = m.Array(m.Var, N, lb=0)  # 购电量（非负）
E_gain = m.Array(m.Var, N, lb=0)  # 售电量（非负）

# 约束条件
# 1. 仅能启动一次设备
m.Equation(sum(is_start) == 1)

# 2. 禁止在无法完成全周期运行的时段启动
m.Equations([is_start[i] == 0 for i in range(N - l + 1, N)])

# 3. 设备运行状态与启动逻辑关联
for i in range(N):
    # 计算当前时刻的启动窗口：[i-l+1, i]（不小于0）
    start_window = range(max(0, i - l + 1), i + 1)
    # 运行状态等于窗口内启动变量的和（仅一个启动变量为1，结果为0或1）
    m.Equation(is_running[i] == sum(is_start[k] for k in start_window))
    # 计算设备功率贡献：启动时刻对应的功率曲线值
    m.Equation(device_contribution[i] == sum(is_start[k] * device_profile[i - k] 
                                             for k in start_window if (i - k) < l))

# 4. 功率平衡约束：家庭功率+设备功率 = 购电量 - 售电量
for i in range(N):
    m.Equation(power_balance[i] + device_contribution[i] == E_spent[i] - E_gain[i])

# 目标函数：最小化总电费（购电成本 - 售电收益）
m.Minimize(sum(E_spent[i] * priceBuy - E_gain[i] * priceSell for i in range(N)))

# 求解模型
m.solve(disp=True)

# 输出结果
print('Results')  
print('总电费开销: ' + str(m.options.objfcnval))
# 提取启动时刻
start_time = np.argmax([var.value[0] for var in is_start])
print(f'设备最佳启动时刻: 第{start_time}小时')

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