from gekko import GEKKO
import numpy as np
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)
g1  = m.Var(d1[0]); g2 = m.Intermediate(g1*2)

J1  = m.CV(0); J1.STATUS=1; J1.SPHI=J1.SPLO=0; J1.WSPHI=1000; J1.WSPLO=1
J2  = m.CV(0); J2.STATUS=1; J2.SPHI=J2.SPLO=0; J2.WSPHI=1000; J2.WSPLO=1
r = m.MV(0,lb=-1,ub=1); r.STATUS=1

m.Equations([g1.dt()==r, J1==d1-g1, J2==d2-g2])

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

plt.figure(figsize=(8,5))
plt.subplot(3,1,1)
plt.plot(t,d1,'r-',label='Demand 1')
plt.plot(t,g1,'b:',label='Prod 1')
plt.grid(); plt.legend()
plt.subplot(3,1,2)
plt.plot(t,d2,'r-',label='Demand 2')
plt.plot(t,g2,'b:',label='Prod 2')
plt.grid(); plt.legend()
plt.subplot(3,1,3)
plt.plot(t,r,'k--',label='Ramp Rate')
plt.grid(); plt.legend(); plt.xlabel('Time')
plt.show()