我如何用Gekko来表述我的优化问题?

Question

我想表达目标函数(最小化问题):sum[sum[Ri*{Pi²+ (Qi - Qcj*Xij)²}for j in range(Nc)] for I in range(N)]， p和Q是常数，Qc是建议解的列表，X是我们的决策变量(二元变量)。我想求出使目标函数最小的向量X。下面是我的尝试:

from gekko import GEKKO
import numpy as np
P=[13.10511598922975,11.2611396806742,10.103920431906348,8.199519500182628,6.411296067052755,4.753519719147589,3.8977762462825973,2.6593092284662734,1.6399999999854893]
Q=[5.06643685386732,4.4344047044589585,3.8082608015186405,3.2626022579039584,1.2568869621197523,0.6152693459109657,0.46237064874523776,0.35226399840832523,0.20000000001140983]
R=[0.1233, 0.014, 0.7463, 0.6984, 1.9831, 0.9053, 2.0552, 4.7953, 5.3434]
Qc=[150, 300, 450, 600,750, 900,1050, 1200,1350,1500,1650,1800,1950,2100,2250,2400,2550,2700,2850,3000,3150,3300,3450,3600,3750,3900,4050]
N=len(Q)
Nc=len(Qc)
m = GEKKO(remote=False)
X = m.Array(m.Var,(N,Nc),integer=True,lb=0,ub=1,value=0)
#convirtir P et Q en KW
for i in range(N):
Q[i]=Q[i]*1000
P[i]=P[i]*1000
#constrainte ## one per line
for i in range(N):
m.Equation(m.sum([X[i][j]for j in range(Nc)])<=1)    
b=m.sum([m.sum([R[i]*((P[i]**2)+((Q[i])-Qc[j]*X[i][j])**2) for j in range(Nc)]) for i in range(N)])
m.Minimize(b)

我尝试了3种方法:方法1:

m.options.SOLVER = 1
m.solve()

方法2:

bv = np.array([[0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1],
[0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0, 0, 0, 0, 0, 0, 1, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0, 0, 0, 0, 1, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0, 0, 0, 0, 1, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]])
for i in range(N):
for j in range(Nc):
X[i,j].value = bv[i,j]
m.options.SOLVER = 1
m.solve()

方法3:

m.options.SOLVER = 3
m.solve(debug=0, disp=True)
m.options.SOLVER = 1
m.solve(debug=0, disp=True)

下面是我的尝试:我尝试了3种方法:方法1:方法2方法3:这三种方法都不能给我最优的解决方案。

Answer 1

使用求解器选项来获得更好的解决方案，当gap_tol满足1e-3时不终止(默认)。gap_tol是一种早期终止准则，有助于更快地获得MINLP解，但不太优解。将gap_tol设置为零，并将minlp_max_iter_with_int_sol设置为大数，将遍历所有剩余的可能解。计算时间增加，所以我建议使用较小的gap_tol，如1e-5和minlp_max_iter_with_int_sol 2000。

m.solver_options = ['minlp_gap_tol 1e-5',
'minlp_maximum_iterations 10000',
'minlp_max_iter_with_int_sol 2000']
m.options.SOLVER=1
m.solve(disp=True)

给出了一个目标为9.535e9的解决方案。

---------------------------------------------------
Solver         :  APOPT (v1.0)
Solution time  :    170.673799999990      sec
Objective      :    9535331689.96189     
Successful solution
---------------------------------------------------

初始猜测固定的目标是9.541e9。

bv = np.array([[0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1],
[0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0, 0, 0, 0, 0, 0, 1, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0, 0, 0, 0, 1, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0, 0, 0, 0, 1, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]])
for i in range(N):
for j in range(Nc):
X[i,j].value = bv[i,j]
m.Equation(X[i,j]==bv[i,j])

---------------------------------------------------
Solver         :  APOPT (v1.0)
Solution time  :   2.180000001681037E-002 sec
Objective      :    9540896947.56266     
Successful solution
---------------------------------------------------

还有一些建议对解决方案的准确性没有帮助，但可能是将来需要考虑的事情:

1. 用这种方法表示目标函数可以提高求解速度。

b=[sum([R[i]*((P[i]**2)+((Q[i])-Qc[j]*X[i][j])**2)
for j in range(Nc)]) for i in range(N)]
[m.Minimize(bi) for bi in b]

2. 目标函数是相当高的>1e9在解决方案。您可以考虑通过删除这些线来增加求解器公差，并将问题保持在MW与kW之间。

#convirtir P et Q en KW
#for i in range(N):
#    Q[i]=Q[i]*1000
#    P[i]=P[i]*1000