import numpy as np import matplotlib.pyplot as plt from gekko import GEKKO m = GEKKO() nt = 101 m.time = np.linspace(0,1,nt) # Parameters T = m.MV(value=362,ub=398,lb=298) T.STATUS = 1 T.DCOST = 0 # Variables x1 = m.Var(value=1) x2 = m.Var(value=0) p = np.zeros(nt) p[-1] = 1.0 final = m.Param(value=p) # Intermediates k1 = m.Intermediate(4000*m.exp(-2500/T)) k2 = m.Intermediate(6.2e5*m.exp(-5000/T)) # Equations m.Equation(x1.dt()==-k1*x1**2) m.Equation(x2.dt()==k1*x1**2 - k2*x2) # Objective Function m.Obj(-x2*final) m.options.IMODE = 6 m.solve() print('Objective: ' + str(x2[-1])) plt.figure(1) plt.subplot(2,1,1) plt.plot(m.time,x1.value,'k:',lw=2,label=r'$x_1$') plt.plot(m.time,x2.value,'b-',lw=2,label=r'$x_2$') plt.ylabel('Value') plt.legend(loc='best') plt.subplot(2,1,2) plt.plot(m.time,T.value,'r--',lw=2,label=r'$T$') plt.legend(loc='best') plt.xlabel('Time') plt.ylabel('Value') plt.show()