Is it possible in GEKKO to use an adjusted Fixed Variable (FV) as a bound for a Manipulated Variable (MV)?

Question

I've been learning python and trying to get a handle on GEKKO by modelling a simple dynamic system where there is a storage vessel with an inlet and outlet mass flow. I am trying to optimize for minimum required storage given time-series inlet flow data, and a controlled outlet flowrate.

The inlet flowrate is defined on an hourly basis as parameter feed_rate, while the outlet flowrate is set via the control variable outlet_frac_flow which gives a fractional output, multiplied by fixed flowrate outlet_rate_ref, to maintain the mass within the storage vessel stored_mass between defined upper and lower bounds.

I have got the model to work fine if I fix the minimum and maximum bounds for stored_mass as constant values, with GEKKO solving over the time series to adjust the outlet rate. What I am trying to do now is add a layer of optimization to minimize the maximum storage volume required via new FV max_storage.

Ideally, I would like to use this FV as the upper bound of MV stored_mass, but it does not seem to be working.

When I run the below code, I get the output that Max storage: 73.94 tonnes vs limit of 60.0 tonnes. So the solver is minimizing the 'max_storage' variable, but it does not seem to be carried through into the upper bound set for stored_mass.

I thought maybe I had to use an Intermediate variable to pass max_storage into stored_mass as upper bound but when I tried that then I couldn't even get the model to solve.

Can anyone advise if what I am trying to do is even feasible for the solver that GEKKO uses? And if so where I might be going wrong?

Much appreciated!

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

# Initialise GEKKO model
m = GEKKO(remote=False)

# set initial, minimum and maximum tank storage (tonnes)
min_storage_init = 30
max_storage_init = 120

# initial storage set at 10 tonnes above minimum
init_storage = min_storage_init + 10

# set reference outlet flowrate (to be adjusted by solver)
outlet_rate_ref = 20

# set hourly feed rate into tank in tonnes/hr
feed_rate_input = [4.2, 11.2, 10.2, 5.64, 0.0, 16.92, 18.31, 12.67, 18.31, 18.49, 18.49, 18.49, 18.49, 18.49, 17.95, 15.13, 12.25, 6.61, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

  
# set time
m.time = np.linspace(0,len(feed_rate_input)-1,len(feed_rate_input))

# input feed_rate into model as parameter
feed_rate = m.Param(value=feed_rate_input, name = 'feed rate')

# add maximum storage as fixed variable

max_storage = m.FV(value = max_storage_init, lb = min_storage_init*2, ub = max_storage_init)
max_storage.STATUS = 1

# add outlet_frac_flow as manipulated variable, to adjust the fractional reference outlet rate

outlet_frac_flow = m.MV(value=0.6, lb = 0.3, ub = 1)

outlet_frac_flow.STATUS = 1
outlet_frac_flow.DCOST = 1 # give some cost to changes in rate
outlet_frac_flow.DMAX = 0.3 # max of 30% change per step
outlet_frac_flow.MV_STEP_HOR = 5 # max of 3 changes in any 24hr period


# add stored mass as controlled variable
# initial value = init_storage, upper bound = fixed variable max_storage
stored_mass = m.CV(value=init_storage, lb = min_storage_init, ub = max_storage, name = 'mass stored')
stored_mass.STATUS = 1

# set high setpoint to be 95% of fixed variable max_storage
stored_mass.SPHI = max_storage *0.95
# set low setpoint to be 5% higher than min_storage_init
stored_mass.SPLO = min_storage_init / 0.95

# define mass balance differential
m.Equation(stored_mass.dt() == feed_rate - outlet_frac_flow * outlet_rate_ref)

# try to minimise max_storage

m.Minimize(max_storage)

#solver options
m.options.CV_TYPE=1 
m.options.IMODE = 6 
m.options.NODES = 3 
m.options.SOLVER = 3

#solve model

m.solve(disp=False)

flowrate=[i for i in outlet_frac_flow]
for i in range (0,len(feed_rate_input)):
    flowrate[i]=outlet_frac_flow[i] * outlet_rate_ref

print('Max storage: ' + str(round(max(stored_mass),2)) + ' tonnes vs limit of ' + str(round(max_storage[0],2)) + ' tonnes')    
    
#Plot graphs

#plt.figure(1) # plot results
#plt.set_figheight(13)
#plt.set_figwidth(18)
plt.figure(figsize=(16, 8)) 
plt.subplot(2,1,1)
plt.plot(m.time ,feed_rate.value,'r-',label='Feed rate (te/hr)')
plt.plot(m.time ,stored_mass.value,'b-',label='Mass stored (te)')
plt.legend()
plt.xticks(np.arange(min(m.time), max(m.time)+1, 1))
plt.subplot(2,1,2)
plt.plot(m.time,flowrate,'k--',label='outlet flowrate (te/hr)')
plt.xlabel('Time (hours)')
plt.xticks(np.arange(min(m.time), max(m.time)+1, 1))
plt.legend()
plt.show()

Answer 1

Add an inequality equation to change the max_storage to be greater than the stored_mass value:

m.Equation(stored_mass<=max_storage)

The initialization of the CV upper bound or the stored_mass.SPHI only use the values of max_storage at initialization. They do not update when max_storage is updated by the solver. To include the constraint, add it as an inequality constraint.

If you need the SPHI to also use the updated value of 0.95*max_storage then define a new CV such as ratio = m.CV() with ratio.STATUS=1, ratio.SPHI=0.95, and m.Equation(ratio*max_storage==stored_mass). By not writing the equation as ratio==stored_mass/max_storage, it avoids a potential divide-by-zero problem with the solver.