Question: 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?

Answer: 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()