GEKKO中MILP优化问题的启停行为实现技术问询

问题场景 #

• 2台带上下运行边界的发电机
• 其中一台发电机（gen1）存在30%的功率损耗
• 总发电量必须满足外部用电需求
• 优化目标为最小化能量损耗

原始实现代码 #

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

timesim = 24*1      # hours
timesteps = 1   # steps/hour
n = np.int64(timesim*timesteps + 1) # this is the vector length in GEKKO

Edh_demand =np.ones(timesim+1)*80; 
Edh_demand[int(7*n/24):int(10*n/24)] = 200
Edh_demand[int(3*n/4)-1:int(3*n/4)+2] = 300


m = GEKKO(remote=False)             # initialize gekko, körs lokalt
m.time = np.linspace(0, 24, 25)      # Gekko-tid

Egen1 = m.Var(value = 30, lb = 30, ub = 200)
Egen2 = m.Var(value = 30, lb = 30, ub = 200)
Eloss = m.Var(value = 30*.3)
Eprod = m.Var(value = 30)
Edemand = m.Param(value = Edh_demand)

cost = m.Param(value = 10)

m.Equation(Eloss == Egen1*.3)
m.Equation(Eprod == Egen2 + Egen1 - Eloss)

m.Minimize(Eloss)

m.Equation(Eprod >= Edemand)

m.options.SOLVER = 3
m.options.IMODE = 6
# m.options.COLDSTART=2
m.options.COLDSTART=0

m.solve(disp=True)

plt.figure(5)
plt.subplot(2, 1, 1)
plt.plot(m.time, Egen1, label='gen 1')
plt.plot(m.time, Egen2, label='gen 2')
plt.legend()

plt.subplot(2, 1, 2)
plt.plot(m.time, Edemand, '--r')
plt.plot(m.time, Eprod)

plt.show()

技术问询 #

1. 当前实现中gen1会运行至下限，如何约束gen1仅能处于下限或0状态，不可在0与下限区间运行？
2. 如何对gen1的启停操作设置惩罚机制？
3. 如何约束gen1每次停机后，需间隔2个采样间隔才可再次启动？
4. 是否存在算法效率更高、可读性更好的该问题实现方案？

初步思路 #

可创建整数变量gen1_onoff，将其与gen1的所有引用项相乘，同时添加约束m.Equation(gen1 * gen1_onoff <= gen1_lb)；可通过在目标函数中惩罚gen1_onoff的变化来实现启停惩罚，示例代码如下：

dt_gen1_onoff = Var(value = 0) 
m.Equation(dt_gen1_onoff == gen1_onoff.dt())
m.Minimize(dt_gen1_onoff)

但对于第3点约束的实现方式完全不清楚，猜测可使用m.if3()函数。

1. 约束gen1仅处于0或下限状态 #

引入二进制整数变量gen1_onoff（取值0或1），通过直接绑定变量值实现二选一状态：

• 当gen1_onoff=1时，Egen1强制等于下限值；当gen1_onoff=0时，Egen1为0
• 替换原Egen1的上下限设置，改为如下约束：

gen1_onoff = m.Var(value=1, lb=0, ub=1, integer=True)
Egen1 = m.Var(value=30, lb=0, ub=200)
m.Equation(Egen1 == 30 * gen1_onoff)  # 强制Egen1只能是0或30

2. 启停操作惩罚机制 #

通过捕捉gen1_onoff的状态变化（0→1或1→0），将变化量乘以惩罚系数加入目标函数：

• 用辅助变量计算状态变化的绝对值，避免正负抵消
• 实现代码：

# 计算状态变化的绝对值
delta_onoff = m.Var(lb=0)
m.Equation(delta_onoff >= gen1_onoff - m.previous(gen1_onoff))
m.Equation(delta_onoff >= m.previous(gen1_onoff) - gen1_onoff)
# 添加启停惩罚，系数可按需调整（示例用10）
m.Minimize(10 * delta_onoff)

3. 停机后间隔2个采样间隔才能启动 #

通过逻辑约束限制启动时机：确保当gen1_onoff从0切换到1时，前两个采样间隔必须处于停机状态，直接用以下约束即可实现：

# 约束：当前启动的前提是前两个时刻至少有一个处于运行状态，否则必须停机满2步才能启动
m.Equation(gen1_onoff <= m.previous(gen1_onoff) + m.previous(m.previous(gen1_onoff)))

也可以通过跟踪停机时长的方式实现，适合更复杂的间隔需求：

off_duration = m.Var(value=0, lb=0)
# 停机时累计时长，启动时重置为0
m.Equation(off_duration == m.previous(off_duration) * (1 - gen1_onoff) + (1 - gen1_onoff))
# 启动时必须满足停机时长≥2
m.Equation((1 - m.previous(gen1_onoff)) * gen1_onoff * (off_duration - 2) <= 0)

4. 高效可读的实现方案 #

推荐采用模块化结构+混合整数优化器，提升求解效率和代码可读性：

1. 用二进制变量明确标记启停状态，避免模糊的区间约束
2. 切换到SOLVER=1（APOPT），该求解器更擅长处理混合整数动态优化问题
3. 拆分约束模块，将发电机状态、损耗、需求、启停逻辑分别封装
4. 启用COLDSTART=2帮助求解器快速定位可行解

优化后的核心代码框架：

m = GEKKO(remote=False)
m.time = np.linspace(0,24,25)

# 定义参数
gen1_lb = 30
gen1_ub = 200
gen2_lb = 30
gen2_ub = 200
loss_rate = 0.3
startup_penalty = 10

# 变量定义
gen1_on = m.Var(value=1, lb=0, ub=1, integer=True)
Egen1 = m.Var(value=gen1_lb, lb=0, ub=gen1_ub)
Egen2 = m.Var(value=gen2_lb, lb=gen2_lb, ub=gen2_ub)
Eloss = m.Var()
Eprod = m.Var()
Edemand = m.Param(value=Edh_demand)

# 核心约束
m.Equation(Egen1 == gen1_lb * gen1_on)
m.Equation(Eloss == Egen1 * loss_rate)
m.Equation(Eprod == Egen1 + Egen2 - Eloss)
m.Equation(Eprod >= Edemand)

# 启停惩罚
delta_gen1 = m.Var(lb=0)
m.Equation(delta_gen1 >= gen1_on - m.previous(gen1_on))
m.Equation(delta_gen1 >= m.previous(gen1_on) - gen1_on)
m.Minimize(startup_penalty * delta_gen1)

# 停机后启动间隔约束
m.Equation(gen1_on <= m.previous(gen1_on) + m.previous(m.previous(gen1_on)))

# 求解设置
m.options.SOLVER = 1  # APOPT适配混合整数问题
m.options.IMODE = 6
m.options.COLDSTART = 2
m.solve(disp=True)

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