title: Having Optimizer switch power storage based on Capacities and net load
body: <p>I am trying to have an optimizer use two power storages for a cyclic power generation and consumption cycle. The goal is to have it turn off the primary energy storage(in this case electric battery storage) when it has reached full capacity, and then to discharge first until it is empty. The secondary storage is to charge after the primary and discharge after the primary storage. I would like the optimizer to solve as it goes based on the system. I have tried using a series of switches, but it isn't quite working. I know that using if statements are tricky for gradient based solvers so if there is any help that would be great thanks!</p>
  
  <pre><code>Capacity = 76.2
  EStored = m.SV(value=0,lb=0,ub=Capacity)
  batteryeff = .95
  batteff = m.if3((Enuc - Cons),1/batteryeff,batteryeff)
  
  #Energy Balance 
  Cost = m.Var()
  eneed = m.sign2(Enuc-Cons)  #gives sign of energy need for energy storage or removal from storage
  eswitch = m.if3(eneed,0, 1) #Turns eneed into a binary switch
  switch = m.if3(eswitch*(Capacity-EStored)+(1-eswitch)*(EStored),0,1) #supposed to charge battery until at capacity and then discharge until EStored is 0. Then use thermal energy second
  
  
  m.Equation(EStored.dt() == (switch)*batteff*(Enuc - Cons))  #Energy balance for Battery
  m.Equation(T.dt() == (1-switch)*thermeff*(Enuc - Cons)/(mass*Cp)) #Energy balance for Thermal Storage
  m.Equation(Cost == Enuc*1000 )
  m.Obj(Cost)
  
  m.options.IMODE = 5
  m.options.SOLVER = 3
  m.solve()
  </code></pre>
  
  <p>This is the section of the code that I am working with however if more details are necissary here is the simplified version of total code.</p>
  
  <pre><code>from gekko import GEKKO
  import numpy as np
  import pandas as pd
  import matplotlib.pyplot as plt
  from scipy.integrate import quad
  
  #Set up basic power consumption data
  n = 24
  t = np.linspace(0,n,n)
  def load(t):
      return  -10*np.sin(2*np.pi*t/24)+40
  Load = load(t)
  Gen = np.ones(n)*40
  def need(t):
      return 10*np.sin(2*np.pi*t/24)+10
  
  #Set up Model
  m = GEKKO()
  m.time = t
  
  Cons = m.Param(value=Load)
  Enuc = m.FV(value=45, lb=0) #nuclear power
  Enuc.STATUS = 1
  
  #Thermal Energy Storage
  T = m.SV(value=300,ub=500,lb=300)
  mass = m.FV(value=.0746,lb=0)
  mass.STATUS=0
  Cp = m.Param(value=5)
  thermaleff = .8   #80%efficient
  thermeff = m.if3((Enuc - Cons)/(mass*Cp),1/thermaleff,thermaleff)
  
  #Battery Electrical storage
  Capacity = 76.2
  EStored = m.SV(value=0,lb=0,ub=Capacity)
  batteryeff = .95
  batteff = m.if3((Enuc - Cons),1/batteryeff,batteryeff)
  
  #Energy Balance 
  Cost = m.Var()
  eneed = m.sign2(Enuc-Cons)  #gives sign of energy need for energy storage or removal from storage
  eswitch = m.if3(eneed,0, 1) #Turns eneed into a binary switch
  switch = m.if3(eswitch*(Capacity-EStored)+(1-eswitch)*(EStored),0,1) #supposed to charge battery until at capacity and then discharge until EStored is 0. Then use thermal energy second
  
  
  m.Equation(EStored.dt() == (switch)*batteff*(Enuc - Cons))  #Energy balance for Battery
  m.Equation(T.dt() == (1-switch)*thermeff*(Enuc - Cons)/(mass*Cp)) #Energy balance for Thermal Storage
  m.Equation(Cost == Enuc*1000 )
  m.Obj(Cost)
  
  m.options.IMODE = 5
  m.options.SOLVER = 3
  m.solve()
  
  #plot
  plt.subplot(3,1,1)
  plt.plot(t,Load)
  plt.plot(t,Enuc.value, label=f'Enuc = {Enuc.value[-1]}')
  plt.ylabel("Energy")
  plt.legend()
  
  plt.subplot(3,1,2)
  plt.plot(t,EStored.value, label=f'Capacity = {np.max(EStored.value):.03}')
  plt.title("Battery Storage")
  plt.ylabel("Energy")
  plt.legend()
  
  plt.subplot(3,1,3)
  plt.plot(t,T.value,label=f'mass = {mass.value[-1]:.03}')
  plt.title("Thermal Storage")
  plt.ylabel("Temperature(K)")
  plt.legend()
  plt.show()```
  </code></pre>