python - 在 GEKKO 中执行动态模拟时出错

Question

我尝试为apmonitor 网站上动态优化课程的建模（ https://apmonitor.com/do/index.php/Main/ModelFormulation ）部分中可用的油藏模拟编写动态模拟。代码如下：

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

m = GEKKO()

#time variable initialization
m.time = np.linspace(0,1,13)

Vf_in= [0, 0, 0, 0]
#define constants
#Area of Reservoir / Lake (km2)
areas =np.array([13.4, 12.0, 384.5, 4400])

#Initial Volume of Reservoir / Lake (km3)
volumes = np.array([0.26, 0.18, 0.68, 22.0])

#flow coefficient
c_out = np.array([0.03, 0.015, 0.06, 0])

h0 = 1000*(volumes/areas)
Vout0 = c_out*np.sqrt(h0)
print("Vout0", Vout0)

#rate of consumption
Vuse  = [0.03, 0.05, 0.02,0.00]

#inlet flow rate
Vfi_1 = [0.13, 0.13, 0.13, 0.21, 0.21, 0.21, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13]

#parameter
Vf_in[0] = m.Param(value = V_fi1)
A = [m.Param(value =i) for i in areas]
c_ev = [m.Param(value = 1e-5) for i in range(4)]
c_ev[3] = 0.5e-5
print("C_ev", c_ev)
c_f = [m.Param (value = i) for i in c_out]

#variables
h = [m.Var(value = i) for i in h0]
V = [m.Var(value = i) for i in volumes]
Vout = [m.Var(value = i) for i in Vout0]
print("volumes", V)
print("Heights", h)

#intermediates
Vf_in[1:4] = [m.Intermediate(Vout[i]) for i in range(3)]
V_ev = [m.Intermediate(c_ev[i]*A[i]) for i in range(4)]

#equations
m.Equation([V[i].dt() == Vf_in[i] - Vout[i] - V_ev[i] - Vuse[i] for i in range(4)])
m.Equation([h[i] == 1000*(V[i]/A[i]) for i in range(4)])
m.Equation(Vout[i]**2 == c_f[i]**2 *h[i] for i in range(4))
m.options.IMODE = 4
m.solve(disp = False)

time = [x*12 for x in m.time]
plt.figure()
plt.plot(time, h1.value, label = 'h1')
plt.plot(time, h2.value, label = 'h2')
plt.plot(time, h3.value, label = 'h3')
plt.plot(time, h4.value, label = 'h4')
plt.xlabel('time (in months)')
plt.ylabel('height in metres')
plt.legend(loc = 'best')
plt.show()

我遇到了下面给出的错误：

Exception:  @error: Inequality Definition
 invalid inequalities: z > x < y
 <generatorobject<genexpr>at0x000001655c7ad900>
 STOPPING . . .

任何帮助，将不胜感激。谢谢 ！

标签： pythonnumpygekko

Answer 1

尝试使用Get Code页面上的链接，而不是从网站复制文本。有时网站代码包含可能导致问题的特殊字符。如果这不能解决问题，请尝试升级 gekko：

pip install gekko --upgrade

from __future__ import division
from gekko import GEKKO
import numpy as np

#Initial conditions
c = np.array([0.03,0.015,0.06,0])
areas = np.array([13.4, 12, 384.5, 4400])
V0 = np.array([0.26, 0.18, 0.68, 22])
h0 = 1000 * V0 / areas
Vout0 = c * np.sqrt(h0)
vin = [0.13,0.13,0.13,0.21,0.21,0.21,0.13,\
       0.13,0.13,0.13,0.13,0.13,0.13]
Vin = [0,0,0,0]

#Initialize model
m = GEKKO(remote=False)

#time array 
m.time = np.linspace(0,1,13)
#define constants
c = m.Array(m.Const,4,value=0)
c[0].value = 0.03
c[1].value = c[0] / 2
c[2].value = c[0] * 2
c[3].value = 0
Vuse = [0.03,0.05,0.02,0.00]

#Parameters
evap_c = m.Array(m.Param,4,value=1e-5)
evap_c[-1].value = 0.5e-5

A = [m.Param(value=i) for i in areas]

Vin[0] = m.Param(value=vin)

#Variables
V = [m.Var(value=i) for i in V0]
h = [m.Var(value=i) for i in h0]
Vout = [m.Var(value=i) for i in Vout0]

#Intermediates
Vin[1:4] = [m.Intermediate(Vout[i]) for i in range(3)]
Vevap = [m.Intermediate(evap_c[i] * A[i]) for i in range(4)]

#Equations
m.Equations([V[i].dt() == \
             Vin[i] - Vout[i] - Vevap[i] - Vuse[i] \
             for i in range(4)])
m.Equations([1000*V[i] == h[i]*A[i] for i in range(4)])
m.Equations([Vout[i]**2 == c[i]**2 * h[i] for i in range(4)])


#Set to simulation mode
m.options.imode = 4

#Solve
m.solve()

#%% Plot results
time = [x * 12 for x in m.time] 

# plot results
import matplotlib.pyplot as plt
plt.figure(1)

plt.subplot(311)
plt.plot(time,h[0].value,'r-')
plt.plot(time,h[1].value,'b--')
plt.ylabel('Level (m)')
plt.legend(['Jordanelle Reservoir','Deer Creek Reservoir'])

plt.subplot(312)
plt.plot(time,h[3].value,'g-')
plt.plot(time,h[2].value,'k:')
plt.ylabel('Level (m)')
plt.legend(['Great Salt Lake','Utah Lake'])

plt.subplot(313)
plt.plot(time,Vin[0].value,'k-')
plt.plot(time,Vout[0].value,'r-')
plt.plot(time,Vout[1].value,'b--')
plt.plot(time,Vout[2].value,'g-')
plt.xlabel('Time (month)')
plt.ylabel('Flow (km3/yr)')
plt.legend(['Supply Flow','Upper Provo River', \
            'Lower Provo River','Jordan River'])
plt.show()

EXIT: Optimal Solution Found.

 The solution was found.

 The final value of the objective function is  0.
 
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :  0.0534 sec
 Objective      :  0.
 Successful solution
 ---------------------------------------------------