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 (°F)'); plt.legend(), plt.grid()
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.grid(); plt.tight_layout()
plt.savefig('furnace_control.png',dpi=300)
plt.show()