pythondynamic-programmingsimulationodegekko
No solution found in GEKKO when modelling a tank level
I'm trying to simulate the level of a tank that has two inlet flows and one outlet. The idea is that after this will be used in a control problem. However, I can't get it to work with gekko but it did work with scypy odeint.
#set time
tmax=60*6
i = 0
t = np.linspace(i,tmax,int(tmax/10)+1) # minutes
#Assign values
d2_h0 = 13.377 #initial height m
inf1 = np.array([32.6354599 , 32.41882451, 32.08460871, 32.11487071, 32.71570587,
32.59923999, 31.66669464, 30.11240896, 29.31222725, 29.35761197,
29.62183634, 29.67505582, 29.24057325, 29.13853518, 29.48321724,
29.61703173, 29.49874306, 28.99679947, 29.24003156, 29.40070153,
29.70169004, 29.2913545 , 29.47371801, 29.91566467, 31.31636302,
31.6771698 , 31.65268326, 31.06637255, 31.39147377, 31.88083331,
32.59566625, 32.70952861, 32.78859075, 32.87391027, 32.97800064,
32.99872208, 33.02946218])
inf2 = np.array([66.91262309, 67.16797638, 67.77143351, 66.85663605, 67.43820954,
67.96041107, 68.7215627 , 68.91900635, 69.20062764, 68.29413096,
68.56461334, 67.67184957, 68.84806824, 67.61451467, 69.58069102,
71.284935 , 75.60562642, 74.83906555, 74.06419373, 71.20425161,
69.60981496, 69.45553589, 70.35860697, 71.17754873, 72.16390737,
72.0528005 , 72.49635569, 73.09021505, 72.7195816 , 71.9975001 ,
70.13828532, 71.11123403, 72.16157023, 73.27675883, 71.9024353 ,
71.17524719, 70.34394582])
eff = np.array(([110.97348786, 108.6726354 , 109.4272232 , 110.57080078,
114.20512136, 114.84948222, 113.96173604, 110.81165822,
110.4366506 , 111.61210887, 112.75804393, 111.23046112,
108.35852305, 108.21724955, 110.47168223, 112.10458374,
109.28511048, 107.31727092, 108.55026245, 111.30213165,
111.88119253, 110.62695313, 111.76373037, 115.09386699,
115.75547282, 113.47773488, 107.95795441, 106.46175893,
105.83562978, 109.9902064 , 110.59869131, 110.49962108,
109.35623678, 108.35690053, 107.0867513 , 104.34462484,
103.1198527 ]))
from gekko import GEKKO
#Create gekko model
m = GEKKO(remote=False)
m.time = t
qin1 = m.Param(value=inf1)
qin2 = m.Param(value=inf2)
Ac=m.Const(value = 226.98) # m2,Cross section Area
qout = m.Param(value=eff)
h1 = m.Var(value=d2_h0, lb=0)
m.Equation(h1.dt() == (qin1 + qin2 - qout)/Ac)
m.options.IMODE = 4
m.options.SOLVER = 3
# Solve the model
m.solve()
What I get is: Exception: @error: Solution Not Found
If I remove the lowerbound it is able to find a solution but it is not correct. Below I have a comparison with the real value and what I get with odeint.
What am I doing wrong?
Solution
There is no solution because the lower bound on tank level height is set to 0 with:
h1 = m.Var(value=d2_h0, lb=0)
When this constraint is removed, it is solved successfully but with an unrealistic solution with a negative height.
Overflow or complete drainage can be included in the model by adding slack variables:
# slack variables
s_in = m.Var(value=0,lb=0)
s_out = m.Var(value=0,lb=0) # slack variable
m.Minimize(1e-3*(s_in+s_out))
They are normally zero but can be used by the optimizer to maintain feasibility for overflow or complete drainage.
Here is the complete code with the slack variables and an upper and lower limit on tank fluid height.
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
#set time
tmax=60*6
i = 0
t = np.linspace(i,tmax,int(tmax/10)+1) # minutes
#Assign values
d2_h0 = 13.377 #initial height m
inf1 = np.array([32.6354599 , 32.41882451, 32.08460871, 32.11487071, 32.71570587,
32.59923999, 31.66669464, 30.11240896, 29.31222725, 29.35761197,
29.62183634, 29.67505582, 29.24057325, 29.13853518, 29.48321724,
29.61703173, 29.49874306, 28.99679947, 29.24003156, 29.40070153,
29.70169004, 29.2913545 , 29.47371801, 29.91566467, 31.31636302,
31.6771698 , 31.65268326, 31.06637255, 31.39147377, 31.88083331,
32.59566625, 32.70952861, 32.78859075, 32.87391027, 32.97800064,
32.99872208, 33.02946218])
inf2 = np.array([66.91262309, 67.16797638, 67.77143351, 66.85663605, 67.43820954,
67.96041107, 68.7215627 , 68.91900635, 69.20062764, 68.29413096,
68.56461334, 67.67184957, 68.84806824, 67.61451467, 69.58069102,
71.284935 , 75.60562642, 74.83906555, 74.06419373, 71.20425161,
69.60981496, 69.45553589, 70.35860697, 71.17754873, 72.16390737,
72.0528005 , 72.49635569, 73.09021505, 72.7195816 , 71.9975001 ,
70.13828532, 71.11123403, 72.16157023, 73.27675883, 71.9024353 ,
71.17524719, 70.34394582])
eff = np.array(([110.97348786, 108.6726354 , 109.4272232 , 110.57080078,
114.20512136, 114.84948222, 113.96173604, 110.81165822,
110.4366506 , 111.61210887, 112.75804393, 111.23046112,
108.35852305, 108.21724955, 110.47168223, 112.10458374,
109.28511048, 107.31727092, 108.55026245, 111.30213165,
111.88119253, 110.62695313, 111.76373037, 115.09386699,
115.75547282, 113.47773488, 107.95795441, 106.46175893,
105.83562978, 109.9902064 , 110.59869131, 110.49962108,
109.35623678, 108.35690053, 107.0867513 , 104.34462484,
103.1198527 ]))
from gekko import GEKKO
#Create gekko model
m = GEKKO(remote=False)
m.time = t
qin1 = m.Param(value=inf1)
qin2 = m.Param(value=inf2)
Ac=m.Const(value = 226.98) # m2,Cross section Area
qout = m.Param(value=eff)
h1 = m.Var(value=d2_h0, lb=0, ub=20)
# slack variables
s_in = m.Var(value=0,lb=0)
s_out = m.Var(value=0,lb=0) # slack variable
m.Minimize(1e-3*(s_in+s_out))
m.Equation(h1.dt() == (qin1 + qin2 - qout + s_in - s_out)/Ac)
m.options.IMODE = 6
m.options.SOLVER = 3
# Solve the model
m.solve()
plt.figure(figsize=(6,3.5))
plt.subplot(2,1,1)
plt.plot(m.time,h1,'k-',label='height')
plt.grid(); plt.legend(); plt.ylabel('height (m)')
plt.subplot(2,1,2)
plt.plot(m.time,qin1,'b-',label=r'$q_{in,1}$')
plt.plot(m.time,qin2,'k:',label=r'$q_{in,2}$')
plt.plot(m.time,qout,'r--',label=r'$q_{out}$')
plt.grid(); plt.legend(); plt.ylabel(r'flow ($m^3$/min)')
plt.xlabel('Time (min)')
plt.tight_layout()
plt.savefig('tank.png',dpi=300)
plt.show()