求助：解决Gekko滚动优化中‘Solution Not Found’错误

程序背景 #

我开发了一套用于优化家用电池充放电的程序，结合光伏（PV）系统以最小化年度电费。核心逻辑基于净取电判断：

• 当净取电>0时，家庭用电大于光伏发电量，需从电池放电或电网购电；
• 当净取电<0时，家庭用电小于光伏发电量，可为电池充电或向电网售电。

家庭用电数据每15分钟采集一次，每天共96个数据点。程序采用滚动优化策略，每次对连续2天的电池充放电进行优化，让第一天的决策能考虑次日数据。

控制器代码 #

from gekko import GEKKO
from simulationModel2_2d_1 import getSimulation2
from exportModel2 import exportToExcelModel2
import numpy as np
#import matplotlib.pyplot as plt
import pandas as pd
import time
import math

# ------------------------ Import and read input data ------------------------

file = r'Data Sim 2.xlsx'

data = pd.read_excel(file, sheet_name='Input', na_values='NaN')

dataRead = pd.DataFrame(data, columns= ['Timestep','Verbruik woning (kWh)','Netto afname (kWh)','Prijs afname (€/kWh)',
                                            'Prijs injectie (€/kWh)','Capaciteit batterij (kW)',
                                            'Capaciteit batterij (kWh)','Rendement (%)',
                                            'Verbruikersprofiel','Capaciteit PV (kWp)','Aantal dagen'])

timestep = dataRead['Timestep'].to_numpy()                                  
usage_home = dataRead['Verbruik woning (kWh)'].to_numpy()
net_offtake = dataRead['Netto afname (kWh)'].to_numpy()
price_offtake = dataRead['Prijs afname (€/kWh)'].to_numpy()
price_injection = dataRead['Prijs injectie (€/kWh)'].to_numpy()
cap_batt_kW = dataRead['Capaciteit batterij (kW)'].iloc[0]              
cap_batt_kWh = dataRead['Capaciteit batterij (kWh)'].iloc[0]
efficiency = dataRead['Rendement (%)'].iloc[0]
usersprofile = dataRead['Verbruikersprofiel'].iloc[0]
days = dataRead['Aantal dagen'].iloc[0]
pv = dataRead['Capaciteit PV (kWp)'].iloc[0]

# ------------- Optimization model & Rolling principle (2 days) --------------

# Initialise model
m = GEKKO()

# Output data
ts = []
charging = []                       # Amount to charge/decharge batterij
e_batt = []                         # Amount of energy in the battery
usage_net = []                      # Usage after home, battery and pv
p_paid = []                         # Price paid for energy of 15min

# Energy in battery to pass
energy = 0

# Iterate each day for one year
for d in range(int(days)-1):
    d1_timestep = []
    d1_net_offtake = []
    d1_price_offtake = []
    d1_price_injection = []
    d2_timestep = []
    d2_net_offtake = []
    d2_price_offtake = []
    d2_price_injection = []
    
    # Iterate timesteps
    for i in range(96):
        d1_timestep.append(timestep[d*96+i])
        d2_timestep.append(timestep[d*96+i+96])
        
        d1_net_offtake.append(net_offtake[d*96+i])
        d2_net_offtake.append(net_offtake[d*96+i+96])
        
        d1_price_offtake.append(price_offtake[d*96+i])
        d2_price_offtake.append(price_offtake[d*96+i+96])
        
        d1_price_injection.append(price_injection[d*96+i])
        d2_price_injection.append(price_injection[d*96+i+96])
    
    # Input data simulation of 2 days
    ts_temp = np.concatenate((d1_timestep, d2_timestep))
    net_offtake_temp = np.concatenate((d1_net_offtake, d2_net_offtake))
    price_offtake_temp = np.concatenate((d1_price_offtake, d2_price_offtake))
    price_injection_temp = np.concatenate((d1_price_injection, d2_price_injection))
    if(d == 7):
        print(ts_temp)
        print(energy)
    # Simulatie uitvoeren
    charging_temp, e_batt_temp, usage_net_temp, p_paid_temp, energy_temp = getSimulation2(ts_temp, net_offtake_temp, price_offtake_temp, price_injection_temp,                                                                                     cap_batt_kW, cap_batt_kWh, efficiency, energy, pv)
    
    # Take over output first day, unless last 2 days
    energy = energy_temp
    if(d == (days-2)):
        for t in range(1,len(ts_temp)):
            ts.append(ts_temp[t])
            charging.append(charging_temp[t])
            e_batt.append(e_batt_temp[t])
            usage_net.append(usage_net_temp[t])
            p_paid.append(p_paid_temp[t])
    elif(d == 0):
        for t in range(int(len(ts_temp)/2)+1):
            ts.append(ts_temp[t])
            charging.append(charging_temp[t])
            e_batt.append(e_batt_temp[t])
            usage_net.append(usage_net_temp[t])
            p_paid.append(p_paid_temp[t])
    else:
        for t in range(1,int(len(ts_temp)/2)+1):
            ts.append(ts_temp[t])
            charging.append(charging_temp[t])
            e_batt.append(e_batt_temp[t])
            usage_net.append(usage_net_temp[t])
            p_paid.append(p_paid_temp[t])
            
    print('Simulation day '+str(d+1)+' complete.')

# ------------------------ Export output data to Excel -----------------------

