含if2/if3条件的光伏电池管理问题Gekko求解无可行解求助

问题场景 #

我正在研究两栋配备光伏板与电池的住宅的能量管理优化问题，核心逻辑如下：

• 每时间步决策电池充放电动作（action），动作范围受限；
• 根据充放电动作计算能量平衡，值为正则缺能，为负则盈余；
• 盈余住宅可与缺能住宅进行能量交换，交换成本低于电网购电；
• 每时间步的总供给为各住宅盈余之和，总需求为各住宅缺能之和，住宅间交换价格由供需决定。

我用Gekko实现了模型，但求解遇到问题，以下是可复现代码：

import numpy as np
from gekko import gekko as ge

pv_gen = [
    [
        0.        , 0.        , 0.        , 0.        , 0.        ,
        0.        , 0.21628821, 0.45248806, 0.6605544 , 0.82755066,
        0.94309383, 1.        , 0.99473101, 0.92761446, 0.80282334,
        0.62811655, 0.41435652, 0.1748338 , 0.        , 0.        ,
        0.        , 0.        , 0.        , 0.        
    ],
    [
        0.        , 0.        , 0.        , 0.        , 0.        ,
        0.        , 0.21628821, 0.45248806, 0.6605544 , 0.82755066,
        0.94309383, 1.        , 0.99473101, 0.92761446, 0.80282334,
        0.62811655, 0.41435652, 0.1748338 , 0.        , 0.        ,
        0.        , 0.        , 0.        , 0.
    ],

]

load = [
    [
        0.25091601, 0.26052733, 0.25137809, 0.26614229, 0.25835184,
        0.44062643, 0.86531955, 1.        , 0.44139085, 0.25921018,
        0.24908302, 0.26011672, 0.26001465, 0.26010379, 0.26612532,
        0.50450865, 0.5575299 , 0.55895827, 0.74800929, 0.86155478,
        0.74744062, 0.44534278, 0.31157294, 0.25111228
    ],
    [
        0.63107366, 0.62365488, 0.63318315, 0.62548297, 0.63442512,
        0.72115607, 0.94075205, 1.        , 0.72381292, 0.62741716,
        0.63580871, 0.63004974, 0.62598764, 0.62935597, 0.61852053,
        0.74551924, 0.78166073, 0.77483145, 0.87314643, 0.93239358,
        0.86424222, 0.71593412, 0.93499904, 0.62906957
    ],
]

price = [
    0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.12, 0.18, 0.26, 0.18, 0.12,
    0.12, 0.12, 0.12, 0.18, 0.26, 0.18, 0.12, 0.06, 0.06, 0.06, 0.06,
    0.06, 0.06
]

m = ge(remote=False, name='microgrid')

n_houses = 2
steps = 24

# Variables

# Initializing this way allows to put different ub and lb for each house
battery = [m.Array(f=m.Var, dim=(steps + 1), value=0.1, lb=0.1, ub=0.9) for _ in range(n_houses)] 
action = [m.Array(f=m.Var, dim=(steps), value=0, lb=-0.8, ub=0.8) for _ in range(n_houses)]

net_energy = m.Array(f=m.Var, dim=(n_houses,steps), value=0)
l1_import = m.Array(f=m.Var, dim=(n_houses,steps), value=0, lb=0)
l1_export = m.Array(f=m.Var, dim=(n_houses,steps), value=0, lb=0)

# Parameters

surplus_energy = m.Array(f=m.Var, dim=(n_houses,steps), value=0)
shortage_energy = m.Array(f=m.Var, dim=(n_houses,steps), value=0)
offer = m.Array(f=m.Var, dim=(steps), value=0)
demand = m.Array(f=m.Var, dim=(steps), value=0)
l1_import_rate = m.Array(f=m.Var, dim=(steps), value=0)
l1_export_rate = m.Array(f=m.Var, dim=(steps), value=0)
step_transacted_energy = m.Array(f=m.Var, dim=(steps), value=0)
real_net_energy = m.Array(f=m.Var, dim=(n_houses,steps), value=0)

# Cost function initialization

cost = [m.Array(f=m.Var, dim=steps, value=0) for _ in range(n_houses)]

equations = []

# Define constraints for each variable

