Gekko：MINLP - 找不到错误 options.json 文件
Gekko: MINLP - Error options.json file not found
我正在尝试解决 MINLP 问题，首先使用 IPOPT 求解器获得初始解，然后使用 APOPT 获得混合整数解。但是，调用 APOPT 求解器时出现以下错误：
Error: Exception: Access Violation At line 359 of file ./f90/cqp.f90 Traceback: not available, compile with -ftrace=frame or -ftrace=full Error: 'results.json' not found. Check above for additional error details Traceback (most recent call last): File "optimisation_stack.py", line 244, in Optimise_G(t,ob, jofbuses, q, qc, s, oa, k, l, T, G_previous, C, Y, G_previous, G_max, G_min) File "optimisation_stack.py", line 134, in Optimise_G sol = MINLP(xinit, A, B, A_eq, B_eq, LB ,UB, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_previous, C, Y, G_previous) File "optimisation_stack.py", line 215, in MINLP m_APOPT.solve(disp = False) File "C:\Users\Zineb\AppData\Local\Programs\Python\Python37\lib\site-packages\gekko\gekko.py", line 2227, in solve self.load_JSON() File "C:\Users\Zineb\AppData\Local\Programs\Python\Python37\lib\site-packages\gekko\gk_post_solve.py", line 13, in load_JSON f = open(os.path.join(self._path,'options.json')) FileNotFoundError: [Errno 2] No such file or directory: 'C:\Users\Zineb\AppData\Local\Temp\tmptdgafg1zgk_model1\options.json'
我的代码如下，我尽量简化了:
import numpy as np
from gekko import GEKKO
# Define matrices A,A_eq, and vectors b, b_eq for the optimization
def Optimise_G(t,ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous, G_max, G_min):
Mbig_1 = T*C
Mbig_2 = C
nb_phases = len(G_next)
b_max = len(t)
no_lanegroups = len(q)
A_eq = np.zeros(((nb_phases+1)*b_max + 1, (3*nb_phases+3)*b_max+nb_phases))
for i in range(nb_phases):
A_eq[0][i] = 1
#B_eq = np.zeros(((nb_phases+1)*b_max + 1, 1))
B_eq = np.zeros((nb_phases+1)*b_max + 1)
B_eq[0] = C - sum(Y[0:nb_phases])
counter_eq = 0
# G(i)=Ga(i,b)+Gb(i,b)+Gc(i,b)
for b in range(b_max):
for i in range(nb_phases):
counter_eq = counter_eq + 1
A_eq[counter_eq][i] = 1
A_eq[counter_eq][nb_phases*(b+1)+ i] = -1
A_eq[counter_eq][nb_phases*b_max + nb_phases*(b+1) + i] = -1
A_eq[counter_eq][2*nb_phases*b_max + nb_phases*(b+1) + i] = -1
# ya(b)+y(b)+y(c)=1
for b in range(b_max):
counter_eq = counter_eq + 1
A_eq[counter_eq][3*nb_phases*b_max + nb_phases + b] = 1
A_eq[counter_eq][(3*nb_phases+1)*b_max + nb_phases + b] = 1
A_eq[counter_eq][(3*nb_phases+2)*b_max + nb_phases + b] = 1
B_eq[counter_eq] = 1
A = np.zeros((no_lanegroups + (2*3*nb_phases+4)*b_max, (3*nb_phases+3)*b_max+nb_phases))
B = np.zeros(no_lanegroups + (2*3*nb_phases+4)*b_max)
counter = -1
# Sum Gi (i in Ij)>=Gj,min
for j in range(no_lanegroups):
counter = counter + 1
for i in range(k[j], l[j]+1):
A[counter][i-1] = -1
B[counter] = -C*qc[j]/s[j]
# ya(b)G_lb(i)<=Ga(i,b), yb(b)G_lb(i)<=Gb(i,b), yc(b)G_lb(i)<=Gc(i,b)
for b in range(b_max):
for i in range(nb_phases):
counter = counter + 1
A[counter][nb_phases*(b+1)+i] = -1
A[counter][3*nb_phases*b_max + nb_phases + b] = G_min[i]
B[counter] = 0
counter = counter + 1
A[counter][nb_phases*b_max + nb_phases*(b+1) + i] = -1
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = G_min[i]
B[counter] = 0
counter = counter + 1
A[counter][2*nb_phases*b_max + nb_phases*(b+1) +i] = -1
A[counter][(3*nb_phases+2)*b_max + nb_phases + b] = G_min[i]
B[counter] = 0
# ya(b)Gmax(i)>=Ga(i,b), yb(b)Gmax(i)>=Gb(i,b), yc(b)Gmax(i)>=Gc(i,b)
for b in range(b_max):
for i in range(nb_phases):
counter = counter + 1
A[counter][nb_phases*(b+1) +i] = 1
A[counter][3*nb_phases*b_max + nb_phases + b] = -G_max[i]
B[counter] = 0
counter = counter + 1
A[counter][nb_phases*b_max + nb_phases*(b+1) + i] = 1
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = -G_max[i]
B[counter] = 0
counter = counter + 1
A[counter][2*nb_phases*b_max + nb_phases*(b+1) +i] = 1
A[counter][(3*nb_phases+2)*b_max + nb_phases + b] = -G_max[i]
B[counter] = 0
# (1-yc(b))t(b)<=(T-1)C+sum(Gi(1:l(jofbuses(b))))+sum(Y(1:l(jofbuses(b))-1))
for b in range(b_max):
counter = counter + 1
A[counter][0:l[jofbuses[b]-1]] = -np.ones((1,l[jofbuses[b]-1]))
A[counter][(3*nb_phases+2)*b_max+nb_phases+b] = -t[b]
B[counter] = -t[b] + (T-1)*C + sum(Y[0:l[jofbuses[b]-1]-1])
# (T-1)C+sum(Gi(1:l(jofbuses(b))))+sum(Y(1:l(jofbuses(b))-1))<=yc(b)t(b)+(1-yc(b))Mbig_1
for b in range(b_max):
counter = counter + 1
A[counter][0:l[jofbuses[b]-1]] = np.ones((1,l[jofbuses[b]-1]))
A[counter][(3*nb_phases+2)*b_max+nb_phases+b] = -t[b] + Mbig_1
B[counter] = Mbig_1 - (T-1)*C - sum(Y[0:l[jofbuses[b]-1]-1])
# -Mbig_2(1-yb(b))<=db(b)=right-hand side of Equation (6)
for b in range(b_max):
counter = counter + 1
constant = q[jofbuses[b]-1]/s[jofbuses[b]-1]*(t[b] - (T-1)*C + sum(G_previous[l[jofbuses[b]-1]:nb_phases]) + sum(Y[l[jofbuses[b]-1] -1:nb_phases]))+ (T-1)*C + sum(Y[0:k[jofbuses[b]-1]-1]) - t[b]
A[counter][0:k[jofbuses[b]-1]-1] = -np.ones((1,k[jofbuses[b]-1]-1))
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = Mbig_2
B[counter] = constant + Mbig_2
# db(b)<=Mbig_2 yb(b)
for b in range(b_max):
counter = counter + 1
constant = q[jofbuses[b]-1]/s[jofbuses[b]-1]*(t[b] - (T-1)*C +sum(G_previous[l[jofbuses[b]-1]:nb_phases]) + sum(Y[l[jofbuses[b]-1] -1:nb_phases]))+ (T-1)*C + sum(Y[0:k[jofbuses[b]-1]-1]) - t[b]
A[counter][0:k[jofbuses[b]-1]-1] = np.ones((1,k[jofbuses[b]-1]-1))
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = -Mbig_2
B[counter] = -constant
#Lower Bound LB
LB_zeros = np.zeros(3*b_max*(nb_phases+1))
G_min = np.array(G_min)
LB = np.append(G_min, LB_zeros)
#Upper Bound UB
UB = np.ones(3*b_max)
G_max = np.array(G_max)
for i in range(3*b_max+1):
UB = np.concatenate((G_max,UB))
xinit = np.array([(a+b)/2 for a, b in zip(UB, LB)])
sol = MINLP(xinit, A, B, A_eq, B_eq, LB ,UB, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_previous, C, Y, G_previous)
objective函数：
def objective_fun(x, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous):
nb_phases = len(G_next)
b_max = len(t)
no_lanegroups = len(q)
obj = 0
obj_a = 0
obj_b = 0
G = x[0:nb_phases]
for j in range(no_lanegroups):
delay_a = 0.5*q[j]/(1-q[j]/s[j]) * (pow((sum(G_previous[l[j]:nb_phases]) + sum(G[0:k[j]-1]) + sum(Y[l[j]-1:nb_phases]) + sum(Y[0:k[j]-1])),2) + pow(sum(G[l[j]:nb_phases]) + sum(G_next[0:k[j]-1]) + sum(Y[l[j]-1:nb_phases]) + sum(Y[0:k[j]-1]),2))
obj = obj + oa*delay_a
obj_a = obj_a + oa*delay_a
for b in range(b_max):
delay_b1 = x[(3*nb_phases+1)*b_max + nb_phases + b]*(q[jofbuses[b]-1]/s[jofbuses[b]-1] * (t[b] - (T-1)*C + sum(G_previous[l[jofbuses[b]-1]:nb_phases]) + sum(Y[l[jofbuses[b]-1] -1:nb_phases])) + (T-1)*C - t[b] + sum(Y[0:k[jofbuses[b]-1]-1]))
delay_b2 = x[(3*nb_phases+2)*b_max + nb_phases + b-1]*(q[jofbuses[b]-1]/s[jofbuses[b]-1] * (t[b] - (T-1)*C - sum(Y[0:l[jofbuses[b]-1]-1])) + T*C + sum(G_next[0:k[jofbuses[b]-1]-1]) + sum(Y[0:k[jofbuses[b]-1]-1]) - t[b])
delay_b3 = sum(x[nb_phases*b_max + nb_phases*b:nb_phases*b_max + nb_phases*b+k[jofbuses[b]-1]-1]) - q[jofbuses[b]-1]/s[jofbuses[b]-1]*sum(x[2*nb_phases*b_max + nb_phases*b:2*nb_phases*b_max + nb_phases*b +l[jofbuses[b]-1]])
delay_b = delay_b1+delay_b2 +delay_b3
obj = obj + delay_b*ob[b]
obj_b = obj_b + delay_b*ob[b]
return obj
MINLP 求解器：
def MINLP(xinit, A, B, A_eq, B_eq, LB ,UB, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous):
nb_phases = len(G_next)
b_max = len(t)
## First Solver: IPOPT to get an initial guess
m_IPOPT = GEKKO(remote = False)
m_IPOPT.options.SOLVER = 3 #(IPOPT)
# Array Variable
rows = nb_phases + 3*b_max*(nb_phases+1)#48
x_initial = np.empty(rows,dtype=object)
for i in range(3*nb_phases*b_max+nb_phases+1):
x_initial[i] = m_IPOPT.Var(value = xinit[i], lb = LB[i], ub = UB[i])#, integer = False)
for i in range(3*nb_phases*b_max+nb_phases+1, (3*nb_phases+3)*b_max+nb_phases):
x_initial[i] = m_IPOPT.Var(value = xinit[i], lb = LB[i], ub = UB[i])#, integer = True)
# Constraints
m_IPOPT.axb(A,B,x_initial,etype = '<=',sparse=False)
m_IPOPT.axb(A_eq,B_eq,x_initial,etype = '=',sparse=False)
# Objective Function
f = objective_fun(x_initial, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous)
m_IPOPT.Obj(f)
#Solver
m_IPOPT.solve(disp = True)
####################################################################################################
## Second Solver: APOPT to solve MINLP
m_APOPT = GEKKO(remote = False)
m_APOPT.options.SOLVER = 1 #(APOPT)
x = np.empty(rows,dtype=object)
for i in range(3*nb_phases*b_max+nb_phases+1):
x[i] = m_APOPT.Var(value = x_initial[i], lb = LB[i], ub = UB[i], integer = False)
for i in range(3*nb_phases*b_max+nb_phases+1, (3*nb_phases+3)*b_max+nb_phases):
x[i] = m_APOPT.Var(value = x_initial[i], lb = LB[i], ub = UB[i], integer = True)
# Constraints
m_APOPT.axb(A,B,x,etype = '<=',sparse=False)
m_APOPT.axb(A_eq,B_eq,x,etype = '=',sparse=False)
# Objective Function
f = objective_fun(x, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous)
m_APOPT.Obj(f)
#Solver
m_APOPT.solve(disp = False)
return x
定义参数并调用Optimise_G:
#Define Parameters
C = 120
T = 31
b_max = 2
G_base = [12,31,12,11,1,41]
G_max = [106, 106, 106, 106, 106, 106]
G_min = [7,3,7,10,0,7]
G_previous = [7.3333333, 34.16763, 7.1333333, 10.0, 2.2008602e-16, 47.365703]
Y = [2, 2, 3, 2, 2, 3]
jofbuses = [9,3]
k = [3,1,6,4,1,2,5,4,2,3]
l = [3,2,6,5,1,3,6,4,3,4]
nb_phases = 6
oa = 1.25
ob = [39,42]
t = [3600.2, 3603.5]
q = [0.038053758888888886, 0.215206065, 0.11325116416666667, 0.06299876472222223,0.02800455611111111,0.18878488361111112,0.2970903402777778, 0.01876728472222222, 0.2192723663888889, 0.06132227222222222]
qc = [0.04083333333333333, 0.2388888888888889, 0.10555555555555556, 0.0525, 0.030555555555555555, 0.20444444444444446,0.31083333333333335, 0.018333333333333333, 0.12777777777777777, 0.07138888888888889]
s = [1.0, 1.0, 1.0, 1.0, 0.5, 1.0, 1.0, 0.5, 0.5, 0.5]
nb_phases = len(G_base)
G_max = []
for i in range(nb_phases):
G_max.append(C - sum(Y[0:nb_phases]))
Optimise_G(t,ob, jofbuses, q, qc, s, oa, k, l, T, G_previous, C, Y, G_previous, G_max, G_min)
有办法解决这个问题吗？
非常感谢！
Windows 版本的 APOPT 求解器崩溃，无法找到解决方案。不过，在线Linux版本的APOPT是可以找到解决办法的。获取 Gekko (v1.0.0 pre-release) available on GitHub 的最新版本。当新版本发布时，这将在 pip install gekko --upgrade 中可用，但现在您需要将源代码复制到 Lib\site-packages\gekko\gekko.py。更新gekko后，切换到remote=True如下图。
import numpy as np
from gekko import GEKKO
# Define matrices A,A_eq, and vectors b, b_eq for the optimization
def Optimise_G(t,ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous, G_max, G_min):
Mbig_1 = T*C
Mbig_2 = C
nb_phases = len(G_next)
b_max = len(t)
no_lanegroups = len(q)
A_eq = np.zeros(((nb_phases+1)*b_max + 1, (3*nb_phases+3)*b_max+nb_phases))
for i in range(nb_phases):
A_eq[0][i] = 1
#B_eq = np.zeros(((nb_phases+1)*b_max + 1, 1))
B_eq = np.zeros((nb_phases+1)*b_max + 1)
B_eq[0] = C - sum(Y[0:nb_phases])
counter_eq = 0
# G(i)=Ga(i,b)+Gb(i,b)+Gc(i,b)
for b in range(b_max):
for i in range(nb_phases):
counter_eq = counter_eq + 1
A_eq[counter_eq][i] = 1
A_eq[counter_eq][nb_phases*(b+1)+ i] = -1
A_eq[counter_eq][nb_phases*b_max + nb_phases*(b+1) + i] = -1
A_eq[counter_eq][2*nb_phases*b_max + nb_phases*(b+1) + i] = -1
# ya(b)+y(b)+y(c)=1
for b in range(b_max):
counter_eq = counter_eq + 1
A_eq[counter_eq][3*nb_phases*b_max + nb_phases + b] = 1
A_eq[counter_eq][(3*nb_phases+1)*b_max + nb_phases + b] = 1
A_eq[counter_eq][(3*nb_phases+2)*b_max + nb_phases + b] = 1
B_eq[counter_eq] = 1
A = np.zeros((no_lanegroups + (2*3*nb_phases+4)*b_max, (3*nb_phases+3)*b_max+nb_phases))
B = np.zeros(no_lanegroups + (2*3*nb_phases+4)*b_max)
counter = -1
# Sum Gi (i in Ij)>=Gj,min
for j in range(no_lanegroups):
counter = counter + 1
for i in range(k[j], l[j]+1):
A[counter][i-1] = -1
B[counter] = -C*qc[j]/s[j]
# ya(b)G_lb(i)<=Ga(i,b), yb(b)G_lb(i)<=Gb(i,b), yc(b)G_lb(i)<=Gc(i,b)
for b in range(b_max):
for i in range(nb_phases):
counter = counter + 1
A[counter][nb_phases*(b+1)+i] = -1
A[counter][3*nb_phases*b_max + nb_phases + b] = G_min[i]
B[counter] = 0
counter = counter + 1
A[counter][nb_phases*b_max + nb_phases*(b+1) + i] = -1
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = G_min[i]
B[counter] = 0
counter = counter + 1
A[counter][2*nb_phases*b_max + nb_phases*(b+1) +i] = -1
A[counter][(3*nb_phases+2)*b_max + nb_phases + b] = G_min[i]
B[counter] = 0
# ya(b)Gmax(i)>=Ga(i,b), yb(b)Gmax(i)>=Gb(i,b), yc(b)Gmax(i)>=Gc(i,b)
for b in range(b_max):
for i in range(nb_phases):
counter = counter + 1
A[counter][nb_phases*(b+1) +i] = 1
A[counter][3*nb_phases*b_max + nb_phases + b] = -G_max[i]
B[counter] = 0
counter = counter + 1
A[counter][nb_phases*b_max + nb_phases*(b+1) + i] = 1
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = -G_max[i]
B[counter] = 0
counter = counter + 1
A[counter][2*nb_phases*b_max + nb_phases*(b+1) +i] = 1
A[counter][(3*nb_phases+2)*b_max + nb_phases + b] = -G_max[i]
B[counter] = 0
# (1-yc(b))t(b)<=(T-1)C+sum(Gi(1:l(jofbuses(b))))+sum(Y(1:l(jofbuses(b))-1))
for b in range(b_max):
counter = counter + 1
A[counter][0:l[jofbuses[b]-1]] = -np.ones((1,l[jofbuses[b]-1]))
A[counter][(3*nb_phases+2)*b_max+nb_phases+b] = -t[b]
B[counter] = -t[b] + (T-1)*C + sum(Y[0:l[jofbuses[b]-1]-1])
# (T-1)C+sum(Gi(1:l(jofbuses(b))))+sum(Y(1:l(jofbuses(b))-1))<=yc(b)t(b)+(1-yc(b))Mbig_1
for b in range(b_max):
counter = counter + 1
A[counter][0:l[jofbuses[b]-1]] = np.ones((1,l[jofbuses[b]-1]))
A[counter][(3*nb_phases+2)*b_max+nb_phases+b] = -t[b] + Mbig_1
B[counter] = Mbig_1 - (T-1)*C - sum(Y[0:l[jofbuses[b]-1]-1])
# -Mbig_2(1-yb(b))<=db(b)=right-hand side of Equation (6)
for b in range(b_max):
counter = counter + 1
constant = q[jofbuses[b]-1]/s[jofbuses[b]-1]*(t[b] - (T-1)*C + sum(G_previous[l[jofbuses[b]-1]:nb_phases]) + sum(Y[l[jofbuses[b]-1] -1:nb_phases]))+ (T-1)*C + sum(Y[0:k[jofbuses[b]-1]-1]) - t[b]
A[counter][0:k[jofbuses[b]-1]-1] = -np.ones((1,k[jofbuses[b]-1]-1))
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = Mbig_2
B[counter] = constant + Mbig_2
# db(b)<=Mbig_2 yb(b)
for b in range(b_max):
counter = counter + 1
constant = q[jofbuses[b]-1]/s[jofbuses[b]-1]*(t[b] - (T-1)*C +sum(G_previous[l[jofbuses[b]-1]:nb_phases]) + sum(Y[l[jofbuses[b]-1] -1:nb_phases]))+ (T-1)*C + sum(Y[0:k[jofbuses[b]-1]-1]) - t[b]
A[counter][0:k[jofbuses[b]-1]-1] = np.ones((1,k[jofbuses[b]-1]-1))
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = -Mbig_2
B[counter] = -constant
#Lower Bound LB
LB_zeros = np.zeros(3*b_max*(nb_phases+1))
G_min = np.array(G_min)
LB = np.append(G_min, LB_zeros)
#Upper Bound UB
UB = np.ones(3*b_max)
G_max = np.array(G_max)
for i in range(3*b_max+1):
UB = np.concatenate((G_max,UB))
xinit = np.array([(a+b)/2 for a, b in zip(UB, LB)])
sol = MINLP(xinit, A, B, A_eq, B_eq, LB ,UB, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_previous, C, Y, G_previous)
def objective_fun(x, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous):
nb_phases = len(G_next)
b_max = len(t)
no_lanegroups = len(q)
obj = 0
obj_a = 0
obj_b = 0
G = x[0:nb_phases]
for j in range(no_lanegroups):
delay_a = 0.5*q[j]/(1-q[j]/s[j]) * (pow((sum(G_previous[l[j]:nb_phases]) + sum(G[0:k[j]-1]) + sum(Y[l[j]-1:nb_phases]) + sum(Y[0:k[j]-1])),2) + pow(sum(G[l[j]:nb_phases]) + sum(G_next[0:k[j]-1]) + sum(Y[l[j]-1:nb_phases]) + sum(Y[0:k[j]-1]),2))
obj = obj + oa*delay_a
obj_a = obj_a + oa*delay_a
for b in range(b_max):
delay_b1 = x[(3*nb_phases+1)*b_max + nb_phases + b]*(q[jofbuses[b]-1]/s[jofbuses[b]-1] * (t[b] - (T-1)*C + sum(G_previous[l[jofbuses[b]-1]:nb_phases]) + sum(Y[l[jofbuses[b]-1] -1:nb_phases])) + (T-1)*C - t[b] + sum(Y[0:k[jofbuses[b]-1]-1]))
delay_b2 = x[(3*nb_phases+2)*b_max + nb_phases + b-1]*(q[jofbuses[b]-1]/s[jofbuses[b]-1] * (t[b] - (T-1)*C - sum(Y[0:l[jofbuses[b]-1]-1])) + T*C + sum(G_next[0:k[jofbuses[b]-1]-1]) + sum(Y[0:k[jofbuses[b]-1]-1]) - t[b])
delay_b3 = sum(x[nb_phases*b_max + nb_phases*b:nb_phases*b_max + nb_phases*b+k[jofbuses[b]-1]-1]) - q[jofbuses[b]-1]/s[jofbuses[b]-1]*sum(x[2*nb_phases*b_max + nb_phases*b:2*nb_phases*b_max + nb_phases*b +l[jofbuses[b]-1]])
delay_b = delay_b1+delay_b2 +delay_b3
obj = obj + delay_b*ob[b]
obj_b = obj_b + delay_b*ob[b]
return obj
def MINLP(xinit, A, B, A_eq, B_eq, LB ,UB, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous):
nb_phases = len(G_next)
b_max = len(t)
## First Solver: IPOPT to get an initial guess
m_IPOPT = GEKKO(remote = True)
m_IPOPT.options.SOLVER = 3 #(IPOPT)
# Array Variable
rows = nb_phases + 3*b_max*(nb_phases+1)#48
x_initial = np.empty(rows,dtype=object)
for i in range(3*nb_phases*b_max+nb_phases+1):
x_initial[i] = m_IPOPT.Var(value = xinit[i], lb = LB[i], ub = UB[i])#, integer = False)
for i in range(3*nb_phases*b_max+nb_phases+1, (3*nb_phases+3)*b_max+nb_phases):
x_initial[i] = m_IPOPT.Var(value = xinit[i], lb = LB[i], ub = UB[i])#, integer = True)
# Constraints
m_IPOPT.axb(A,B,x_initial,etype = '<=',sparse=False)
m_IPOPT.axb(A_eq,B_eq,x_initial,etype = '=',sparse=False)
# Objective Function
f = objective_fun(x_initial, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous)
m_IPOPT.Obj(f)
#Solver
m_IPOPT.solve(disp = True)
####################################################################################################
## Second Solver: APOPT to solve MINLP
m_APOPT = GEKKO(remote = True)
m_APOPT.options.SOLVER = 1 #(APOPT)
x = np.empty(rows,dtype=object)
for i in range(3*nb_phases*b_max+nb_phases+1):
x[i] = m_APOPT.Var(value = x_initial[i], lb = LB[i], ub = UB[i], integer = False)
for i in range(3*nb_phases*b_max+nb_phases+1, (3*nb_phases+3)*b_max+nb_phases):
x[i] = m_APOPT.Var(value = x_initial[i], lb = LB[i], ub = UB[i], integer = True)
# Constraints
m_APOPT.axb(A,B,x,etype = '<=',sparse=False)
m_APOPT.axb(A_eq,B_eq,x,etype = '=',sparse=False)
# Objective Function
f = objective_fun(x, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous)
m_APOPT.Obj(f)
#Solver
m_APOPT.solve(disp = True)
return x
#Define Parameters
C = 120
T = 31
b_max = 2
G_base = [12,31,12,11,1,41]
G_max = [106, 106, 106, 106, 106, 106]
G_min = [7,3,7,10,0,7]
G_previous = [7.3333333, 34.16763, 7.1333333, 10.0, 2.2008602e-16, 47.365703]
Y = [2, 2, 3, 2, 2, 3]
jofbuses = [9,3]
k = [3,1,6,4,1,2,5,4,2,3]
l = [3,2,6,5,1,3,6,4,3,4]
nb_phases = 6
oa = 1.25
ob = [39,42]
t = [3600.2, 3603.5]
q = [0.038053758888888886, 0.215206065, 0.11325116416666667, 0.06299876472222223,0.02800455611111111,0.18878488361111112,0.2970903402777778, 0.01876728472222222, 0.2192723663888889, 0.06132227222222222]
qc = [0.04083333333333333, 0.2388888888888889, 0.10555555555555556, 0.0525, 0.030555555555555555, 0.20444444444444446,0.31083333333333335, 0.018333333333333333, 0.12777777777777777, 0.07138888888888889]
s = [1.0, 1.0, 1.0, 1.0, 0.5, 1.0, 1.0, 0.5, 0.5, 0.5]
nb_phases = len(G_base)
G_max = []
for i in range(nb_phases):
G_max.append(C - sum(Y[0:nb_phases]))
Optimise_G(t,ob, jofbuses, q, qc, s, oa, k, l, T, G_previous, C, Y, G_previous, G_max, G_min)
APOPT 求解器成功，returns 一个整数解。
APMonitor, Version 1.0.0
APMonitor Optimization Suite
----------------------------------------------------------------
--------- APM Model Size ------------
Each time step contains
Objects : 2
Constants : 0
Variables : 48
Intermediates: 0
Connections : 96
Equations : 1
Residuals : 1
Number of state variables: 48
Number of total equations: - 105
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : -57
* Warning: DOF <= 0
----------------------------------------------
Steady State Optimization with APOPT Solver
----------------------------------------------
Iter: 1 I: 0 Tm: 0.00 NLPi: 1 Dpth: 0 Lvs: 3 Obj: 1.64E+04 Gap: NaN
Iter: 2 I: -1 Tm: 0.00 NLPi: 2 Dpth: 1 Lvs: 2 Obj: 1.64E+04 Gap: NaN
Iter: 3 I: 0 Tm: 0.00 NLPi: 3 Dpth: 1 Lvs: 3 Obj: 1.81E+04 Gap: NaN
Iter: 4 I: -1 Tm: 0.01 NLPi: 6 Dpth: 1 Lvs: 2 Obj: 1.64E+04 Gap: NaN
--Integer Solution: 2.15E+04 Lowest Leaf: 1.81E+04 Gap: 1.68E-01
Iter: 5 I: 0 Tm: 0.00 NLPi: 4 Dpth: 2 Lvs: 1 Obj: 2.15E+04 Gap: 1.68E-01
Iter: 6 I: -1 Tm: 0.01 NLPi: 4 Dpth: 2 Lvs: 0 Obj: 1.81E+04 Gap: 1.68E-01
No additional trial points, returning the best integer solution
Successful solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 4.670000000623986E-002 sec
Objective : 21455.9882666580
Successful solution
---------------------------------------------------
如果您需要在 Linux 上使用 remote=False（本地解决方案），那么新的 Gekko 版本将提供新的可执行文件。存在一些编译器差异，创建本地 apm.exe 的 Windows FORTRAN 编译器似乎不如创建本地 apm 可执行文件的 Linux FORTRAN 编译器好。这些可执行文件位于 bin 文件夹中：Python 文件夹的 Lib\site-packages\gekko\bin。
我正在尝试解决 MINLP 问题，首先使用 IPOPT 求解器获得初始解，然后使用 APOPT 获得混合整数解。但是，调用 APOPT 求解器时出现以下错误：
Error: Exception: Access Violation At line 359 of file ./f90/cqp.f90 Traceback: not available, compile with -ftrace=frame or -ftrace=full Error: 'results.json' not found. Check above for additional error details Traceback (most recent call last): File "optimisation_stack.py", line 244, in Optimise_G(t,ob, jofbuses, q, qc, s, oa, k, l, T, G_previous, C, Y, G_previous, G_max, G_min) File "optimisation_stack.py", line 134, in Optimise_G sol = MINLP(xinit, A, B, A_eq, B_eq, LB ,UB, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_previous, C, Y, G_previous) File "optimisation_stack.py", line 215, in MINLP m_APOPT.solve(disp = False) File "C:\Users\Zineb\AppData\Local\Programs\Python\Python37\lib\site-packages\gekko\gekko.py", line 2227, in solve self.load_JSON() File "C:\Users\Zineb\AppData\Local\Programs\Python\Python37\lib\site-packages\gekko\gk_post_solve.py", line 13, in load_JSON f = open(os.path.join(self._path,'options.json')) FileNotFoundError: [Errno 2] No such file or directory: 'C:\Users\Zineb\AppData\Local\Temp\tmptdgafg1zgk_model1\options.json'
我的代码如下，我尽量简化了:
import numpy as np
from gekko import GEKKO
# Define matrices A,A_eq, and vectors b, b_eq for the optimization
def Optimise_G(t,ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous, G_max, G_min):
Mbig_1 = T*C
Mbig_2 = C
nb_phases = len(G_next)
b_max = len(t)
no_lanegroups = len(q)
A_eq = np.zeros(((nb_phases+1)*b_max + 1, (3*nb_phases+3)*b_max+nb_phases))
for i in range(nb_phases):
A_eq[0][i] = 1
#B_eq = np.zeros(((nb_phases+1)*b_max + 1, 1))
B_eq = np.zeros((nb_phases+1)*b_max + 1)
B_eq[0] = C - sum(Y[0:nb_phases])
counter_eq = 0
# G(i)=Ga(i,b)+Gb(i,b)+Gc(i,b)
for b in range(b_max):
for i in range(nb_phases):
counter_eq = counter_eq + 1
A_eq[counter_eq][i] = 1
A_eq[counter_eq][nb_phases*(b+1)+ i] = -1
A_eq[counter_eq][nb_phases*b_max + nb_phases*(b+1) + i] = -1
A_eq[counter_eq][2*nb_phases*b_max + nb_phases*(b+1) + i] = -1
# ya(b)+y(b)+y(c)=1
for b in range(b_max):
counter_eq = counter_eq + 1
A_eq[counter_eq][3*nb_phases*b_max + nb_phases + b] = 1
A_eq[counter_eq][(3*nb_phases+1)*b_max + nb_phases + b] = 1
A_eq[counter_eq][(3*nb_phases+2)*b_max + nb_phases + b] = 1
B_eq[counter_eq] = 1
A = np.zeros((no_lanegroups + (2*3*nb_phases+4)*b_max, (3*nb_phases+3)*b_max+nb_phases))
B = np.zeros(no_lanegroups + (2*3*nb_phases+4)*b_max)
counter = -1
# Sum Gi (i in Ij)>=Gj,min
for j in range(no_lanegroups):
counter = counter + 1
for i in range(k[j], l[j]+1):
A[counter][i-1] = -1
B[counter] = -C*qc[j]/s[j]
# ya(b)G_lb(i)<=Ga(i,b), yb(b)G_lb(i)<=Gb(i,b), yc(b)G_lb(i)<=Gc(i,b)
for b in range(b_max):
for i in range(nb_phases):
counter = counter + 1
A[counter][nb_phases*(b+1)+i] = -1
A[counter][3*nb_phases*b_max + nb_phases + b] = G_min[i]
B[counter] = 0
counter = counter + 1
A[counter][nb_phases*b_max + nb_phases*(b+1) + i] = -1
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = G_min[i]
B[counter] = 0
counter = counter + 1
A[counter][2*nb_phases*b_max + nb_phases*(b+1) +i] = -1
A[counter][(3*nb_phases+2)*b_max + nb_phases + b] = G_min[i]
B[counter] = 0
# ya(b)Gmax(i)>=Ga(i,b), yb(b)Gmax(i)>=Gb(i,b), yc(b)Gmax(i)>=Gc(i,b)
for b in range(b_max):
for i in range(nb_phases):
counter = counter + 1
A[counter][nb_phases*(b+1) +i] = 1
A[counter][3*nb_phases*b_max + nb_phases + b] = -G_max[i]
B[counter] = 0
counter = counter + 1
A[counter][nb_phases*b_max + nb_phases*(b+1) + i] = 1
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = -G_max[i]
B[counter] = 0
counter = counter + 1
A[counter][2*nb_phases*b_max + nb_phases*(b+1) +i] = 1
A[counter][(3*nb_phases+2)*b_max + nb_phases + b] = -G_max[i]
B[counter] = 0
# (1-yc(b))t(b)<=(T-1)C+sum(Gi(1:l(jofbuses(b))))+sum(Y(1:l(jofbuses(b))-1))
for b in range(b_max):
counter = counter + 1
A[counter][0:l[jofbuses[b]-1]] = -np.ones((1,l[jofbuses[b]-1]))
A[counter][(3*nb_phases+2)*b_max+nb_phases+b] = -t[b]
B[counter] = -t[b] + (T-1)*C + sum(Y[0:l[jofbuses[b]-1]-1])
# (T-1)C+sum(Gi(1:l(jofbuses(b))))+sum(Y(1:l(jofbuses(b))-1))<=yc(b)t(b)+(1-yc(b))Mbig_1
for b in range(b_max):
counter = counter + 1
A[counter][0:l[jofbuses[b]-1]] = np.ones((1,l[jofbuses[b]-1]))
A[counter][(3*nb_phases+2)*b_max+nb_phases+b] = -t[b] + Mbig_1
B[counter] = Mbig_1 - (T-1)*C - sum(Y[0:l[jofbuses[b]-1]-1])
# -Mbig_2(1-yb(b))<=db(b)=right-hand side of Equation (6)
for b in range(b_max):
counter = counter + 1
constant = q[jofbuses[b]-1]/s[jofbuses[b]-1]*(t[b] - (T-1)*C + sum(G_previous[l[jofbuses[b]-1]:nb_phases]) + sum(Y[l[jofbuses[b]-1] -1:nb_phases]))+ (T-1)*C + sum(Y[0:k[jofbuses[b]-1]-1]) - t[b]
A[counter][0:k[jofbuses[b]-1]-1] = -np.ones((1,k[jofbuses[b]-1]-1))
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = Mbig_2
B[counter] = constant + Mbig_2
# db(b)<=Mbig_2 yb(b)
for b in range(b_max):
counter = counter + 1
constant = q[jofbuses[b]-1]/s[jofbuses[b]-1]*(t[b] - (T-1)*C +sum(G_previous[l[jofbuses[b]-1]:nb_phases]) + sum(Y[l[jofbuses[b]-1] -1:nb_phases]))+ (T-1)*C + sum(Y[0:k[jofbuses[b]-1]-1]) - t[b]
A[counter][0:k[jofbuses[b]-1]-1] = np.ones((1,k[jofbuses[b]-1]-1))
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = -Mbig_2
B[counter] = -constant
#Lower Bound LB
LB_zeros = np.zeros(3*b_max*(nb_phases+1))
G_min = np.array(G_min)
LB = np.append(G_min, LB_zeros)
#Upper Bound UB
UB = np.ones(3*b_max)
G_max = np.array(G_max)
for i in range(3*b_max+1):
UB = np.concatenate((G_max,UB))
xinit = np.array([(a+b)/2 for a, b in zip(UB, LB)])
sol = MINLP(xinit, A, B, A_eq, B_eq, LB ,UB, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_previous, C, Y, G_previous)
objective函数：
def objective_fun(x, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous):
nb_phases = len(G_next)
b_max = len(t)
no_lanegroups = len(q)
obj = 0
obj_a = 0
obj_b = 0
G = x[0:nb_phases]
for j in range(no_lanegroups):
delay_a = 0.5*q[j]/(1-q[j]/s[j]) * (pow((sum(G_previous[l[j]:nb_phases]) + sum(G[0:k[j]-1]) + sum(Y[l[j]-1:nb_phases]) + sum(Y[0:k[j]-1])),2) + pow(sum(G[l[j]:nb_phases]) + sum(G_next[0:k[j]-1]) + sum(Y[l[j]-1:nb_phases]) + sum(Y[0:k[j]-1]),2))
obj = obj + oa*delay_a
obj_a = obj_a + oa*delay_a
for b in range(b_max):
delay_b1 = x[(3*nb_phases+1)*b_max + nb_phases + b]*(q[jofbuses[b]-1]/s[jofbuses[b]-1] * (t[b] - (T-1)*C + sum(G_previous[l[jofbuses[b]-1]:nb_phases]) + sum(Y[l[jofbuses[b]-1] -1:nb_phases])) + (T-1)*C - t[b] + sum(Y[0:k[jofbuses[b]-1]-1]))
delay_b2 = x[(3*nb_phases+2)*b_max + nb_phases + b-1]*(q[jofbuses[b]-1]/s[jofbuses[b]-1] * (t[b] - (T-1)*C - sum(Y[0:l[jofbuses[b]-1]-1])) + T*C + sum(G_next[0:k[jofbuses[b]-1]-1]) + sum(Y[0:k[jofbuses[b]-1]-1]) - t[b])
delay_b3 = sum(x[nb_phases*b_max + nb_phases*b:nb_phases*b_max + nb_phases*b+k[jofbuses[b]-1]-1]) - q[jofbuses[b]-1]/s[jofbuses[b]-1]*sum(x[2*nb_phases*b_max + nb_phases*b:2*nb_phases*b_max + nb_phases*b +l[jofbuses[b]-1]])
delay_b = delay_b1+delay_b2 +delay_b3
obj = obj + delay_b*ob[b]
obj_b = obj_b + delay_b*ob[b]
return obj
MINLP 求解器：
def MINLP(xinit, A, B, A_eq, B_eq, LB ,UB, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous):
nb_phases = len(G_next)
b_max = len(t)
## First Solver: IPOPT to get an initial guess
m_IPOPT = GEKKO(remote = False)
m_IPOPT.options.SOLVER = 3 #(IPOPT)
# Array Variable
rows = nb_phases + 3*b_max*(nb_phases+1)#48
x_initial = np.empty(rows,dtype=object)
for i in range(3*nb_phases*b_max+nb_phases+1):
x_initial[i] = m_IPOPT.Var(value = xinit[i], lb = LB[i], ub = UB[i])#, integer = False)
for i in range(3*nb_phases*b_max+nb_phases+1, (3*nb_phases+3)*b_max+nb_phases):
x_initial[i] = m_IPOPT.Var(value = xinit[i], lb = LB[i], ub = UB[i])#, integer = True)
# Constraints
m_IPOPT.axb(A,B,x_initial,etype = '<=',sparse=False)
m_IPOPT.axb(A_eq,B_eq,x_initial,etype = '=',sparse=False)
# Objective Function
f = objective_fun(x_initial, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous)
m_IPOPT.Obj(f)
#Solver
m_IPOPT.solve(disp = True)
####################################################################################################
## Second Solver: APOPT to solve MINLP
m_APOPT = GEKKO(remote = False)
m_APOPT.options.SOLVER = 1 #(APOPT)
x = np.empty(rows,dtype=object)
for i in range(3*nb_phases*b_max+nb_phases+1):
x[i] = m_APOPT.Var(value = x_initial[i], lb = LB[i], ub = UB[i], integer = False)
for i in range(3*nb_phases*b_max+nb_phases+1, (3*nb_phases+3)*b_max+nb_phases):
x[i] = m_APOPT.Var(value = x_initial[i], lb = LB[i], ub = UB[i], integer = True)
# Constraints
m_APOPT.axb(A,B,x,etype = '<=',sparse=False)
m_APOPT.axb(A_eq,B_eq,x,etype = '=',sparse=False)
# Objective Function
f = objective_fun(x, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous)
m_APOPT.Obj(f)
#Solver
m_APOPT.solve(disp = False)
return x
定义参数并调用Optimise_G:
#Define Parameters
C = 120
T = 31
b_max = 2
G_base = [12,31,12,11,1,41]
G_max = [106, 106, 106, 106, 106, 106]
G_min = [7,3,7,10,0,7]
G_previous = [7.3333333, 34.16763, 7.1333333, 10.0, 2.2008602e-16, 47.365703]
Y = [2, 2, 3, 2, 2, 3]
jofbuses = [9,3]
k = [3,1,6,4,1,2,5,4,2,3]
l = [3,2,6,5,1,3,6,4,3,4]
nb_phases = 6
oa = 1.25
ob = [39,42]
t = [3600.2, 3603.5]
q = [0.038053758888888886, 0.215206065, 0.11325116416666667, 0.06299876472222223,0.02800455611111111,0.18878488361111112,0.2970903402777778, 0.01876728472222222, 0.2192723663888889, 0.06132227222222222]
qc = [0.04083333333333333, 0.2388888888888889, 0.10555555555555556, 0.0525, 0.030555555555555555, 0.20444444444444446,0.31083333333333335, 0.018333333333333333, 0.12777777777777777, 0.07138888888888889]
s = [1.0, 1.0, 1.0, 1.0, 0.5, 1.0, 1.0, 0.5, 0.5, 0.5]
nb_phases = len(G_base)
G_max = []
for i in range(nb_phases):
G_max.append(C - sum(Y[0:nb_phases]))
Optimise_G(t,ob, jofbuses, q, qc, s, oa, k, l, T, G_previous, C, Y, G_previous, G_max, G_min)
有办法解决这个问题吗？ 非常感谢！
Windows 版本的 APOPT 求解器崩溃，无法找到解决方案。不过，在线Linux版本的APOPT是可以找到解决办法的。获取 Gekko (v1.0.0 pre-release) available on GitHub 的最新版本。当新版本发布时，这将在 pip install gekko --upgrade 中可用，但现在您需要将源代码复制到 Lib\site-packages\gekko\gekko.py。更新gekko后，切换到remote=True如下图。
import numpy as np
from gekko import GEKKO
# Define matrices A,A_eq, and vectors b, b_eq for the optimization
def Optimise_G(t,ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous, G_max, G_min):
Mbig_1 = T*C
Mbig_2 = C
nb_phases = len(G_next)
b_max = len(t)
no_lanegroups = len(q)
A_eq = np.zeros(((nb_phases+1)*b_max + 1, (3*nb_phases+3)*b_max+nb_phases))
for i in range(nb_phases):
A_eq[0][i] = 1
#B_eq = np.zeros(((nb_phases+1)*b_max + 1, 1))
B_eq = np.zeros((nb_phases+1)*b_max + 1)
B_eq[0] = C - sum(Y[0:nb_phases])
counter_eq = 0
# G(i)=Ga(i,b)+Gb(i,b)+Gc(i,b)
for b in range(b_max):
for i in range(nb_phases):
counter_eq = counter_eq + 1
A_eq[counter_eq][i] = 1
A_eq[counter_eq][nb_phases*(b+1)+ i] = -1
A_eq[counter_eq][nb_phases*b_max + nb_phases*(b+1) + i] = -1
A_eq[counter_eq][2*nb_phases*b_max + nb_phases*(b+1) + i] = -1
# ya(b)+y(b)+y(c)=1
for b in range(b_max):
counter_eq = counter_eq + 1
A_eq[counter_eq][3*nb_phases*b_max + nb_phases + b] = 1
A_eq[counter_eq][(3*nb_phases+1)*b_max + nb_phases + b] = 1
A_eq[counter_eq][(3*nb_phases+2)*b_max + nb_phases + b] = 1
B_eq[counter_eq] = 1
A = np.zeros((no_lanegroups + (2*3*nb_phases+4)*b_max, (3*nb_phases+3)*b_max+nb_phases))
B = np.zeros(no_lanegroups + (2*3*nb_phases+4)*b_max)
counter = -1
# Sum Gi (i in Ij)>=Gj,min
for j in range(no_lanegroups):
counter = counter + 1
for i in range(k[j], l[j]+1):
A[counter][i-1] = -1
B[counter] = -C*qc[j]/s[j]
# ya(b)G_lb(i)<=Ga(i,b), yb(b)G_lb(i)<=Gb(i,b), yc(b)G_lb(i)<=Gc(i,b)
for b in range(b_max):
for i in range(nb_phases):
counter = counter + 1
A[counter][nb_phases*(b+1)+i] = -1
A[counter][3*nb_phases*b_max + nb_phases + b] = G_min[i]
B[counter] = 0
counter = counter + 1
A[counter][nb_phases*b_max + nb_phases*(b+1) + i] = -1
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = G_min[i]
B[counter] = 0
counter = counter + 1
A[counter][2*nb_phases*b_max + nb_phases*(b+1) +i] = -1
A[counter][(3*nb_phases+2)*b_max + nb_phases + b] = G_min[i]
B[counter] = 0
# ya(b)Gmax(i)>=Ga(i,b), yb(b)Gmax(i)>=Gb(i,b), yc(b)Gmax(i)>=Gc(i,b)
for b in range(b_max):
for i in range(nb_phases):
counter = counter + 1
A[counter][nb_phases*(b+1) +i] = 1
A[counter][3*nb_phases*b_max + nb_phases + b] = -G_max[i]
B[counter] = 0
counter = counter + 1
A[counter][nb_phases*b_max + nb_phases*(b+1) + i] = 1
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = -G_max[i]
B[counter] = 0
counter = counter + 1
A[counter][2*nb_phases*b_max + nb_phases*(b+1) +i] = 1
A[counter][(3*nb_phases+2)*b_max + nb_phases + b] = -G_max[i]
B[counter] = 0
# (1-yc(b))t(b)<=(T-1)C+sum(Gi(1:l(jofbuses(b))))+sum(Y(1:l(jofbuses(b))-1))
for b in range(b_max):
counter = counter + 1
A[counter][0:l[jofbuses[b]-1]] = -np.ones((1,l[jofbuses[b]-1]))
A[counter][(3*nb_phases+2)*b_max+nb_phases+b] = -t[b]
B[counter] = -t[b] + (T-1)*C + sum(Y[0:l[jofbuses[b]-1]-1])
# (T-1)C+sum(Gi(1:l(jofbuses(b))))+sum(Y(1:l(jofbuses(b))-1))<=yc(b)t(b)+(1-yc(b))Mbig_1
for b in range(b_max):
counter = counter + 1
A[counter][0:l[jofbuses[b]-1]] = np.ones((1,l[jofbuses[b]-1]))
A[counter][(3*nb_phases+2)*b_max+nb_phases+b] = -t[b] + Mbig_1
B[counter] = Mbig_1 - (T-1)*C - sum(Y[0:l[jofbuses[b]-1]-1])
# -Mbig_2(1-yb(b))<=db(b)=right-hand side of Equation (6)
for b in range(b_max):
counter = counter + 1
constant = q[jofbuses[b]-1]/s[jofbuses[b]-1]*(t[b] - (T-1)*C + sum(G_previous[l[jofbuses[b]-1]:nb_phases]) + sum(Y[l[jofbuses[b]-1] -1:nb_phases]))+ (T-1)*C + sum(Y[0:k[jofbuses[b]-1]-1]) - t[b]
A[counter][0:k[jofbuses[b]-1]-1] = -np.ones((1,k[jofbuses[b]-1]-1))
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = Mbig_2
B[counter] = constant + Mbig_2
# db(b)<=Mbig_2 yb(b)
for b in range(b_max):
counter = counter + 1
constant = q[jofbuses[b]-1]/s[jofbuses[b]-1]*(t[b] - (T-1)*C +sum(G_previous[l[jofbuses[b]-1]:nb_phases]) + sum(Y[l[jofbuses[b]-1] -1:nb_phases]))+ (T-1)*C + sum(Y[0:k[jofbuses[b]-1]-1]) - t[b]
A[counter][0:k[jofbuses[b]-1]-1] = np.ones((1,k[jofbuses[b]-1]-1))
A[counter][(3*nb_phases+1)*b_max + nb_phases + b] = -Mbig_2
B[counter] = -constant
#Lower Bound LB
LB_zeros = np.zeros(3*b_max*(nb_phases+1))
G_min = np.array(G_min)
LB = np.append(G_min, LB_zeros)
#Upper Bound UB
UB = np.ones(3*b_max)
G_max = np.array(G_max)
for i in range(3*b_max+1):
UB = np.concatenate((G_max,UB))
xinit = np.array([(a+b)/2 for a, b in zip(UB, LB)])
sol = MINLP(xinit, A, B, A_eq, B_eq, LB ,UB, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_previous, C, Y, G_previous)
def objective_fun(x, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous):
nb_phases = len(G_next)
b_max = len(t)
no_lanegroups = len(q)
obj = 0
obj_a = 0
obj_b = 0
G = x[0:nb_phases]
for j in range(no_lanegroups):
delay_a = 0.5*q[j]/(1-q[j]/s[j]) * (pow((sum(G_previous[l[j]:nb_phases]) + sum(G[0:k[j]-1]) + sum(Y[l[j]-1:nb_phases]) + sum(Y[0:k[j]-1])),2) + pow(sum(G[l[j]:nb_phases]) + sum(G_next[0:k[j]-1]) + sum(Y[l[j]-1:nb_phases]) + sum(Y[0:k[j]-1]),2))
obj = obj + oa*delay_a
obj_a = obj_a + oa*delay_a
for b in range(b_max):
delay_b1 = x[(3*nb_phases+1)*b_max + nb_phases + b]*(q[jofbuses[b]-1]/s[jofbuses[b]-1] * (t[b] - (T-1)*C + sum(G_previous[l[jofbuses[b]-1]:nb_phases]) + sum(Y[l[jofbuses[b]-1] -1:nb_phases])) + (T-1)*C - t[b] + sum(Y[0:k[jofbuses[b]-1]-1]))
delay_b2 = x[(3*nb_phases+2)*b_max + nb_phases + b-1]*(q[jofbuses[b]-1]/s[jofbuses[b]-1] * (t[b] - (T-1)*C - sum(Y[0:l[jofbuses[b]-1]-1])) + T*C + sum(G_next[0:k[jofbuses[b]-1]-1]) + sum(Y[0:k[jofbuses[b]-1]-1]) - t[b])
delay_b3 = sum(x[nb_phases*b_max + nb_phases*b:nb_phases*b_max + nb_phases*b+k[jofbuses[b]-1]-1]) - q[jofbuses[b]-1]/s[jofbuses[b]-1]*sum(x[2*nb_phases*b_max + nb_phases*b:2*nb_phases*b_max + nb_phases*b +l[jofbuses[b]-1]])
delay_b = delay_b1+delay_b2 +delay_b3
obj = obj + delay_b*ob[b]
obj_b = obj_b + delay_b*ob[b]
return obj
def MINLP(xinit, A, B, A_eq, B_eq, LB ,UB, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous):
nb_phases = len(G_next)
b_max = len(t)
## First Solver: IPOPT to get an initial guess
m_IPOPT = GEKKO(remote = True)
m_IPOPT.options.SOLVER = 3 #(IPOPT)
# Array Variable
rows = nb_phases + 3*b_max*(nb_phases+1)#48
x_initial = np.empty(rows,dtype=object)
for i in range(3*nb_phases*b_max+nb_phases+1):
x_initial[i] = m_IPOPT.Var(value = xinit[i], lb = LB[i], ub = UB[i])#, integer = False)
for i in range(3*nb_phases*b_max+nb_phases+1, (3*nb_phases+3)*b_max+nb_phases):
x_initial[i] = m_IPOPT.Var(value = xinit[i], lb = LB[i], ub = UB[i])#, integer = True)
# Constraints
m_IPOPT.axb(A,B,x_initial,etype = '<=',sparse=False)
m_IPOPT.axb(A_eq,B_eq,x_initial,etype = '=',sparse=False)
# Objective Function
f = objective_fun(x_initial, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous)
m_IPOPT.Obj(f)
#Solver
m_IPOPT.solve(disp = True)
####################################################################################################
## Second Solver: APOPT to solve MINLP
m_APOPT = GEKKO(remote = True)
m_APOPT.options.SOLVER = 1 #(APOPT)
x = np.empty(rows,dtype=object)
for i in range(3*nb_phases*b_max+nb_phases+1):
x[i] = m_APOPT.Var(value = x_initial[i], lb = LB[i], ub = UB[i], integer = False)
for i in range(3*nb_phases*b_max+nb_phases+1, (3*nb_phases+3)*b_max+nb_phases):
x[i] = m_APOPT.Var(value = x_initial[i], lb = LB[i], ub = UB[i], integer = True)
# Constraints
m_APOPT.axb(A,B,x,etype = '<=',sparse=False)
m_APOPT.axb(A_eq,B_eq,x,etype = '=',sparse=False)
# Objective Function
f = objective_fun(x, t, ob, jofbuses, q, qc, s, oa, k, l, T, G_next, C, Y, G_previous)
m_APOPT.Obj(f)
#Solver
m_APOPT.solve(disp = True)
return x
#Define Parameters
C = 120
T = 31
b_max = 2
G_base = [12,31,12,11,1,41]
G_max = [106, 106, 106, 106, 106, 106]
G_min = [7,3,7,10,0,7]
G_previous = [7.3333333, 34.16763, 7.1333333, 10.0, 2.2008602e-16, 47.365703]
Y = [2, 2, 3, 2, 2, 3]
jofbuses = [9,3]
k = [3,1,6,4,1,2,5,4,2,3]
l = [3,2,6,5,1,3,6,4,3,4]
nb_phases = 6
oa = 1.25
ob = [39,42]
t = [3600.2, 3603.5]
q = [0.038053758888888886, 0.215206065, 0.11325116416666667, 0.06299876472222223,0.02800455611111111,0.18878488361111112,0.2970903402777778, 0.01876728472222222, 0.2192723663888889, 0.06132227222222222]
qc = [0.04083333333333333, 0.2388888888888889, 0.10555555555555556, 0.0525, 0.030555555555555555, 0.20444444444444446,0.31083333333333335, 0.018333333333333333, 0.12777777777777777, 0.07138888888888889]
s = [1.0, 1.0, 1.0, 1.0, 0.5, 1.0, 1.0, 0.5, 0.5, 0.5]
nb_phases = len(G_base)
G_max = []
for i in range(nb_phases):
G_max.append(C - sum(Y[0:nb_phases]))
Optimise_G(t,ob, jofbuses, q, qc, s, oa, k, l, T, G_previous, C, Y, G_previous, G_max, G_min)
APOPT 求解器成功，returns 一个整数解。
APMonitor, Version 1.0.0
APMonitor Optimization Suite
----------------------------------------------------------------
--------- APM Model Size ------------
Each time step contains
Objects : 2
Constants : 0
Variables : 48
Intermediates: 0
Connections : 96
Equations : 1
Residuals : 1
Number of state variables: 48
Number of total equations: - 105
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : -57
* Warning: DOF <= 0
----------------------------------------------
Steady State Optimization with APOPT Solver
----------------------------------------------
Iter: 1 I: 0 Tm: 0.00 NLPi: 1 Dpth: 0 Lvs: 3 Obj: 1.64E+04 Gap: NaN
Iter: 2 I: -1 Tm: 0.00 NLPi: 2 Dpth: 1 Lvs: 2 Obj: 1.64E+04 Gap: NaN
Iter: 3 I: 0 Tm: 0.00 NLPi: 3 Dpth: 1 Lvs: 3 Obj: 1.81E+04 Gap: NaN
Iter: 4 I: -1 Tm: 0.01 NLPi: 6 Dpth: 1 Lvs: 2 Obj: 1.64E+04 Gap: NaN
--Integer Solution: 2.15E+04 Lowest Leaf: 1.81E+04 Gap: 1.68E-01
Iter: 5 I: 0 Tm: 0.00 NLPi: 4 Dpth: 2 Lvs: 1 Obj: 2.15E+04 Gap: 1.68E-01
Iter: 6 I: -1 Tm: 0.01 NLPi: 4 Dpth: 2 Lvs: 0 Obj: 1.81E+04 Gap: 1.68E-01
No additional trial points, returning the best integer solution
Successful solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 4.670000000623986E-002 sec
Objective : 21455.9882666580
Successful solution
---------------------------------------------------
如果您需要在 Linux 上使用 remote=False（本地解决方案），那么新的 Gekko 版本将提供新的可执行文件。存在一些编译器差异，创建本地 apm.exe 的 Windows FORTRAN 编译器似乎不如创建本地 apm 可执行文件的 Linux FORTRAN 编译器好。这些可执行文件位于 bin 文件夹中：Python 文件夹的 Lib\site-packages\gekko\bin。