Как использовать Gekko для оптимизации траектории в дискретном времени

Я пытаюсь использовать Gekko для оптимизации (разряда) зарядки аккумуляторной системы хранения энергии. Цены на электроэнергию в час EP, производство энергии с помощью солнечных панелей PV и спрос на энергию Dem учитываются на всем горизонте (0-24 часа), чтобы минимизировать общие затраты TC. Арбитраж должен происходить по мере того, как аккумулятор (Pbat_ch & Pbat_dis) заряжается (разряжается) в / из сети (Pgrid_in & Pgrid_out) в оптимальные моменты.

В отличие от большинства примеров в Интернете, проблема не формулируется как модель в пространстве состояний, а в основном опирается на экзогенные данные о цене, потреблении и производстве. Ниже описаны 3 конкретных проблемы со ссылкой на Gurobi, весь код, приводящий к следующей ошибке, можно найти внизу этого сообщения.

Exception:  @error: Inequality Definition
 invalid inequalities: z > x < y
 at0x0000016c6b214040>
 STOPPING . . .
  1. Целевая функция - это сумма затрат, связанных с покупкой / продажей электроэнергии в сеть за весь горизонт. Я привык к Gurobi, который позволяет ссылаться на управляемые переменные (PowerGridOut и PowerGridIn = m.MV(...)) в определенных временных интервалах таким образом ([t]).
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()

person Robert_RP    schedule 20.09.2020    source источник


Ответы (1)


Хорошее приложение! Вы можете либо записать все свои дискретные уравнения самостоятельно с помощью 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,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()
person John Hedengren    schedule 21.09.2020
comment
Дорогой Джон, Большое спасибо за ваш быстрый и подробный отзыв, он действительно решает проблемы, которые я не мог бы решить иначе. Если можно, у меня осталась одна проблема: применительно к моим измеренным данным PV и Dem (которые являются плавающими, а не целыми числами) ваша модель приводит к ошибке (проблема может быть недопустимой ... код ошибки - 2). Решение кажется возможным, только если sum(Dem) = sum(PV). Однако это по определению не так: модель должна гарантировать, что спрос удовлетворяется предложением, путем изменения _2 _ / _ 3_, _4 _ / _ 5_ (для каждой пары каждый должен быть равен нулю в каждый момент времени). - person Robert_RP; 22.09.2020
comment
Я понимаю, что последующий вопрос выше несколько отличается от основного, поэтому я разместил его как новый вопрос, используя другой, более интуитивный пример: stackoverflow.com/questions/64033190/ В любом случае еще раз большое спасибо за первый ответ ! - person Robert_RP; 23.09.2020