Gekko中B_SOC动态更新及二进制变量I_disch求解异常咨询

问题1：B_SOC动态更新异常 #

问题描述 #

用户新增B_SOC_test变量后，发现电池充放电时只有B_SOC_test能正常动态更新，而B_SOC被赋值为B_Bound后无法实现动态更新。需求是让B_SOC始终取计算值与电池容量的最小值，且后续计算基于最新的B_SOC。

原因分析 #

1. 代码中先将B_SOC定义为GEKKO优化变量，后续又通过B_SOC=[B_Bound[t] for t in range(row)]将其重新赋值为B_Bound的列表，覆盖了原有的优化变量属性，导致B_SOC不再参与优化过程。
2. 原有的B_Bound约束仅设置了B_Bound[t]<=B_SOC[t]和B_Bound[t]<=Battery_capacity，没有实现“取计算值与容量最小值”的逻辑。

解决方法 #

直接对B_SOC设置上限约束，或者用GEKKO的m.min2函数定义B_Bound为B_SOC计算值与容量的最小值，再让B_SOC等于B_Bound，保证动态更新：

# 替换原B_Bound相关约束
for t in range(row):
    # 用min2实现取最小值：B_Bound[t] = min(B_SOC_calc, Battery_capacity)
    m.Equation(B_Bound[t] == m.min2(B_SOC[t], Battery_capacity))
# 让B_SOC等于B_Bound，确保不超过容量
m.Equations(B_SOC[t] == B_Bound[t] for t in range(row))

同时删除B_SOC=[B_Bound[t] for t in range(row)]这一行，避免覆盖优化变量。

问题2：二进制变量非严格0/1值 #

问题描述 #

二进制变量I_disch在第三组数据中的输出值为1.4295102702e-21，接近0但并非严格整数0，询问是否需要调整求解器参数。

原因分析 #

这是MINLP求解器的常见数值精度问题：求解器在迭代过程中，由于浮点运算误差，可能会输出非常接近整数边界的非严格整数值，本质上代表的是0。

解决方法 #

1. 调整求解器整数容忍度：当前设置的minlp_integer_tol 1e-100过于严格，反而可能导致求解器难以收敛到严格整数。建议调整为1e-6，让求解器将接近0/1的值视为整数：

   m.solver_options = ['minlp_gap_tol 1.0e-2',\
                       'minlp_maximum_iterations 10000',\
                       'minlp_max_iter_with_int_sol 500',\
                       'minlp_branch_method 1',\
                       'minlp_integer_tol 1e-6', \  # 修改此处
                       'minlp_integer_leaves 2']
   

2. 结果后处理：在输出结果时，对二进制变量进行四舍五入，比如：

   print(f'I_disch{t}={round(I_disch[t].value[0], 6)}')
   

修改后的完整代码 #

from gekko import GEKKO
import pandas as pd
import sys

Inp_Data=pd.DataFrame([[3.61,4.55,13],[3.6,5.47,11],[3.47,6.02,11],[3.99,6.12,12],[4.45,3.69,14],[3.78,2.89,15],[3.72,2.02,18]])
col_name=list(Inp_Data)
row=len(Inp_Data)
Demand=Inp_Data.loc[:,col_name[0]].tolist()
RE=Inp_Data.loc[:,col_name[1]].tolist()
Price=Inp_Data.loc[:,col_name[2]].tolist()

#Battery Specifications
Battery_capacity=6 #kW-h
Battery_Initial_capacity=2.4

Big_M=1000000000000
epsilon=sys.float_info.epsilon

m = GEKKO(remote=False)
m.options.SOLVER=1

# 调整整数容忍度
m.solver_options = ['minlp_gap_tol 1.0e-2',\
                    'minlp_maximum_iterations 10000',\
                    'minlp_max_iter_with_int_sol 500',\
                    'minlp_branch_method 1',\
                    'minlp_integer_tol 1e-6', \
                    'minlp_integer_leaves 2']

B_ch=[m.Var(lb=0) for t in range(row)]
B_disch=[m.Var(lb=0) for t in range(row)]
B_Bound=[m.Var(lb=0) for t in range(row)]
B_SOC=[m.Var(lb=0) for t in range(row)]
# 可保留B_SOC_test用于对比，或删除
B_SOC_test=[m.Var(lb=0) for t in range(row)]

I_ch=[m.Var(integer=True,lb=0,ub=1) for t in range(row)]
I_disch=[m.Var(integer=True,lb=0,ub=1) for t in range(row)]

#Objective function
c = [None]*row
for t in range(row):
    c[t] = (Demand[t]*Price[t])

m.Equations(RE[t]-I_ch[t]*B_ch[t]+I_disch[t]*B_disch[t]==Demand[t] for t in range(row))

#Energy Storage constraints - 先计算未受限制的SOC
m.Equations(B_SOC[t]==Battery_Initial_capacity+I_ch[t]*B_ch[t]-I_disch[t]*B_disch[t] for t in range(1))
m.Equations(B_SOC[t]==B_SOC[t-1]+I_ch[t]*B_ch[t]-I_disch[t]*B_disch[t] for t in range(1,row))

#B_SOC_test的约束（可保留用于对比）
m.Equations(B_SOC_test[t]==Battery_Initial_capacity+I_ch[t]*B_ch[t]-I_disch[t]*B_disch[t] for t in range(1))
m.Equations(B_SOC_test[t]==B_SOC_test[t-1]+I_ch[t]*B_ch[t]-I_disch[t]*B_disch[t] for t in range(1,row))

m.Equations((I_ch[t]+I_disch[t])<=1 for t in range(row))
m.Equations((B_ch[t]*B_disch[t])==0 for t in range(row))

# 修正B_SOC的上限约束：B_Bound取B_SOC计算值和容量的最小值，再赋值给B_SOC
for t in range(row):
    m.Equation(B_Bound[t] == m.min2(B_SOC[t], Battery_capacity))
m.Equations(B_SOC[t] == B_Bound[t] for t in range(row))

m.Maximize(m.sum(c))

m.solve(disp=True)

print('Results')
for t in range(row):
    # 对二进制变量四舍五入处理
    i_ch_val = round(I_ch[t].value[0], 6)
    i_disch_val = round(I_disch[t].value[0], 6)
    print(f'I_ch{t}={i_ch_val} I_disch{t}={i_disch_val} B_ch{t}={B_ch[t].value[0]:.4f}  B_disch{t}={B_disch[t].value[0]:.4f}  B_SOC{t}={B_SOC[t].value[0]:.4f} B_Bound{t}={B_Bound[t].value[0]:.4f} B_SOC_test{t}={B_SOC_test[t].value[0]:.4f}')

print('Objective: ' + str(m.options.objfcnval))

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