如何使用 Gekko 进行离散时间的轨迹优化
How to use Gekko for trajectory optimization in discrete time
我正在尝试使用 Gekko 来优化电池储能系统的(放电)充电。每小时电价 EP
、太阳能电池板的能源生产 PV
和能源需求 Dem
在整个时间范围内 (0-24h) 都被考虑在内,以最小化总成本 TC
。当电池在最佳时刻充电 (Pbat_ch
& Pbat_dis
) to/from 电网 (Pgrid_in
& Pgrid_out
) 时,应该进行套利。
与在线的大多数示例相反,该问题并未表述为状态-space 模型,而是主要依赖于价格、消费和生产的外生数据。参考 Gurobi 的 3 个具体问题概述如下,导致以下错误的整个代码可以在 post.
的底部找到
Exception: @error: Inequality Definition
invalid inequalities: z > x < y
at0x0000016c6b214040>
STOPPING . . .
- objective 函数是整个范围内 buying/selling 电力输送到电网所产生的成本总和。我习惯了 Gurobi,它允许以这种方式 ([t]) 在特定时间步长引用操纵变量 (
PowerGridOut
和 PowerGridIn = m.MV(...)
)。
m.Obj(sum(ElectricityPrice[t]*PowerGridOut[t] - ElectricityPrice[t]*PowerGridIn[t]) for t in range(25))
这在 Gekko 中是否也可行,还是应该将此求和重铸为积分?以下代码正确吗?
ElectricityPrice = m.Param([..])
.
.
.
TotalCosts = m.integral(ElectricityPrice*(PowerGridOut - PowerGridIn))
m.Obj(TotalCosts)
m.options.IMODE = 6
m.solve()
- Gurobi 允许这种对电池充电状态变化的约束公式:
m.addConstrs(SoC[t+1] == (SoC[t] - ((1/(DischargeEfficiency*BatteryCapacity)) * (PowerBattery
Discharge[t+1]) * Delta_t - ChargeEfficiency/BatteryCapacity * (PowerBatteryCharge[t+1]) * Delta_t)) for t in range(24))
基于关于类似问题的 Whosebug 问题,我以连续的方式将其重新表述为:
m.Equation(SoC.dt() == SoC - 1/(DischargeEfficiency*BatteryCapacity) * Pbattdis - (ChargeEfficiency/BatteryCapacity) * Pbattch)
- 最终的关键约束应该是权力平衡,其中
Demand[t]
& PV[t]
是外生向量,而其他变量是m.MV()
:
m.Equation(((Demand[t] + Pbat_ch + Pgrid_in) == (PV[t] + Pgrid_out + Pbat_dis)) for t in range(25))
不幸的是,到目前为止所有这些都没有奏效。如果有人能给我一些提示,我将不胜感激。理想情况下,我想用离散的术语来表达 objective 函数和约束。
完整代码
m = GEKKO()
# horizon
m.time = list(range(0,25))
# data vectors
EP = m.Param(list(Eprice))
Dem = m.Param(list(demand))
PV = m.Param(list(production))
# constants
bat_cap = 13.5
ch_eff = 0.94
dis_eff = 0.94
# manipulated variables
Pbat_ch = m.MV(lb=0, ub=4)
Pbat_ch.DCOST = 0
Pbat_ch.STATUS = 1
Pbat_dis = m.MV(lb=0, ub=4)
Pbat_dis.DCOST = 0
Pbat_dis.STATUS = 1
Pgrid_in = m.MV(lb=0, ub=3)
Pgrid_in.DCOST = 0
Pgrid_in.STATUS = 1
Pgrid_out = m.MV(lb=0, ub=3)
Pgrid_out.DCOST = 0
Pgrid_out.STATUS = 1
#State of Charge Battery
SoC = m.Var(value=0.5, lb=0.2, ub=1)
#Battery Balance
m.Equation(SoC.dt() == SoC - 1/(dis_eff*bat_cap) * Pbat_dis - (ch_eff/bat_cap) * Pbat_ch)
#Energy Balance
m.Equation(((Dem[t] + Pbat_ch + Pgrid_in) == (PV[t] + Pbat_dis + Pgrid_out)) for t in range(0,25))
#Objective
TC = m.Var()
m.Equation(TC == sum(EP[t]*(Pgrid_out-Pgrid_in) for t in range(0,25)))
m.Obj(TC)
m.options.IMODE=6
m.options.NODES=3
m.options.SOLVER=3
m.solve()
不错的应用!您可以使用 m.options.IMODE=3
自己写出所有离散方程,或者让 Gekko 为您管理时间维度。当您包含 objective 或约束时,它会将它们应用于您指定的所有时间点。使用 m.options.IMODE=6
,就不需要在 Gekko 中添加集合索引,例如 [t]
。这是一个简化的模型:
from gekko import GEKKO
import numpy as np
m = GEKKO()
# horizon
m.time = np.linspace(0,3,4)
# data vectors
EP = m.Param([0.1,0.05,0.2,0.25])
Dem = m.Param([10,12,9,8])
PV = m.Param([10,11,8,10])
# constants
bat_cap = 13.5
ch_eff = 0.94
dis_eff = 0.94
# manipulated variables
Pbat_ch = m.MV(lb=0, ub=4)
Pbat_ch.DCOST = 0
Pbat_ch.STATUS = 1
Pbat_dis = m.MV(lb=0, ub=4)
Pbat_dis.DCOST = 0
Pbat_dis.STATUS = 1
Pgrid_in = m.MV(lb=0, ub=3)
Pgrid_in.DCOST = 0
Pgrid_in.STATUS = 1
Pgrid_out = m.MV(lb=0, ub=3)
Pgrid_out.DCOST = 0
Pgrid_out.STATUS = 1
#State of Charge Battery
SoC = m.Var(value=0.5, lb=0.2, ub=1)
#Battery Balance
m.Equation(bat_cap * SoC.dt() == -dis_eff*Pbat_dis + ch_eff*Pbat_ch)
#Energy Balance
m.Equation(Dem + Pbat_ch + Pgrid_in == PV + Pbat_dis + Pgrid_out)
#Objective
m.Minimize(EP*Pgrid_in)
# sell power at 90% of purchase (in) price
m.Maximize(0.9*EP*Pgrid_out)
m.options.IMODE=6
m.options.NODES=3
m.options.SOLVER=3
m.solve()
我修改了你的微分方程以从右侧删除 SoC
,否则你将得到指数增长。能量平衡微分方程为Accumulation=In-Out
。下面是一些额外的代码来可视化解决方案。
import matplotlib.pyplot as plt
plt.subplot(3,1,1)
plt.plot(m.time,SoC.value,'b--',label='State of Charge')
plt.ylabel('SoC')
plt.legend()
plt.subplot(3,1,2)
plt.plot(m.time,Dem.value,'r--',label='Demand')
plt.plot(m.time,PV.value,'k:',label='PV Production')
plt.legend()
plt.subplot(3,1,3)
plt.plot(m.time,Pbat_ch.value,'g--',label='Battery Charge')
plt.plot(m.time,Pbat_dis.value,'r:',label='Battery Discharge')
plt.plot(m.time,Pgrid_in.value,'k--',label='Grid Power In')
plt.plot(m.time,Pgrid_in.value,':',color='orange',label='Grid Power Out')
plt.ylabel('Power')
plt.legend()
plt.xlabel('Time')
plt.show()
我正在尝试使用 Gekko 来优化电池储能系统的(放电)充电。每小时电价 EP
、太阳能电池板的能源生产 PV
和能源需求 Dem
在整个时间范围内 (0-24h) 都被考虑在内,以最小化总成本 TC
。当电池在最佳时刻充电 (Pbat_ch
& Pbat_dis
) to/from 电网 (Pgrid_in
& Pgrid_out
) 时,应该进行套利。
与在线的大多数示例相反,该问题并未表述为状态-space 模型,而是主要依赖于价格、消费和生产的外生数据。参考 Gurobi 的 3 个具体问题概述如下,导致以下错误的整个代码可以在 post.
的底部找到Exception: @error: Inequality Definition
invalid inequalities: z > x < y
at0x0000016c6b214040>
STOPPING . . .
- objective 函数是整个范围内 buying/selling 电力输送到电网所产生的成本总和。我习惯了 Gurobi,它允许以这种方式 ([t]) 在特定时间步长引用操纵变量 (
PowerGridOut
和PowerGridIn = m.MV(...)
)。
m.Obj(sum(ElectricityPrice[t]*PowerGridOut[t] - ElectricityPrice[t]*PowerGridIn[t]) for t in range(25))
这在 Gekko 中是否也可行,还是应该将此求和重铸为积分?以下代码正确吗?
ElectricityPrice = m.Param([..])
.
.
.
TotalCosts = m.integral(ElectricityPrice*(PowerGridOut - PowerGridIn))
m.Obj(TotalCosts)
m.options.IMODE = 6
m.solve()
- Gurobi 允许这种对电池充电状态变化的约束公式:
m.addConstrs(SoC[t+1] == (SoC[t] - ((1/(DischargeEfficiency*BatteryCapacity)) * (PowerBattery
Discharge[t+1]) * Delta_t - ChargeEfficiency/BatteryCapacity * (PowerBatteryCharge[t+1]) * Delta_t)) for t in range(24))
基于关于类似问题的 Whosebug 问题,我以连续的方式将其重新表述为:
m.Equation(SoC.dt() == SoC - 1/(DischargeEfficiency*BatteryCapacity) * Pbattdis - (ChargeEfficiency/BatteryCapacity) * Pbattch)
- 最终的关键约束应该是权力平衡,其中
Demand[t]
&PV[t]
是外生向量,而其他变量是m.MV()
:
m.Equation(((Demand[t] + Pbat_ch + Pgrid_in) == (PV[t] + Pgrid_out + Pbat_dis)) for t in range(25))
不幸的是,到目前为止所有这些都没有奏效。如果有人能给我一些提示,我将不胜感激。理想情况下,我想用离散的术语来表达 objective 函数和约束。
完整代码
m = GEKKO()
# horizon
m.time = list(range(0,25))
# data vectors
EP = m.Param(list(Eprice))
Dem = m.Param(list(demand))
PV = m.Param(list(production))
# constants
bat_cap = 13.5
ch_eff = 0.94
dis_eff = 0.94
# manipulated variables
Pbat_ch = m.MV(lb=0, ub=4)
Pbat_ch.DCOST = 0
Pbat_ch.STATUS = 1
Pbat_dis = m.MV(lb=0, ub=4)
Pbat_dis.DCOST = 0
Pbat_dis.STATUS = 1
Pgrid_in = m.MV(lb=0, ub=3)
Pgrid_in.DCOST = 0
Pgrid_in.STATUS = 1
Pgrid_out = m.MV(lb=0, ub=3)
Pgrid_out.DCOST = 0
Pgrid_out.STATUS = 1
#State of Charge Battery
SoC = m.Var(value=0.5, lb=0.2, ub=1)
#Battery Balance
m.Equation(SoC.dt() == SoC - 1/(dis_eff*bat_cap) * Pbat_dis - (ch_eff/bat_cap) * Pbat_ch)
#Energy Balance
m.Equation(((Dem[t] + Pbat_ch + Pgrid_in) == (PV[t] + Pbat_dis + Pgrid_out)) for t in range(0,25))
#Objective
TC = m.Var()
m.Equation(TC == sum(EP[t]*(Pgrid_out-Pgrid_in) for t in range(0,25)))
m.Obj(TC)
m.options.IMODE=6
m.options.NODES=3
m.options.SOLVER=3
m.solve()
不错的应用!您可以使用 m.options.IMODE=3
自己写出所有离散方程,或者让 Gekko 为您管理时间维度。当您包含 objective 或约束时,它会将它们应用于您指定的所有时间点。使用 m.options.IMODE=6
,就不需要在 Gekko 中添加集合索引,例如 [t]
。这是一个简化的模型:
from gekko import GEKKO
import numpy as np
m = GEKKO()
# horizon
m.time = np.linspace(0,3,4)
# data vectors
EP = m.Param([0.1,0.05,0.2,0.25])
Dem = m.Param([10,12,9,8])
PV = m.Param([10,11,8,10])
# constants
bat_cap = 13.5
ch_eff = 0.94
dis_eff = 0.94
# manipulated variables
Pbat_ch = m.MV(lb=0, ub=4)
Pbat_ch.DCOST = 0
Pbat_ch.STATUS = 1
Pbat_dis = m.MV(lb=0, ub=4)
Pbat_dis.DCOST = 0
Pbat_dis.STATUS = 1
Pgrid_in = m.MV(lb=0, ub=3)
Pgrid_in.DCOST = 0
Pgrid_in.STATUS = 1
Pgrid_out = m.MV(lb=0, ub=3)
Pgrid_out.DCOST = 0
Pgrid_out.STATUS = 1
#State of Charge Battery
SoC = m.Var(value=0.5, lb=0.2, ub=1)
#Battery Balance
m.Equation(bat_cap * SoC.dt() == -dis_eff*Pbat_dis + ch_eff*Pbat_ch)
#Energy Balance
m.Equation(Dem + Pbat_ch + Pgrid_in == PV + Pbat_dis + Pgrid_out)
#Objective
m.Minimize(EP*Pgrid_in)
# sell power at 90% of purchase (in) price
m.Maximize(0.9*EP*Pgrid_out)
m.options.IMODE=6
m.options.NODES=3
m.options.SOLVER=3
m.solve()
我修改了你的微分方程以从右侧删除 SoC
,否则你将得到指数增长。能量平衡微分方程为Accumulation=In-Out
。下面是一些额外的代码来可视化解决方案。
import matplotlib.pyplot as plt
plt.subplot(3,1,1)
plt.plot(m.time,SoC.value,'b--',label='State of Charge')
plt.ylabel('SoC')
plt.legend()
plt.subplot(3,1,2)
plt.plot(m.time,Dem.value,'r--',label='Demand')
plt.plot(m.time,PV.value,'k:',label='PV Production')
plt.legend()
plt.subplot(3,1,3)
plt.plot(m.time,Pbat_ch.value,'g--',label='Battery Charge')
plt.plot(m.time,Pbat_dis.value,'r:',label='Battery Discharge')
plt.plot(m.time,Pgrid_in.value,'k--',label='Grid Power In')
plt.plot(m.time,Pgrid_in.value,':',color='orange',label='Grid Power Out')
plt.ylabel('Power')
plt.legend()
plt.xlabel('Time')
plt.show()