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()