a = exportToExcelModel2(ts, usage_home, net_offtake, price_offtake, price_injection, charging, e_batt, usage_net, p_paid, cap_batt_kW, cap_batt_kWh, efficiency, usersprofile, pv)
print(a)

Gekko优化模型代码 #

from gekko import GEKKO

def getSimulation2(timestep, net_offtake, price_offtake, price_injection, 
                   cap_batt_kW, cap_batt_kWh, efficiency, start_energy, pv):
    
    # ---------------------------- Optimization model ----------------------------
       
    # Initialise model
    m = GEKKO(remote = False)
    
    # Global options
    m.options.SOLVER = 1
    m.options.IMODE  = 6
    
    # Constants
    speed_charging = cap_batt_kW/4
    m.time = timestep
        
    max_cap_batt = m.Const(value = cap_batt_kWh)
    min_cap_batt = m.Const(value = 0)
    max_charge = m.Const(value = speed_charging)                                    # max battery can charge in 15min
    max_decharge = m.Const(value = -speed_charging)                                 # max battery can decharge in 15min
    
    # Parameters
    usage_home = m.Param(net_offtake)
    price_offtake = m.Param(price_offtake)
    price_injection = m.Param(price_injection)
    
    # Variables
    e_batt = m.Var(value=start_energy, lb = min_cap_batt, ub = max_cap_batt)  # energy in battery
    price_paid = m.Var()                                              # price paid each 15min
    charging = m.Var(lb = max_decharge, ub = max_charge)              # amount charge/decharge each 15min
    usage_net = m.Var(lb=min_cap_batt)
    
    # Equations
    m.Equation(e_batt==(m.integral(charging)+start_energy)*efficiency)
    m.Equation(-charging <= e_batt)
    m.Equation(usage_net==usage_home + charging)
    price = m.Intermediate(m.if2(usage_net*1e6, price_injection, price_offtake))
    price_paid = m.Intermediate(usage_net * price / 100)
    
    # Objective 
    m.Minimize(price_paid)
    
    # Solve problem
    m.options.COLDSTART=2
    m.solve()
    
    m.options.TIME_SHIFT=0
    m.options.COLDSTART=0
    m.solve()
          
    # Energy to pass
    energy_left = e_batt[95]
    
    #m.cleanup()
    
    return charging, e_batt, usage_net, price_paid, energy_left

问题描述 #

程序运行至第17天时，持续出现「Solution Not Found」错误。已尝试以下优化手段但均无效：

• 将默认迭代上限提升至500
• 更换求解器
• 启用COLDSTART预求解（未使用预求解时，程序在第8天就会报错）

但单独求解出错的天数时，相同代码能立即得到解。

原因分析 #

1. 滚动优化的状态传递矛盾：前一天优化结束后传递的电池剩余能量energy_left，可能在后续两天的优化场景中触发约束冲突。比如当传递的初始电量接近电池上下限，结合当天的用电/发电数据，模型无法找到满足所有约束的可行解；但单独求解时初始电量为合理值（如0），因此能正常求解。
2. m.if2函数的阈值设置不合理：代码中使用usage_net*1e6作为if2的触发阈值，过大的数值可能导致边界场景下的逻辑判断失效，引发求解器的数值不稳定。
3. 积分方程的累计误差：电池能量的积分方程e_batt==(m.integral(charging)+start_energy)*efficiency在滚动迭代中可能产生累计误差，导致约束条件无法满足。
4. 求解流程的状态不一致：两次solve()调用（COLDSTART预求解和正常求解）的衔接可能存在状态残留，尤其是在滚动迭代中，前一次求解的状态会干扰下一次的初始化。

解决方案 #

1. 松弛状态传递的边界约束 #

在传递电池初始能量时，添加微小松弛量，避免初始值卡在上下限：

# 在调用getSimulation2前修改energy值
energy = max(min(energy, cap_batt_kWh - 0.01), 0.01)

2. 修正if2的逻辑判断 #

直接基于usage_net的符号判断购电/售电价，避免放大数值引发的问题：

# 替换原price的Intermediate定义
price = m.Intermediate(m.if2(usage_net, price_injection, price_offtake))

注：m.if2(a,b,c)的逻辑是当a>0时返回b，否则返回c，这里usage_net>0对应购电（取price_offtake），usage_net<0对应售电（取price_injection），符合业务逻辑。

3. 优化积分方程写法 #

改用Gekko推荐的动态建模方式，明确导数关系，减少累计误差：

# 替换原积分方程
m.Equation(e_batt.dt() == charging * efficiency)
m.fix_initial(e_batt, val=start_energy)

4. 简化求解流程并增加迭代次数 #

合并求解步骤，确保状态重置，同时增加迭代上限：

# 替换原两次solve调用
m.options.COLDSTART = 2
m.options.MAX_ITER = 500
m.solve(disp=False)
m.options.COLDSTART = 0
m.options.MAX_ITER = 500
m.solve(disp=False)

disp=False可减少输出干扰，便于观察迭代过程。

5. 验证异常批次的初始状态 #

打印第16天结束后的energy_left值，检查是否处于异常范围（如负数或超过电池容量）：

# 在循环中添加打印
print(f"Day {d+1} end energy: {energy_temp}")

若发现异常值，回溯排查前几天的求解结果，确认是否存在隐性错误。

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