Search code examples
gekko

Contraint on variable output sum

I am new to Gekko and I am trying to model a production process. The model is working as expected, but I am failing to understand how to implement an equation that constraints the minimum output of my output variable. My code:

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

maxspeed = 30 # rev/hr
minspeed = 2 #rev/hr
lowrev = 1/(maxspeed/3600) # how fast can we go as sec/rev
highrev = 1/(minspeed/3600) # ow slow can we go as sec/rev
batchsize = 0.5

m = GEKKO(remote=False)
m.time = np.linspace(0,200,201)

tfr = m.Param(value=54.2)
subf = m.Param(value=13.47)
drumd = m.Param(value=2.9)
a1m = m.Param(value=1814.28)
dl = m.Param(value=0.42) 

drumspeed = m.MV(value=400, lb=lowrev, ub=highrev) 
drumspeed.STATUS = 1 
drumspeed.DMAX = 2

filtrationvol = m.Var(value=0, lb=0)

m.Equation(filtrationvol == tfr*(m.sqrt(4*a1m*m.acos(1- 
2*dl/drumd)*drumspeed/drumd))/(2*a1m*subf*drumspeed)) # output model
m.Maximize(filtrationvol)
m.options.IMODE = 6 # control mode
m.solve(disp=True)

I want to limit the total output of variable filtrationvol to a maximum of 0.5, meaning I want to maximize the output but only up to a certain point. I have been playing around with m.sum and m.vsum but without success.

Solution

  • Use the integral function to get the total filtration volume.

    m.Equation(total_Vfilt == m.integral(filtrationvol))

    You can either include a hard constraint with total_Vfilt<0.5 or you can set up a soft constraint (objective function) with a CV type. A rearrangement of the equation makes it easier for the solver because it avoids m.sqrt() and a potential divide by zero if the drum speed ever needs to go to zero.

    m.Equation((2*a1m*subf*drumspeed*filtrationvol/tfr)**2 \
           ==(4*a1m*m.acos(1-2*dl/drumd)*drumspeed/drumd))

    An initialization strategy with COLDSTART=1 also helps the solver because it sets the STATUS=0 to simulate and obtain a better initial guess before optimizing.

    m.options.COLDSTART = 1
m.solve(disp=False)

m.options.TIME_SHIFT= 0
m.options.COLDSTART = 0
m.solve(disp=True)

    Filtration plot

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

maxspeed = 30 # rev/hr
minspeed = 2 #rev/hr
lowrev = 1/(maxspeed/3600) # how fast can we go as sec/rev
highrev = 1/(minspeed/3600) # ow slow can we go as sec/rev
batchsize = 0.5

m = GEKKO(remote=False)
m.time = np.linspace(0,200,201)

tfr = m.Param(value=54.2)
subf = m.Param(value=13.47)
drumd = m.Param(value=2.9)
a1m = m.Param(value=1814.28)
dl = m.Param(value=0.42) 

drumspeed = m.MV(value=400, lb=lowrev, ub=highrev) 
drumspeed.STATUS = 1 
drumspeed.DMAX = 2

filtrationvol = m.Var(value=0.0, lb=0)

m.Equation((2*a1m*subf*drumspeed*filtrationvol/tfr)**2 \
           ==(4*a1m*m.acos(1-2*dl/drumd)*drumspeed/drumd))
m.Maximize(filtrationvol)

# summation with an integral
total_Vfilt = m.CV(0,lb=0) # hard constraint with ub=0.5
# total_Vfilt.UPPER = 0.5  #  or adjust .UPPER
m.Equation(total_Vfilt == m.integral(filtrationvol))

# soft constraint with abs() value penalty for exceeding limit
total_Vfilt.SPHI = 0.5; total_Vfilt.SPLO = 0.0
total_Vfilt.WSPHI = 1000
total_Vfilt.STATUS = 1

m.options.IMODE = 6  # control mode
m.options.SOLVER = 3 # IPOPT solver

m.options.COLDSTART = 1
m.solve(disp=False)

m.options.TIME_SHIFT= 0
m.options.COLDSTART = 0
m.solve(disp=True)

plt.figure(figsize=(12,6))
plt.subplot(4,1,1)
plt.plot(m.time,drumspeed.value,'r-',label='Sdrum')
plt.legend()
plt.subplot(4,1,2)
plt.plot(m.time,drumspeed.value,'r-',label='Sdrum')
plt.plot([0,200],[lowrev,lowrev],'k:',label='Limits')
plt.plot([0,200],[highrev,highrev],'k:')
plt.legend()
plt.subplot(4,1,3)
plt.plot(m.time,filtrationvol.value,'b--',label='Vfilt')
plt.legend()
plt.subplot(4,1,4)
plt.plot(m.time,total_Vfilt.value,'b:',\
         label=r'$\int(Vfilt)$')
plt.plot([0,200],[0.5,0.5],'k-',label='Limit')
plt.legend()
plt.xlabel('Batch / Time')
plt.show()