for t in range(steps):
    
    # Define offer and demand

    equations.append(offer[t] == m.sum([surplus_energy[house_ix][t] for house_ix in range(n_houses)]))
    equations.append(demand[t] == m.sum([shortage_energy[house_ix][t] for house_ix in range(n_houses)]))

    # Define import and export rates

    baseline = price[t] * (1 - 0.25)
    equations.append(l1_import_rate[t] == m.sum([baseline, 0.001 * demand[t], - 0.001 * offer[t]]))
    equations.append(l1_export_rate[t] == m.sum([l1_import_rate[t], - 0.01]))

    # Transaction contraints

    equations.append(step_transacted_energy[t] == m.if2(condition=(m.sum([offer[t], -demand[t]])), x1=offer[t], x2=demand[t]))

    for house_ix in range(n_houses):

        # Constraints to use all the battery energy in the cycle

        equations.append(battery[house_ix][0] == 0.1)
        equations.append(battery[house_ix][steps] == 0.1)

        # Update battery SoC and net energy after action

        equations.append(
            battery[house_ix][t + 1] == m.if3(
                condition = action[house_ix][t],
                x1 = m.sum([battery[house_ix][t], action[house_ix][t] / 0.9 ]),
                x2 = m.sum([battery[house_ix][t], action[house_ix][t] * 0.9 ]),
            )
        )

        equations.append(net_energy[house_ix][t] == m.if3(
            condition = action[house_ix][t],
            x1 = m.sum([load[house_ix][t], -pv_gen[house_ix][t], action[house_ix][t] / 0.9 ]),
            x2 = m.sum([load[house_ix][t], -pv_gen[house_ix][t], action[house_ix][t] * 0.9 ]),
        ))

        equations.append(l1_export[house_ix][t] == m.if3(
            condition=-(offer[t] * demand[t]),
            x1= step_transacted_energy[t] * (surplus_energy[house_ix][t] / offer[t]),
            x2= 0
        ))

        equations.append(l1_import[house_ix][t] == m.if3(
            condition=-(offer[t] * demand[t]),
            x1= step_transacted_energy[t] * (shortage_energy[house_ix][t] / demand[t]),
            x2= 0
        ))
        
        # Define surplus and shortage energy

        equations.append(surplus_energy[house_ix][t] ==  m.max2(-net_energy[house_ix][t],0))
        equations.append(shortage_energy[house_ix][t] == m.max2(net_energy[house_ix][t],0))

        # Cost

        cost[house_ix][t] = m.if3(
            condition=real_net_energy[house_ix][t],
            x1= real_net_energy[house_ix][t] * (price[t] * 0.25) - l1_export[house_ix][t] * l1_export_rate[t] + l1_import[house_ix][t] * l1_import_rate[t],
            x2= real_net_energy[house_ix][t] * (price[t]) + l1_import[house_ix][t] * l1_import_rate[t],
        )

# Define objective function

objs_cost = np.sum(cost, axis=1)

for house_ix, o_c in enumerate(objs_cost):

    m.Minimize(o_c)

m.Equations(equations)

m.options.SOLVER = 1

m.open_folder()

m.solve(disp=True)

原代码问题分析 #

1. 变量与逻辑混淆：real_net_energy未定义关联约束，surplus_energy等变量的推导逻辑存在循环依赖；
2. 冗余约束：电池初始/终态约束在每个时间步重复添加，导致约束爆炸；
3. 多目标冲突：同时最小化两栋住宅的成本，优化器无法找到统一可行解；
4. 数值风险：surplus_energy/offer等除法在分母为0时会引发求解错误；
5. 求解器选择不当：APOPT（SOLVER=1）对部分非线性逻辑的适配性不如IPOPT。

修正后的代码 #

import numpy as np
from gekko import gekko as ge

pv_gen = [
    [0.        , 0.        , 0.        , 0.        , 0.        ,
     0.        , 0.21628821, 0.45248806, 0.6605544 , 0.82755066,
     0.94309383, 1.        , 0.99473101, 0.92761446, 0.80282334,
     0.62811655, 0.41435652, 0.1748338 , 0.        , 0.        ,
     0.        , 0.        , 0.        , 0.        ],
    [0.        , 0.        , 0.        , 0.        , 0.        ,
     0.        , 0.21628821, 0.45248806, 0.6605544 , 0.82755066,
     0.94309383, 1.        , 0.99473101, 0.92761446, 0.80282334,
     0.62811655, 0.41435652, 0.1748338 , 0.        , 0.        ,
     0.        , 0.        , 0.        , 0.        ]
]

load = [
    [0.25091601, 0.26052733, 0.25137809, 0.26614229, 0.25835184,
     0.44062643, 0.86531955, 1.        , 0.44139085, 0.25921018,
     0.24908302, 0.26011672, 0.26001465, 0.26010379, 0.26612532,
     0.50450865, 0.5575299 , 0.55895827, 0.74800929, 0.86155478,
     0.74744062, 0.44534278, 0.31157294, 0.25111228],
    [0.63107366, 0.62365488, 0.63318315, 0.62548297, 0.63442512,
     0.72115607, 0.94075205, 1.        , 0.72381292, 0.62741716,
     0.63580871, 0.63004974, 0.62598764, 0.62935597, 0.61852053,
     0.74551924, 0.78166073, 0.77483145, 0.87314643, 0.93239358,
     0.86424222, 0.71593412, 0.93499904, 0.62906957]
]

price = [0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.12, 0.18, 0.26, 0.18, 0.12,
         0.12, 0.12, 0.12, 0.18, 0.26, 0.18, 0.12, 0.06, 0.06, 0.06, 0.06,
         0.06, 0.06]

m = ge(remote=False, name='microgrid')

n_houses = 2
steps = 24

#