使用Gekko，目标函数不会收敛

Question

我试着运行这段代码，但即使我提高了迭代次数的上限，目标函数也会出现分歧(将其提高到2000次，但仍然没有找到解决方案(。

欢迎任何帮助。

The code:

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
m = GEKKO()
plt.rcParams['font.size'] = 15
# Parameters
nu_float = 0.5
beta_E_float = 10.0
tau_E_float=0.05
delta_E_float=0.03
beta_F_float=0.01
s_h_float = 0.9951
delta_F_float = 0.04
eta_float = 0.95
delta_float = 1.25
F_star = 5106.0
b_float = nu_float*beta_F_float*beta_E_float/(tau_E_float+delta_E_float)#3.125
print('b='+str(b_float))
K_float = F_star/(nu_float*beta_F_float/delta_F_float-nu_float*beta_F_float/b_float)
print('K='+str(K_float))
nu=m.Param(value=nu_float)
beta_E=m.Param(value=beta_E_float)
tau_E=m.Param(value=tau_E_float)
delta_E=m.Param(value=delta_E_float)
beta_F=m.Param(value=beta_F_float)
b=m.Param(value=b_float)
s_h=m.Param(value=s_h_float)
delta_F=m.Param(value=delta_F_float)
eta=m.Param(value=eta_float)
delta=m.Param(value=delta_float)
K=m.Param(value=K_float)
T = 80
C = 1000
U_bar = 150

# Time interval
nt = 500
m.time = np.linspace(0,T,nt)
delta_t = T/nt

# Control
u = m.MV(lb=0,ub=U_bar)
u.STATUS = 1
u.DCOST = 0

# Variables
E_u = m.Var(value=K*(1-delta_F/b))
F_u = m.Var(value=F_star)
E_i = m.Var(value=0.0)
F_i = m.Var(value=0.0)
U = m.Var(value=0.0)

# Vector to extract the final value
final_array = np.zeros(nt)
final_array[-1] = 1.0
final = m.Param(value=final_array)

# Equations
m.Equation(E_u.dt()==beta_E*F_u*(1-s_h*F_i/(F_u+F_i))*(1-(E_u+E_i)/K)-(tau_E+delta_E)*E_u)
m.Equation(F_u.dt()==nu*beta_F*E_u-delta_F*F_u)
m.Equation(E_i.dt()==eta*beta_E*F_i*(1-(E_u+E_i)/K)-(tau_E+delta_E)*E_i)
m.Equation(F_i.dt()==nu*beta_F*E_i-delta*delta_F*F_i+u)
m.Equation(U.dt()==u)
m.Equation(final*U<=C)

# Objective Function
m.Obj(0.5*(final*E_u)**2+0.5*(0.5*(abs(K*(1-delta*delta_F/(b*eta))-final*E_i)+K*(1-delta*delta_F/(b*eta))-final*E_i))**2+0.5*(final*F_u)**2+0.5*(0.5*(abs(K*(nu*beta_F/(delta*delta_F)-nu*beta_F/(b*eta))-final*F_i)+K*(nu*beta_F/(delta*delta_F)-nu*beta_F/(b*eta))-final*F_i))**2)

# Resolution
m.options.IMODE = 6
m.options.NODES = 4
m.options.MV_TYPE = 1
m.options.SOLVER = 3
print('beginning of computations')
m.solve()
print('end of computations')

# Print results
E_u_array=np.transpose(np.matrix(E_u))
F_u_array=np.transpose(np.matrix(F_u))
E_i_array=np.transpose(np.matrix(E_i))
F_i_array=np.transpose(np.matrix(F_i))
u_array=np.transpose(np.matrix(u))
U_array=np.transpose(np.matrix(U))
t_array = np.arange(nt)*T/nt
E_u_T = float(final_array*E_u_array)
F_u_T = float(final_array*F_u_array)
E_i_T = float(final_array*E_i_array)
F_i_T = float(final_array*F_i_array)
print(r'$E_u(T):$ ' + str(E_u_T))
print(r'$F_u(T):$ ' + str(F_u_T))
print(r'$E_i(T):$ ' + str(E_i_T))
print(r'$F_i(T):$ ' + str(F_i_T))

错误

An error occured.
The error code is           -1


---------------------------------------------------
Solver         :  IPOPT (v3.12)
Solution time  :    124.735700000005      sec
Objective      :    11931117306673.5
Unsuccessful with error code            0
---------------------------------------------------

Creating file: infeasibilities.txt
Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
@error: Solution Not Found
Traceback (most recent call last):
File "C:Usersdavetipe 2k22.py", line 86, in <module>
m.solve()
File "c:usersdaveappdatalocalprogramspythonpython37libsite-packagesgekkogekko.py", line 2185, in solve
raise Exception(response)
Exception:  @error: Solution Not Found

Answer 1

目标函数中的abs()函数在x=0处不可微。对于连续可微的二进制变量开关，请尝试切换到m.abs3()。

# Objective Function
m.Minimize(0.5*(final*E_u)**2
+0.5*(0.5*(m.abs3(K*(1-delta*delta_F/(b*eta))-final*E_i)
+K*(1-delta*delta_F/(b*eta))-final*E_i))**2+0.5*(final*F_u)**2
+0.5*(0.5*(m.abs3(K*(nu*beta_F/(delta*delta_F)
-nu*beta_F/(b*eta))-final*F_i)+K*(nu*beta_F/(delta*delta_F)
-nu*beta_F/(b*eta))-final*F_i))**2)

您还需要切换到求解器APOPT，以使用m.options.SOLVER=1计算混合整数非线性规划(MINLP(解决方案。如果省略m.options.SOLVER行，则Gekko会在使用m.abs3()功能时自动选择它。如果需要，可以使用另一个解算器选择来覆盖它，但它不一定会返回整数解。有关abs()函数二进制形式的更多信息，请参阅优化在线课程(优化中的逻辑条件(。

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
m = GEKKO()
plt.rcParams['font.size'] = 15
# Parameters
nu_float = 0.5
beta_E_float = 10.0
tau_E_float=0.05
delta_E_float=0.03
beta_F_float=0.01
s_h_float = 0.9951
delta_F_float = 0.04
eta_float = 0.95
delta_float = 1.25
F_star = 5106.0
b_float = nu_float*beta_F_float*beta_E_float/(tau_E_float+delta_E_float)#3.125
print('b='+str(b_float))
K_float = F_star/(nu_float*beta_F_float/delta_F_float-nu_float*beta_F_float/b_float)
print('K='+str(K_float))
nu=m.Param(value=nu_float)
beta_E=m.Param(value=beta_E_float)
tau_E=m.Param(value=tau_E_float)
delta_E=m.Param(value=delta_E_float)
beta_F=m.Param(value=beta_F_float)
b=m.Param(value=b_float)
s_h=m.Param(value=s_h_float)
delta_F=m.Param(value=delta_F_float)
eta=m.Param(value=eta_float)
delta=m.Param(value=delta_float)
K=m.Param(value=K_float)
T = 80
C = 1000
U_bar = 150
# Time interval
nt = 100
m.time = np.linspace(0,T,nt)
delta_t = T/nt
# Control
u = m.MV(lb=0,ub=U_bar)
u.STATUS = 1
u.DCOST = 0
# Variables
E_u = m.Var(value=K*(1-delta_F/b))
F_u = m.Var(value=F_star)
E_i = m.Var(value=0.0)
F_i = m.Var(value=0.0)
U = m.Var(value=0.0)
# Vector to extract the final value
final_array = np.zeros(nt)
final_array[-1] = 1.0
final = m.Param(value=final_array)
# Equations
m.Equation(E_u.dt()==beta_E*F_u*(1-s_h*F_i/(F_u+F_i))
*(1-(E_u+E_i)/K)-(tau_E+delta_E)*E_u)
m.Equation(F_u.dt()==nu*beta_F*E_u-delta_F*F_u)
m.Equation(E_i.dt()==eta*beta_E*F_i*(1-(E_u+E_i)/K)
-(tau_E+delta_E)*E_i)
m.Equation(F_i.dt()==nu*beta_F*E_i-delta*delta_F*F_i+u)
m.Equation(U.dt()==u)
m.Equation(final*U<=C)
# Objective Function
m.Minimize(0.5*(final*E_u)**2
+0.5*(0.5*(m.abs3(K*(1-delta*delta_F/(b*eta))-final*E_i)
+K*(1-delta*delta_F/(b*eta))-final*E_i))**2+0.5*(final*F_u)**2
+0.5*(0.5*(m.abs3(K*(nu*beta_F/(delta*delta_F)
-nu*beta_F/(b*eta))-final*F_i)+K*(nu*beta_F/(delta*delta_F)
-nu*beta_F/(b*eta))-final*F_i))**2)
# Resolution
m.options.IMODE = 6
m.options.NODES = 3
m.options.MV_TYPE = 1
m.options.SOLVER = 1
print('beginning of computations')
m.solve()
print('end of computations')
# Print results
E_u_array=np.transpose(np.matrix(E_u))
F_u_array=np.transpose(np.matrix(F_u))
E_i_array=np.transpose(np.matrix(E_i))
F_i_array=np.transpose(np.matrix(F_i))
u_array=np.transpose(np.matrix(u))
U_array=np.transpose(np.matrix(U))
t_array = np.arange(nt)*T/nt
E_u_T = float(final_array*E_u_array)
F_u_T = float(final_array*F_u_array)
E_i_T = float(final_array*E_i_array)
F_i_T = float(final_array*F_i_array)
print(r'$E_u(T):$ ' + str(E_u_T))
print(r'$F_u(T):$ ' + str(F_u_T))
print(r'$E_i(T):$ ' + str(E_i_T))
print(r'$F_i(T):$ ' + str(F_i_T))

这产生了一个成功的解决方案：

--------- APM Model Size ------------
Each time step contains
Objects      :            0
Constants    :            0
Variables    :           27
Intermediates:            0
Connections  :            0
Equations    :           13
Residuals    :           13

Number of state variables:           4356
Number of total equations: -         3861
Number of slack variables: -          990
---------------------------------------
Degrees of freedom       :           -495

* Warning: DOF <= 0
----------------------------------------------
Dynamic Control with APOPT Solver
----------------------------------------------
Iter:     1 I:  0 Tm:      5.72 NLPi:    7 Dpth:    0 Lvs:    0 Obj:  1.61E+11 Gap:  0.00E+00
Successful solution

---------------------------------------------------
Solver         :  APOPT (v1.0)
Solution time  :    5.74739999999292      sec
Objective      :    160555572371.779     
Successful solution
---------------------------------------------------

end of computations
$E_u(T):$ 40801.053817
$F_u(T):$ 5092.8119666
$E_i(T):$ 39.66478887
$F_i(T):$ 63.713604115