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.
The results I get:.
What am I doing wrong?
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.
The results I get:.
What am I doing wrong?
Share Improve this question edited Nov 21, 2024 at 21:43 John Hedengren 14.4k1 gold badge23 silver badges31 bronze badges asked Nov 21, 2024 at 14:20 iarimaiarima 453 bronze badges1 Answer
Reset to default 0There 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()