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

# uncertain delay instances
n = 6

m = GEKKO()
m.time = np.linspace(0,20,41)

# manipulated variable
p = m.MV(value=0, lb=0, ub=100)
p.STATUS = 1 
p.DCOST = 0.1  
p.DMAX = 20

# controlled variable
v = m.Array(m.CV,n)
d = m.Array(m.Var,n)
for i in range(n):
    delay = np.random.randint(1,11)
    v[i].STATUS = 1
    v[i].SPHI = 40
    v[i].SPLO = 38
    v[i].TAU = 5
    v[i].TR_INIT = 1
    m.delay(p,d[i],delay)
    m.Equation(2*v[i].dt()==-v[i]+d[i])

# solve optimal control problem
m.options.IMODE = 6
m.options.CV_TYPE = 1
m.solve()

# get additional solution information
import json
with open(m.path+'//results.json') as f:
    results = json.load(f)

# plot results
plt.figure()
plt.subplot(2,1,1)
plt.plot(m.time,p.value,'b-',lw=2)
plt.ylabel('MV')
plt.subplot(2,1,2)
plt.plot(m.time,results['v1.tr_hi'],'k--')
plt.plot(m.time,results['v1.tr_lo'],'k--')
for i in range(n):
    plt.plot(m.time,v[i].value,':',lw=2)
plt.ylabel('CV')
plt.xlabel('Time')
plt.show()