The code runs for a particularly long time
Answered
My optimization problem is as follows, and optimization variables are underlined in red.
The corresponding code is as follows，
from gurobipy import *
import gurobipy as gp
from numpy import *
import numpy as np
import math
import random
#define global variables
global K
K = 10
global M
M = 5
global ATtheta1
ATtheta1=0.42
global ATtheta2
ATtheta2=2.00
global H
H=100
global Pmax
Pmax=10**(-2)
global N0w
#N0w=10**(-11) #to avoid large coefficients numerical issue
N0w=0
global beta0
beta0=10**(-4)
global sumcombination
sumcombination=int(math.factorial(K)/math.factorial(M)/math.factorial(K-M))
#The Positions of K users are given randomly
global upos
upos=[[103.1791349465347, 245.55998960612703], [255.48098415065635, 320.2532665566555], [64.55202835511899, 144.57501761863028], [379.2670054991101, 260.5062095134453], [28.71787081756203, 327.8530626055657], [14.324474423137357, 288.27329202612844], [206.68794921797164, 10.573469015806047], [60.666225853541754, 279.34522636111916], [312.530530067231, 93.03878739105653], [337.8884265853058, 197.4821008557357]]
#define distance
def DS(i,j):
return math.sqrt((upos[i][0]-upos[j][0])**2+(upos[i][1]-upos[j][1])**2)
#The maximum value is obtained in order to normalize the distance
global DSset
DSset= []
global DSmax
for i in range(0,K):
for j in range(0,K):
DSset.append(DS(i,j))
DSmax=max(DSset)
#define function
def cov(i,j):
return math.exp(-((DS(i,j)/DSmax/ATtheta1)**ATtheta2))
global deltaz
sum0=0
for k in range(0,K):
for k1 in range(0,K):
sum0=sum0+cov(k,k1)
deltaz=sum0/K/K
#initialized variables
b0=math.sqrt(Pmax)*np.array(ones((K)))
#Obtain the maximum distance to tighten the range of variable a
global DSset_a
DSset_a= []
global DSmax_a
for i in range(0,K):
DSset_a.append((upos[i][0]-0)**2+(upos[i][1]-0)**2)
DSset_a.append((upos[i][0]-400)**2+(upos[i][1]-0)**2)
DSset_a.append((upos[i][0]-0)**2+(upos[i][1]-400)**2)
DSset_a.append((upos[i][0]-400)**2+(upos[i][1]-400)**2)
DSmax_a=max(DSset_a)
#create model
model1 = gp.Model("MyModel1")
#declare variables
q=model1.addVars(2, lb=0, ub=400, vtype=GRB.CONTINUOUS, name="q")
a=model1.addVars(K, lb=0, ub=DSmax_a, name="a")
v=model1.addVars(K, lb=beta0/(H**2+DSmax_a), ub=beta0/H/H, name="v")
r=model1.addVars(K, lb=math.sqrt(beta0/(H**2+DSmax_a)), ub=math.sqrt(beta0/H/H), name="r")
t=model1.addVars(K, K, lb=beta0/(H**2+DSmax_a), ub=beta0/H/H, name="t")
nu=model1.addVars(sumcombination, lb=0, ub=deltaz, name="nu")
#set objective(The objective function is to choose M numbers from 1 to K, and for each combination, calculate the sum7, and finally obtain the sum value of K!/M!/(K-M)! combinations, which is shown as sumvalue. For how to get the combination, we use the depth-first search approach. In bold below are the optimization variables.
sumvalue=0
for i in range(0, sumcombination):
sumvalue=sumvalue+deltaz-nu[i]
obj=math.factorial(M)*math.factorial(K-M)/math.factorial(K)*sumvalue
model1.setObjective(obj, GRB.MINIMIZE)
model1.addConstrs((a[k]==(q[0]-upos[k][0])**2+(q[1]-upos[k][1])**2 for k in range(0,K)), name="c1")
model1.addConstrs((v[k]*(H**2+a[k])/beta0 == 1 for k in range(0,K)), name="c2")
model1.addConstrs((v[k] == r[k]*r[k] for k in range(0,K)), name="c3")
for k in range(0,K):
for k1 in range(0,K):
model1.addConstr(t[k,k1] == r[k]*r[k1], name="c4")
combination=[]
index=0
def dfs1(last_idx, cnt):
global index
if cnt == M:
print(combination)
sum1=0
for k in range(0,K):
for km in range(0,M):
sum1=sum1+(b0[combination[km]]*r[combination[km]]*cov(k,combination[km]))
sum2=0
for km in range(0,M):
for km1 in range(0,M):
sum2=sum2+(b0[combination[km]]*b0[combination[km1]]*t[combination[km],combination[km1]]*cov(combination[km],combination[km1]))
model1.addConstr(sum1**2 == nu[index]*K*K*(sum2+N0w), name="c5")
index=index+1
return
for idx in range(last_idx+1,K):
combination.append(idx)
dfs1(idx, cnt+1)
combination.pop()
dfs1(-1,0)
model1.params.NonConvex=2
model1.params.NumericFocus=2
model1.params.FeasibilityTol=1e-9
model1.params.OptimalityTol=1e-9
model1.params.ScaleFlag=2
model1.params.Aggregate=0
model1.params.MIPFocus=1
model1.optimize()
model1.printQuality()
model1.feasRelaxS()
print("Obj: ", model1.objVal)
print("q: ", q)
However,
Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (win64)
CPU model: 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 0 rows, 384 columns and 0 nonzeros
Model fingerprint: 0x7dcd91d6
Model has 382 quadratic constraints
Coefficient statistics:
Matrix range [0e+00, 0e+00]
QMatrix range [3e-03, 1e+04]
QLMatrix range [1e+00, 1e+08]
Objective range [4e-03, 4e-03]
Bounds range [4e-10, 2e+05]
RHS range [0e+00, 0e+00]
QRHS range [1e+00, 2e+05]
Continuous model is non-convex -- solving as a MIP
Presolve time: 0.00s
Presolved: 25928 rows, 6686 columns, 61820 nonzeros
Presolved model has 6422 bilinear constraint(s)
Variable types: 6686 continuous, 0 integer (0 binary)
Root relaxation: objective 1.331575e-02, 15360 iterations, 1.05 seconds (2.72 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 0.01332 0 10 - 0.01332 - - 1s
0 0 0.01332 0 10 - 0.01332 - - 1s
0 2 0.01332 0 10 - 0.01332 - - 6s
7 8 0.01332 3 11 - 0.01332 - 7251 50s
11 12 0.01339 4 11 - 0.01332 - 4695 71s
19 20 0.01360 5 15 - 0.01332 - 7102 107s
23 25 0.01337 5 11 - 0.01332 - 8756 114s
30 28 infeasible 6 - 0.01332 - 7198 120s
35 32 infeasible 7 - 0.01332 - 6236 145s
41 43 0.01445 8 12 - 0.01332 - 6454 189s
52 51 0.01435 10 11 - 0.01332 - 6909 201s
92 89 0.01490 14 11 - 0.01332 - 4293 226s
154 139 0.01558 22 10 - 0.01332 - 3135 246s
206 171 0.01355 9 11 - 0.01332 - 2717 250s
503 386 0.01915 28 10 - 0.01332 - 1196 255s
647 539 0.01931 34 11 - 0.01332 - 944 276s
* 741 375 90 0.0255348 0.01332 47.8% 838 276s
746 377 0.01932 50 11 0.02553 0.01332 47.8% 900 325s
* 821 310 77 0.0193228 0.01332 31.0% 933 328s
845 315 0.01332 17 10 0.01932 0.01332 31.0% 921 335s
880 393 0.01383 31 15 0.01932 0.01332 31.0% 916 342s
998 485 0.01332 10 11 0.01932 0.01332 31.0% 841 364s
* 1152 399 112 0.0187885 0.01332 29.1% 781 364s
1193 397 0.01332 21 0 0.01879 0.01332 29.1% 761 366s
1197 400 0.01411 22 0 0.01879 0.01411 24.9% 758 371s
1200 407 0.01507 16 15 0.01879 0.01454 22.6% 769 410s
1204 407 infeasible 17 0.01879 0.01477 21.4% 817 415s
1212 411 infeasible 18 0.01879 0.01494 20.5% 831 433s
1220 413 infeasible 19 0.01879 0.01497 20.3% 852 435s
1232 414 0.01519 20 15 0.01879 0.01497 20.3% 861 450s
1238 421 0.01523 21 15 0.01879 0.01497 20.3% 894 455s
1374 433 0.01641 27 13 0.01879 0.01497 20.3% 837 483s
1399 431 infeasible 29 0.01879 0.01497 20.3% 849 489s
1417 427 infeasible 29 0.01879 0.01497 20.3% 846 513s
1437 425 infeasible 22 0.01879 0.01497 20.3% 877 518s
1479 416 infeasible 28 0.01879 0.01497 20.3% 874 527s
1506 415 0.01573 25 10 0.01879 0.01497 20.3% 880 533s
1538 398 0.01862 27 15 0.01879 0.01497 20.3% 876 551s
1582 383 cutoff 27 0.01879 0.01500 20.2% 913 564s
1600 373 0.01641 26 11 0.01879 0.01500 20.2% 933 588s
1664 348 cutoff 24 0.01879 0.01500 20.2% 924 590s
......
325432 183862 0.01777 130 15 0.01799 0.01776 1.28% 155 11784s
325904 184452 cutoff 123 0.01799 0.01776 1.28% 155 11810s
326696 185019 0.01777 79 15 0.01799 0.01776 1.28% 154 11836s
327447 185364 0.01777 111 15 0.01799 0.01776 1.28% 154 11858s
328024 185763 0.01777 86 15 0.01799 0.01776 1.28% 154 11883s
328611 186067 cutoff 109 0.01799 0.01776 1.28% 154 11900s
329141 186628 0.01777 88 15 0.01799 0.01776 1.28% 154 11922s
329906 187006 infeasible 79 0.01799 0.01776 1.28% 154 11943s
330462 187326 0.01777 100 15 0.01799 0.01776 1.28% 154 11970s
The running time is very long and there is no output result. Thus，we interrupted the code and the output is as follows.
445484 257485 0.01777 64 14 0.01799 0.01777 1.27% 142 16604s 446103 257785 0.01777 161 14 0.01799 0.01777 1.27% 142 16627s Interrupt request received 446665 258081 infeasible 77 0.01799 0.01777 1.27% 142 16644s Explored 447119 nodes (63297580 simplex iterations) in 16644.65 seconds (26060.42 work units) Thread count was 8 (of 8 available processors) Solution count 10: 0.0179948 0.0179955 0.0179955 ... 0.0181281 Solve interrupted Warning: max constraint violation (4.4369e-09) exceeds tolerance Best objective 1.799483059283e-02, best bound 1.776545063259e-02, gap 1.2747% Solution quality statistics for model MyModel1 : Maximum violation: Bound : 0.00000000e+00 Constraint : 4.43693149e-09 (c5) Integrality : 0.00000000e+00
Moreover, we add
model1.feasRelaxS(0, False, False, True) model1.optimize()
or
model1.feasRelaxS(0, False, True, True) model1.optimize()
The output results are all 0, which is incorrect. Besides, I tried some methods. On the one hand, since the lower bound of variables v and t can reach 1e-10, I scale the variables as
v=model1.addVars(K, lb=(10**3)*beta0/(H**2+DSmax_a), ub=(10**3)*beta0/H/H, name="v") t=model1.addVars(K, K, lb=(10**3)*beta0/(H**2+DSmax_a), ub=(10**3)*beta0/H/H, name="t")
and then the corresponding constraints are changed to
model1.addConstrs((v[k]*(H**2+a[k])/beta0 == 1*(10**(3)) for k in range(0,K)), name="c2") model1.addConstrs(v[k] == (10**3)*r[k]*r[k] for k in range(0,K)) for k in range(0,K): for k1 in range(0,K): model1.addConstr(t[k,k1] == (10**3)*r[k]*r[k1], name="c3")
On the other hand, in constraint (e), since the sum1 and sum2 have orders of magnitude 1e-5 and 1e-10, respectively, I multiply both sides of the equation by 1e3. However, the simulation time is still very long and the results cannot be obtained. Thus, how to improve the performance of the code? Looking forward to the reply. Many thanks！
Yali
0
-
-
Gurobi Staff
Hi Yali,
The violations you obtain (4e-9) are very low. If those violations are still too large for your application, then you need to scale the ranges of your coefficients significantly. We usually recommend setting optimality and feasibility tolerances smaller than the smallest coefficients, but this is not possible in your case since 1e-9 is already the lowest tolerance value. Maybe you can find more hints in our guidelines for numerical issues.
It is not surprising that solving the feasibility relaxation leads to values of 0. The feasibility relaxation minimizes the violations in the model and 0 means that there exists a feasible solution subject to tolerances. Provided that your violations are in the range of 4e-9, this outcome is expected.
There are several improvements for your type of model in version 11 which will be released in about 1 month. Please try it out!
Best regards,
Mario0 -
Please sign in to leave a comment.
Comments
1 comment