from gekko import GEKKO import matplotlib.pyplot as plt m = GEKKO() # create GEKKO model m.time = [0,1,2,4,8,12,16,20,25,30,35,40,45,50,\ 55,60,65,70,75,80,85,90,95,100,105,110,\ 115,120,125,130,135,140,145,150] # change solver options m.solver_options = ['minlp_gap_tol 1.0e-2',\ 'minlp_maximum_iterations 10000',\ 'minlp_max_iter_with_int_sol 5000'] m.options.SOLVER = 1 m.options.IMODE = 6 # parameters Kp = 60 taup = 120 T0 = 40 # create heater binary variable u = m.MV(integer=True,lb=0,ub=1) u.dcost = 0.1 u.status = 1 # controlled variable (temperature) T = m.CV(value=T0) T.SPHI = 68 T.SPLO = 70 T.STATUS = 1 T.TR_INIT = 0 # first order equation m.Equation(taup * T.dt() == -(T-T0) + Kp * u) m.solve() # solve MILP # plot solution plt.subplot(2,1,1) plt.plot(m.time,T,'r-',label='Temperature') plt.plot([0,m.time[-1]],[T.SPHI,T.SPHI],'k:',label=r'$T_{sphi}$') plt.plot([0,m.time[-1]],[T.SPLO,T.SPLO],'k:',label=r'$T_{splo}$') plt.ylabel('T (degF)') plt.legend() plt.subplot(2,1,2) plt.plot(m.time,u,'b.-',label='Furnace (On/Off)') plt.ylabel('Heat (On/Off)') plt.xlabel('Time (min)') plt.legend() plt.show()