Mixed-Integer Model Predictive Control using Gekko

Question

Problem statement:

Trying to utilize GEKKO to optimally select the number of appliances within the available capacity while minimizing the energy cost.

Problem model:

I cannot post an image but here is the link to see the formula and also you can copy and paste below to find it visually better:

\min \sum _{t=1}^{24} n x_{i} P_{i}

subject to:

\sum _{i=1}^{24} C_{i} - x_{i} = 0 
C_{min} < C_{i}< C_{max}
x_{min} < x_{i}< x_{max}
n_{min} < n_{i}< n_{max}
\sum_{t=1}^{T} \tau_{i}^{t} = I_{i}

where, x and n are the state variables representing the assigned power and number of appliances in a station. C is the capacity, is the operating state of appliance (the ON/OFF condition) at time slot t, and I_{i} is the permissible interval for each appliances (i.e. 5 hrs)

Work done so far:

Thanks to post here I put the following together:

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

m = GEKKO(remote=False)
m.time = np.linspace(0,23,24)

#initialize variables
cap_low = [55.,55.,55.,55.,55.,55.,55.,68.,68.,68.,68.,\
           55.,55.,68.,68.,68.,68.,55.,55.,55.,55.,55.,55.,55.]
cap_upper = [75.,75.,75.,75.,75.,75.,75.,70.,70.,70.,70.,\
             75.,75.,70.,70.,70.,70.,75.,75.,75.,75.,75.,75.,75.]
TOU = [0.05,0.05,0.05,0.05,0.05,0.05,0.05,50.,50.,50.,50.,50.,\
       50.,50.,50.,50.,50.,50.,50.,50.,50.,0.05,0.05,0.05]

k = 1.1 #Safety factor
CL = m.Param(value=cap_low)
CH = m.Param(value=cap_upper)
TOU = m.Param(value=TOU)

x = m.MV(lb=6, ub=32)
x.STATUS = 1  # allow optimizer to change

n = m.MV(lb=1, ub=10, integer=True)
n.STATUS = 1  # allow optimizer to change

# Controlled Variable
C = m.SV(value=55)
m.Equations([C>=CL, C<=CH])

m.Equation(k*C - n* x == 0)
m.Minimize(TOU*n*x)

m.options.IMODE = 6
m.solve(disp=True,debug=True)

plt.subplot(3,1,1)
plt.plot(m.time,cap_low,'k--')
plt.plot(m.time,cap_upper,'k--')
plt.plot(m.time,C.value,'r-')
plt.ylabel('Demand')
plt.subplot(3,1,2)
plt.step(m.time,x.value,'b:')
plt.ylabel('Load')
plt.xlabel('Time (hr)')
plt.subplot(3,1,3)
plt.step(m.time,n.value,'r:')
plt.ylabel('No. of Appliances')
plt.xlabel('Time (hr)')
plt.show()

Outcome

The above code results in a promising result but there is a lack of last constraint \sum_{t=1}^{T} \tau_{i}^{t} = I_{i} for limiting the operation time of each appliances.

Question:

This problem seems like a combination of Model Predictive Control (tracking the available load capacity) and a Mixed Integer Linear Programing (multiple appliances with multiple linear operational equations). So the questions are:

1. How can I define the permissible interval for each appliances in above code to limit their operation in a day.

2. What is the best way to integrate MPC and MILP in Gekko

Answer 1

Answer to question 2: Switch to the APOPT solver to solve MILP or MINLP problems.

Answer to question 1: Separate n into the individual appliance state with:

# decide when to turn on each appliance
na = 20 # number of appliances
a = m.Array(m.MV,na,lb=0, ub=1, integer=True)
for ai in a:
    ai.STATUS = 1
    ai.DCOST = 0

# number of appliances running
n = m.Var(lb=1, ub=5)
m.Equation(n==m.sum(a))

The bottom subplot shows when the individual appliances are on or off. Limiting the appliance to 5 hours per day is causing problems with convergence:

ax = [m.Intermediate(m.integral(ai)) for ai in a]
for ai in ax:
    # maximum of 5 hours each day
    m.Minimize(1e-4*ax**2)
    m.Equation(ax<=5)

Perhaps you could start with this code and try other strategies to improve convergence such as different APOPT solver options or form of the constraint. Creating a soft constraint (penalize with a financial incentive instead of a hard limit) can also help to improve convergence.

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

m = GEKKO(remote=True)
m.time = np.linspace(0,23,24)

#initialize variables
cap_low = [55.,55.,55.,55.,55.,55.,55.,68.,68.,68.,68.,\
           55.,55.,68.,68.,68.,68.,55.,55.,55.,55.,55.,55.,55.]
cap_upper = [75.,75.,75.,75.,75.,75.,75.,70.,70.,70.,70.,\
             75.,75.,70.,70.,70.,70.,75.,75.,75.,75.,75.,75.,75.]
TOU = [0.05,0.05,0.05,0.05,0.05,0.05,0.05,50.,50.,50.,50.,50.,\
       50.,50.,50.,50.,50.,50.,50.,50.,50.,0.05,0.05,0.05]

k = 1.1 #Safety factor
CL = m.Param(value=cap_low)
CH = m.Param(value=cap_upper)
TOU = m.Param(value=TOU)

x = m.MV(lb=6, ub=32)
x.STATUS = 1  # allow optimizer to change

# decide when to turn on each appliance
na = 20 # number of appliances
a = m.Array(m.MV,na,lb=0, ub=1, integer=True)
for ai in a:
    ai.STATUS = 1
    ai.DCOST = 0

ax = [m.Intermediate(m.integral(ai)) for ai in a]
for ai in ax:
    # maximum of 5 hours each day
for ai in ax:
    # maximum of 5 hours each day
    m.Minimize(1e-4*ai**2)
    m.Equation(ai<=5)

# number of appliances running
n = m.Var(lb=1, ub=5)
m.Equation(n==m.sum(a))

# Controlled Variable
C = m.SV(value=55)
m.Equations([C>=CL, C<=CH])

m.Equation(k*C - n* x == 0)
m.Minimize(TOU*n*x)

m.options.IMODE = 6
m.options.NODES = 2
m.options.SOLVER = 1

m.solver_options = ['minlp_gap_tol 1.0e-2',\
                    'minlp_maximum_iterations 10000',\
                    'minlp_max_iter_with_int_sol 500',\
                    'minlp_branch_method 3']

m.solve(disp=True)

plt.subplot(4,1,1)
plt.plot(m.time,cap_low,'k--')
plt.plot(m.time,cap_upper,'k--')
plt.plot(m.time,C.value,'r-')
plt.ylabel('Demand')
plt.subplot(4,1,2)
plt.step(m.time,x.value,'b:')
plt.ylabel('Load')
plt.xlabel('Time (hr)')
plt.subplot(4,1,3)
plt.step(m.time,n.value,'r:')
plt.ylabel('No. of Appliances')
plt.subplot(4,1,4)
for ai in a:
    plt.plot(m.time,ai.value)
plt.xlabel('Time (hr)')
plt.show()