Including disturbance rejection in GEKKO MPC

Try setting h1.SP = 0 instead of 1 to set the objective to drive it back to 1. Also, these options are not needed if h4.TR_INIT = 0 to create a setpoint that is zero everywhere (not a reference trajectory).

h4.TAU = 1
h4.BIAS = 1
h4.FSTATUS = 1

I also added lower bounds of zero for each of the height CVs with lb=0. There is a move on v2 when the disturbance starts. It appears that the controller can't resolve this offset. Is this the correct response?

from gekko import GEKKO
import numpy 
import matplotlib.pyplot as plt  

t = numpy.linspace(0,300,100)
m = GEKKO(remote=False)
m.time = t

T1 = m.Param(value = 53.97272679974334)
T2 = m.Param(value = 48.06851424706475)
T3 = m.Param(value = 38.48651254747577)
T4 = m.Param(value = 31.018933652439845)
k1 = m.Param(value = 5.51)
k2 = m.Param(value = 6.58)
γ1bar = m.Param(value = 0.333)
γ2bar = m.Param(value = 0.307)
A1 = m.Param(value = 730)
A2 = m.Param(value = 730)
A3 = m.Param(value = 730)
A4 = m.Param(value = 730)
v1bar = m.Param(value = 60)
v2bar = m.Param(value = 60)

# Manipulated variable
v1 = m.MV(value=0, lb=0, ub=100)
v1.STATUS = 1  
#v1.DCOST = 0.1 
#v2.DMAX = 20   

v2 = m.MV(value=0, lb=0, ub=100)
v2.STATUS = 1  
#v2.DCOST = 0.1 
#v2.DMAX = 20   

γ1 = m.MV(value=0, lb=0, ub=1)
γ1.STATUS = 1  
#γ1.DCOST = 0.1 
#v2.DMAX = 20   

γ2 = m.MV(value=0, lb=0, ub=1)
γ2.STATUS = 1  
#γ2.DCOST = 0.1 
#v2.DMAX = 20   

d = numpy.ones(t.shape)
d[0:10] = 0

d22 = numpy.zeros(t.shape)

d1 = m.Param(name='d1', value=d)  # Disturbance
d2 = m.Param(name='d2', value=d22)  # Disturbance

m.options.CV_TYPE = 2 # squared error

# Controlled Variable
h1 = m.CV(value=0,lb=0)
h1.STATUS = 1  # add the SP to the objective
h1.SP = 0    # set point
h1.TR_INIT = 0 # set point trajectory
#h1.TAU = 1   # time constant of trajectory
#h1.BIAS = 1
#h1.FSTATUS = 1

h2 = m.CV(value=0,lb=0)
h2.STATUS = 1  # add the SP to the objective
h2.SP = 0  # set point
h2.TR_INIT = 0 # set point trajectory
#h2.TAU = 1     # time constant of trajectory
#h2.BIAS = 1
#h2.FSTATUS = 1

h3 = m.CV(value=0,lb=0)
h3.STATUS = 1  # add the SP to the objective
h3.SP = 0     # set point
h3.TR_INIT = 0 # set point trajectory
#h3.TAU = 1     # time constant of trajectory
#h3.BIAS = 1
#h3.FSTATUS = 1

h4 = m.CV(value=0,lb=0)
h4.STATUS = 1  # add the SP to the objective
h4.SP = 0     # set point
h4.TR_INIT = 0 # set point trajectory
#h4.TAU = 1     # time constant of trajectory
#h4.BIAS = 1
#h4.FSTATUS = 1

m.Equation(h1.dt() == -(1/T1)*h1 + (A3/(A1*T3))*h3 + (γ1bar*k1*v1)/A1 + (γ1*k1*v1bar)/A1)
m.Equation(h2.dt() == -(1/T2)*h2 + (A4/(A2*T4))*h4 + (γ2bar*k2*v2)/A2 + (γ2*k2*v2bar)/A2)
m.Equation(h3.dt() == -(1/T3)*h3 + ((1-γ2bar)*k2*v2)/A3 - k2*v2bar*γ2/A3 - (k1*d1)/A3) 
m.Equation(h4.dt() == -(1/T4)*h4 + ((1-γ1bar)*k1*v1)/A4 - k1*v1bar*γ1/A4 - (k2*d2)/A4)

m.options.IMODE = 6 # control

m.solve(disp=True)

plt.subplot(3,1,1)
plt.plot(m.time,d1.value,'k-',label='Disturbance')
plt.ylabel('Disturbance'); plt.grid(); plt.legend()

plt.subplot(3,1,2)
plt.plot(m.time,h1.value,'r-',label='Height 1')
plt.plot(m.time,h2.value,'g--',label='Height 2')
plt.plot(m.time,h3.value,'b:',label='Height 3')
plt.plot(m.time,h3.value,'k--',alpha=0.5,label='Height 4')
plt.ylabel('CVs'); plt.grid(); plt.legend()

plt.subplot(3,1,3)
plt.plot(m.time,v1.value,'r-',label='MV 1a')
plt.plot(m.time,v2.value,'g--',label='MV 2a')
plt.plot(m.time,γ1.value,'b:',label='MV 1b')
plt.plot(m.time,γ2.value,'k--',alpha=0.5,label='MV 2b')
plt.ylabel('MVs'); plt.grid(); plt.legend()

plt.show()