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

t = np.linspace(0,1,101)
m = GEKKO(remote=False); m.time=t

d1 = m.Param(np.cos(2*np.pi*t)+3)
d2 = m.Param(1.5*np.sin(2*np.pi*t)+7)
d3 = m.Param(np.clip(-0.2*np.sin(2*np.pi*t),0,None))
d1v = m.Var(d1[0]); d2v=m.Intermediate(d1v*2)
d3v = m.Var(0,lb=0); d3v1=m.Var(0); d3v2=m.Var(0)
t1 = m.Intermediate(d1+d3v1); t2 = m.Intermediate(d2+d3v2)

e  = m.CV(0); e.STATUS=1; e.SPHI=e.SPLO=0; e.WSPHI=1000; e.WSPLO=1
h  = m.CV(0); h.STATUS=1; h.SPHI=h.SPLO=0; h.WSPHI=1000; h.WSPLO=1
z  = m.CV(0); z.STATUS=1; z.SPHI=z.SPLO=0; z.WSPHI=1000; z.WSPLO=1
r1 = m.MV(0,lb=-1,ub=1); r1.STATUS=1; r1.DCOST=0.0
r3 = m.MV(0,lb=-1,ub=1); r3.STATUS=1; r3.DCOST=0.0

m.Equations([d1v.dt()==r1, e==t1-d1v, h==t2-d2v])
m.Equations([d3v.dt()==r3, z==d3-d3v, d3v1==d3v*2, d3v2==d3v*3])
m.Maximize(d3v)

m.options.IMODE=6; m.options.SOLVER=1; m.solve()

plt.figure(figsize=(8,7))
plt.subplot(4,1,1)
plt.plot(t,t1,'r-',label='Demand 1')
plt.plot(t,d1v,'b:',label='Prod 1')
plt.legend(); plt.grid()
plt.subplot(4,1,2)
plt.plot(t,t2,'r-',label='Demand 2')
plt.plot(t,d2v,'b:',label='Prod 2')
plt.legend(); plt.grid()
plt.subplot(4,1,3)
plt.plot(t,d3,'r-',label='Demand 3')
plt.plot(t,d3v,'b:',label='Prod 3')
plt.legend(); plt.grid()
plt.subplot(4,1,4)
plt.plot(t,r1,'k:',label='Ramp Rate 1')
plt.plot(t,r3,'k--',label='Ramp Rate 3')
plt.legend(); plt.grid(); plt.xlabel('Time')
plt.show()