Search code examples
pythongekko

Handling on/off behavior in MILP optimization problems in GEKKO

I'm trying to understand how to implement ON/OFF behavior in optimization problems to be solved in GEKKO effectively.

Consider the following scenario:

  • 2 power generators, with upper and lower bounds
  • one generator has a power loss of 30%
  • the total generated power must meet an external power demand
  • the goal is to minimize the energy loss

The above problem is implemented in Python:

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

timesim = 24*1      # hours
timesteps = 1   # steps/hour
n = np.int64(timesim*timesteps + 1) # this is the vector length in GEKKO

Edh_demand =np.ones(timesim+1)*80; 
Edh_demand[int(7*n/24):int(10*n/24)] = 200
Edh_demand[int(3*n/4)-1:int(3*n/4)+2] = 300


m = GEKKO(remote=False)             # initialize gekko, körs lokalt
m.time = np.linspace(0, 24, 25)      # Gekko-tid

Egen1 = m.Var(value = 30, lb = 30, ub = 200)
Egen2 = m.Var(value = 30, lb = 30, ub = 200)
Eloss = m.Var(value = 30*.3)
Eprod = m.Var(value = 30)
Edemand = m.Param(value = Edh_demand)

cost = m.Param(value = 10)

m.Equation(Eloss == Egen1*.3)
m.Equation(Eprod == Egen2 + Egen1 - Eloss)

m.Minimize(Eloss)

m.Equation(Eprod >= Edemand)

m.options.SOLVER = 3
m.options.IMODE = 6
# m.options.COLDSTART=2
m.options.COLDSTART=0

m.solve(disp=True)

plt.figure(5)
plt.subplot(2, 1, 1)
plt.plot(m.time, Egen1, label='gen 1')
plt.plot(m.time, Egen2, label='gen 2')
plt.legend()

plt.subplot(2, 1, 2)
plt.plot(m.time, Edemand, '--r')
plt.plot(m.time, Eprod)

plt.show()

I have the following questions/problems:

  1. As it is implemented now, gen1 goes to its lower bound as expected. I want to constraint gen1 so that it can operate either at the lower bound or at 0, but not between 0 an gen1_lb.
  2. I want to penalize turning off and on gen1.
  3. I want to constraint gen1 so that every time it is turned off, it takes 2 sampling intervals before it can be turned on again.
  4. Are there better ways (mainly for algoritm efficiency but also readibility) to implement such problem?

I believe one way to attack this problem is by creating an Integer variable (gen1_onoff) and multiply all occurances of gen1 by gen1_onoff. One can also add the following constraint

m.Equation(gen1 * gen1_onoff <= gen1_lb)

One could penalize variations of gen1_onoff in the objective function

dt_gen1_onoff = Var(value = 0) 
m.Equation(dt_gen1_onoff == gen1_onoff.dt())
m.Minimize(dt_gen1_onoff)

I have no idea how to implement the constraint for problem 3, but suspect that m.if3() should be useful.

Solution

  • As you correctly noted, when gen1 can only operate above the lower bound (gen1_lb) or at 0, use a binary variable gen1_onoff to switch the operation mode:

    m.options.SOLVER = 1 # APOPT solver for Mixed Integer solutions
gen1_onoff = m.Var(value=1, integer=True, lb=0, ub=1)  # Binary variable
Egen1 = m.Var(value = 30, lb = 30, ub = 200)
Egen1s = m.Var(value = 30, lb = 0, ub = 200)
m.Equation(Egen1s == gen1_onoff * Egen1)  # Egen1s is either at 0 or 30-200

    To create a penalty for changes in gen1_onoff, tracks the change in state and include it in the objective function with a squared error or m.abs3() as an absolute value.

    dt_gen1_onoff = m.Var(value=1, integer=True, lb=0, ub=1) 
m.Equation(dt_gen1_onoff == gen1_onoff.dt())  # Change in on/off state
m.Minimize(cost * m.abs3(dt_gen1_onoff))  # Add to objective

    The m.abs3() ensures that the penalty applies regardless of the sign of the change. However, this introduces another binary variable and greatly slows down the solution.

    A computationally more efficient way is to declare gen1_onoff as an MV type and add a DCOST value to penalize movement.

    m.options.CV_TYPE = 1 # l1-norm objective
gen1_onoff = m.MV(value=1, integer=True, lb=0, ub=1)
gen1_onoff.STATUS = 1
gen1_onoff.DCOST = 1e-3 # penalty on movement

    To impose a delay so that gen1 must remain off for two time steps before being turned back on, use the GEKKO MV_STEP_HOR option. This requires that the MV is held for 2 cycles before changing.

    gen1_onoff.MV_STEP_HOR = 2

    Using several m.if3() functions is also possible but would greatly increase the computational time. To improve efficiency, minimize the use of integer variables as much as possible since they make the problem more complex.