4

I am trying to run an electricity arbitrage model in Gekko. I have an electricity price array for every hour of a year (8760 total hours), a battery of energy size E, and for every hour I want to decide whether to charge the battery or discharge it, based on minimizing the electricity cost, and keep track of the energy in the battery constrained to hold no less than 0 energy and no more than E.

I've tried many times, most recently getting the error that the optimize equation exceeds the limits

from gekko import Gekko

m = Gekko()

#variables
E_battery = m.Var(lb=0, ub=366.2, value=0) #energy in battery at time t, battery size 366 MWh
Pc = m.Var(lb=0, ub=50) #charge power, 50 MW max
Pd = m.Var(lb=0, ub=36.6)  #discharge power, max 36 MW
E_price = m.Param(electricity_price[:,1])
m.time = np.linspace(0,8759, 8760)

m.Equation(E_battery.dt() == (1-delta)*E_battery + roundtrip_eff*(Pc - Pd))

m.Obj(sum(E_price[i]*Pc for i in range(8760)))
m.options.IMODE = 7
m.solve()
David
  • 51
  • 3

1 Answers1

1

If you are optimizing over the entire time horizon then you'll want to switch to IMODE=6. You may be having a problem with:

m.Obj(sum(E_price[i]*Pc for i in range(8760)))

because it creates an objective expression that is very long. For dynamic optimization problems with IMODE=6 you can use the following instead:

m.Obj(E_price*Pc)

Gekko automatically does the summation of all time points in that expression. The script isn't a complete example so it is hard to verify what isn't working. If you can post a minimal and complete example that shows the issue then it will be easier to provide helpful feedback. There are additional examples in the machine learning and dynamic optimization course, especially with economic optimization or some of the other benchmark problems.

Summation

You can sum a variable such as E_battery with something like the following:

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt  
m = GEKKO(remote=False)
m.time = [0,1,2,3,4,5]
x = m.Param(value=m.time)
E_battery = m.Var()
m.Equation(E_battery.dt()==x)
m.options.IMODE = 4
m.solve(disp=False)
print('x = ' + str(x.value))
print('E_battery = ' + str(E_battery.value))

This produces the solution:

x = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
E_battery = [0.0, 1.0, 3.0, 6.0, 10.0, 15.0]
John Hedengren
  • 12,068
  • 1
  • 21
  • 25
  • 1
    Thank you John! The solver ran with your suggestions and was able to find an exact Hessian! – David May 13 '20 at 15:03
  • 1
    quick question. When I write the "m.Equation(E_battery == E_battery + Pc - Pd", is it taking the previous instance of "E_Battery", adding "Pc", subtracting "Pd" for that new time instance, and using that as the new time instance for "E_battery"? – David May 14 '20 at 20:54
  • All equations are solved simultaneously so that equation is equivalent to `m.Equation(E_battery-E_battery==Pc-Pd)` and that is also equivalent to `m.Equation(Pc==Pd)`. If you need to integrate `Pc-Pd` then I recommend the `m.Integral()` function with `E_battery=m.Integral(Pc-Pd)`. – John Hedengren May 14 '20 at 21:53
  • 1
    Thanks so much for responding. I don't need to integrate. I need the value of E_battery to iterate through the time steps such that E_battery[i] is equal to E_battery[i-1]+Pc[i]-Pd[i]. No wonder my optimization gave such a weird answer! – David May 15 '20 at 01:36
  • The integral will give you the summation but you may need to divide by the time step if you don't want to factor in the time. Another way to do this is to use `E_battery.dt()==Pc[i]-Pd[i]`. – John Hedengren May 15 '20 at 22:15