如何在离散时间内使用Gekko进行轨迹优化

问题描述

我正在尝试使用Gekko来优化（充电）电池储能系统的充电。每小时（0-24h）均应考虑每小时的电价EP，太阳能电池板的发电量PV和能源需求Dem，以降低总成本TC。当电池在最佳时刻从电网（Pbat_ch和Pbat_dis）放电（Pgrid_in和Pgrid_out）时，应进行套利。>

与大多数在线示例相反，该问题并非公式化为状态空间模型，而是主要依赖于价格，消费和生产的外部数据。下面概述了有关Gurobi的3个具体问题，导致以下错误的整个代码可以在本文的底部找到。

Exception:  @error: Inequality Definition
 invalid inequalities: z > x < y
 at0x0000016c6b214040>
 STOPPING . . .

1. 目标函数是在整个范围内因向电网购买/出售电力而产生的成本总和。我习惯了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()

1. 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))

基于关于类似问题的stackoverflow问题，我以连续方式将其重新表述为：

m.Equation(SoC.dt() == SoC - 1/(DischargeEfficiency*BatteryCapacity) * Pbattdis - (ChargeEfficiency/BatteryCapacity) * Pbattch)

1. 最后的键约束应该是力量平衡，其中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))

不幸的是，到目前为止，所有这些都还没有奏效。如果有人可以给我一些提示，我将不胜感激。理想情况下，我想用离散术语表述目标函数和约束条件。

整个代码

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为您管理时间维度。当您包含目标或约束时，它将它们应用于您指定的所有时间点。使用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,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,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,':',color='orange',label='Grid Power Out')
plt.ylabel('Power')
plt.legend()
plt.xlabel('Time')
plt.show()