• Gurobi Staff

Hi Yali,

You code snippet is not executable. The are random $$\texttt{return()}$$ statements with the function definition being commented-out.

As the warning message indicates, the linear part of your quadratic constraints has a very wide range: it includes coefficients as small as 1e-9 and as big as 1e8.

• Does a coefficient as small as 1e-9 have any meaning in your application? If this is just noise and does not have any meaning, you can consider zeroing these tiny coefficients.
• The default feasibility and optimality tolerance values are 1e-6 and can be tightened as small as 1e-9. It is a good practice to set these tolerances to be at least one order of magnitude smaller than the smallest meaningful coefficient in your model.
• We highly recommended defining tight lower and upper bounds for all the variables participating in non-linear terms.
• You can consider experimenting with parameters such as ScaleFlag. Setting this parameter to value 2 is typically helpful for models with numerical challenges.
• Why did you turn the presolve off? Some presolve procedures, in particular, Aggregate can worsen the numerical challenges. You can experiment with Aggregate=0 if the presolve is indeed troublesome for your model.
• You can consider experimenting with higher values of NumericFocus such as 2 and 3.
• You can use the method Model.printQuality() available in Gurobi Python API to find out the constraint name which has the maximum constraint violation.

The suggestions above would not address the root cause of the problem. It might help to alleviate the numerical challenges of the model. The best approach would be to consider scaling your variables to get the model coefficients to be more reasonable.

Best regards,

Maliheh

Hi, Maliheh

We use Gurobi to sove the following problem. First, let me introduce the details. We choose M numbers (denoted as the set ) from 1 to K (denoted as the set ), and for each combination, calculate the MSE, and finally obtain the sum value of K!/M!/(K-M)! combinations.

The optimization problem are shown as follows. Optimization variables are underlined in red.

To facilitate the solution, we introduce auxiliary variables. Then, the optimizaion problem is reformulatd as folllows. Optimization variables are underlined in blue.

According to the your comments, we return the noise power spectral density to 0, rather than a very small number. We add "model.params.NumericFocus=2, model.params.FeasibilityTol=10**(-9), model.params.OptimalityTol=10**(-9), model.params.ScaleFlag=2, model.params.Aggregate=0". We have modified the upper and lower bounds of variables, but for variables in constraint (f), we cannot provide tighter upper and lower bounds. Currently, our code and existing issues are as follows.

from gurobipy import *import gurobipy as gpfrom numpy import *import numpy as npimport mathimport random#define global variablesglobal KK = 10global MM = 5global ATtheta1ATtheta1=0.42global ATtheta2ATtheta2=2.00global CPtheta1CPtheta1=0.11global CPtheta2CPtheta2=0.68global CStheta1CStheta1=0.05global CStheta2CStheta2=0.62global HH=100global PmaxPmax=10**(-2)  global N0w#N0w=10**(-11)N0w=0global beta0beta0=10**(-4)  global thresholdthreshold=10**(-3)global sumcombinationsumcombination=int(math.factorial(K)/math.factorial(M)/math.factorial(K-M))#The Positions of K users are given randomlyglobal uposupos = []for i in range(0,K):    upos.append([400*random.random(),400*random.random()])#define distancedef 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 distanceglobal DSsetDSset= []global DSmaxfor i in range(0,K):    for j in range(0,K):            DSset.append(DS(i,j))DSmax=max(DSset)  #define functiondef cov(i,j):    return math.exp(-((DS(i,j)/DSmax/ATtheta1)**ATtheta2))##global deltazsum0=0for k in range(0,K):    for k1 in range(0,K):            sum0=sum0+cov(k,k1)deltaz=sum0/K/K##initialized variablesb0=math.sqrt(Pmax)*np.array(ones((K)))q0=np.array([250,250])#create modelmodel = gp.Model("MyModel2")#declare variablesq=model.addVars(2, lb=0, ub=400, vtype=GRB.CONTINUOUS, name="q")a=model.addVars(K, lb=0, ub=(400*math.sqrt(2))**2, name="a")v=model.addVars(K, lb=beta0/(H**2+(400*math.sqrt(2))**2), ub=beta0/H/H, name="v")r=model.addVars(K, lb=math.sqrt(beta0/(H**2+(400*math.sqrt(2))**2)), ub=math.sqrt(beta0/H/H), name="r")t=model.addVars(K, K, lb=beta0/(H**2+(400*math.sqrt(2))**2), ub=beta0/H/H, name="t")#varpi=model.addVars(sumcombination, lb=0, name="varpi")nu=model.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=0for i in range(0, sumcombination):      sumvalue=sumvalue+deltaz-nu[i]obj=math.factorial(M)*math.factorial(K-M)/math.factorial(K)*sumvaluemodel.setObjective(obj, GRB.MINIMIZE)model.addConstrs(a[k]==(q[0]-upos[k][0])**2+(q[1]-upos[k][1])**2 for k in range(0,K))model.addConstrs(v[k]*(H**2+a[k])/beta0 == 1 for k in range(0,K))#for k in range(0,K):    #model.addGenConstrPow(v[k], r[k], 0.5, "FuncPieces=-2 FuncPieceError=0.001")model.addConstrs(v[k] == (r[k])**2 for k in range(0,K))for k in range(0,K):    for k1 in range(0,K):        model.addConstr(t[k,k1] == r[k]*r[k1])combination=[]index=0def 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]))          #model.addConstr(varpi[index] == sum1)        #model.addConstr((varpi[index])**2 == nu[index]*K*K*(sum2+N0w))        model.addConstr(sum1**2 == nu[index]*K*K*(sum2+N0w))        #        index=index+1        return    for idx in range(last_idx+1,K):        combination.append(idx)        dfs1(idx, cnt+1)        combination.pop()dfs1(-1,0)model.params.NonConvex=2#model.params.Presolve=0model.params.NumericFocus=2model.params.FeasibilityTol=10**(-9)model.params.OptimalityTol=10**(-9)model.params.ScaleFlag=2model.params.Aggregate=0#model.optimize()print("Obj: ", model.objVal)print("q: ", q)opt=model.objValmodel.printQuality()
Changed value of parameter NonConvex to 2   Prev: -1  Min: -1  Max: 2  Default: -1Changed value of parameter NumericFocus to 2   Prev: 0  Min: 0  Max: 3  Default: 0Changed value of parameter FeasibilityTol to 1e-09   Prev: 1e-06  Min: 1e-09  Max: 0.01  Default: 1e-06Changed value of parameter OptimalityTol to 1e-09   Prev: 1e-06  Min: 1e-09  Max: 0.01  Default: 1e-06Changed value of parameter ScaleFlag to 2   Prev: -1  Min: -1  Max: 3  Default: -1Changed value of parameter Aggregate to 0   Prev: 1  Min: 0  Max: 1  Default: 1Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)Thread count: 4 physical cores, 8 logical processors, using up to 8 threadsOptimize a model with 0 rows, 384 columns and 0 nonzerosModel fingerprint: 0x7abb6cb6Model has 382 quadratic constraintsCoefficient 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     [3e-10, 3e+05]  RHS range        [0e+00, 0e+00]  QRHS range       [1e+00, 3e+05]Continuous model is non-convex -- solving as a MIP.Presolve time: 0.00sPresolved: 25928 rows, 6686 columns, 61820 nonzerosPresolved model has 6422 bilinear constraint(s)Variable types: 6686 continuous, 0 integer (0 binary)Root relaxation: objective 1.297265e-02, 11453 iterations, 0.77 seconds    Nodes    |    Current Node    |     Objective Bounds      |     Work Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time     0     0    0.01297    0   10          -    0.01297      -     -    0s     0     0    0.01298    0   10          -    0.01298      -     -    0s     0     2    0.01298    0   10          -    0.01298      -     -    3s     1     4    0.01314    1   12          -    0.01298      - 42634   20s    15    18    0.01314    4   12          -    0.01314      -  3967   25s    71    70    0.01393   10   12          -    0.01314      -  1097   36s    79    76 infeasible   12               -    0.01314      -  1429   55s    93    83    0.01524   13   10          -    0.01314      -  1916   66s   112    94 infeasible   14               -    0.01314      -  1701   81s   135   109    0.01524   16   11          -    0.01314      -  1641   87s   154   125    0.01524   17   11          -    0.01314      -  1591   98s   185   146 infeasible   20               -    0.01314      -  1475  108s   201   157    0.01533   21   12          -    0.01314      -  1499  111s   226   175 infeasible   24               -    0.01314      -  1351  119s   254   185    0.01342    6   10          -    0.01314      -  1299  121s   272   201 infeasible    7               -    0.01314      -  1230  126s   294   216    0.01314    8  328          -    0.01314      -  1230  130s   313   232    0.01314    9  329          -    0.01314      -  1344  151s   341   249    0.01474   10   12          -    0.01314      -  1479  156s   364   265    0.01314   10   11          -    0.01314      -  1468  161s   390   276 infeasible    9               -    0.01314      -  1391  174s   417   284    0.01358   11   10          -    0.01314      -  1479  181s   443   293    0.01395   12   12          -    0.01314      -  1529  220s   472   331    0.01387   12   10          -    0.01314      -  1467  226s   526   345    0.01314   12   11          -    0.01314      -  1524  282s   562   356    0.01340   13   14          -    0.01314      -  1435  286s   601   408 infeasible   14               -    0.01314      -  1353  422s   673   448 infeasible   15               -    0.01314      -  1666  457s   745   473 infeasible   15               -    0.01314      -  1505  461s   814   485    0.01368   18   12          -    0.01314      -  1429  490s   862   493    0.01323   19   11          -    0.01314      -  1363  502s   878   503    0.01364   21   11          -    0.01314      -  1378  534s   896   511    0.01332   23   10          -    0.01314      -  1554  578s   908   516    0.01435   23   14          -    0.01314      -  1569  587s   923   521    0.01330   26   14          -    0.01314      -  1590  607s   942   606    0.01748   33   12          -    0.01314      -  1605  629s  1058   608    0.01748   95    0          -    0.01362      -  1535  676s  1060   610    0.01402    8    0          -    0.01362      -  1532  695s  1075   619    0.01362   17   19          -    0.01362      -  1570  720s  1083   625    0.01362   18   19          -    0.01362      -  1580  766s  1088   629    0.01457   20   16          -    0.01362      -  1637  795s  1094   633    0.01495   20   37          -    0.01362      -  1666  800s  1100   637    0.01362   19   19          -    0.01362      -  1662  855s  1110   640    0.01362   20   19          -    0.01362      -  1772  891s  1116   644    0.01362   21   19          -    0.01362      -  1854  922s  1126   642    0.01362   23   14          -    0.01362      -  1919  934s  1131   658    0.01362   24   12          -    0.01362      -  1921  959s  1159   687 infeasible   25               -    0.01362      -  1922 1019s  1211   690    0.01362   25  376          -    0.01362      -  1973 1080s  1264   704 infeasible   26               -    0.01362      -  2018 1143s  1339   780    0.01380   25   19          -    0.01362      -  2013 1161s  1472   860    0.01389   20   13          -    0.01362      -  1857 1193s  1601   963    0.01409   25   11          -    0.01362      -  1763 1205s  1821  1019    0.01462   44   10          -    0.01362      -  1561 1266s  2050  1028    0.01399   25   11          -    0.01362      -  1429 1295s* 2057   597             171       0.0259407    0.01362  47.5%  1434 1295s* 2070   564             177       0.0258932    0.01362  47.4%  1425 1295s  2211   442    0.01362   27   14    0.02589    0.01362  47.4%  1367 1335s  2233   447    0.01440   28   11    0.02589    0.01362  47.4%  1396 1368s  2271   461    0.01489   29   14    0.02589    0.01362  47.4%  1407 1395s  2380   455 infeasible   37         0.02589    0.01362  47.4%  1388 1475s  2474   546    0.01411   44   10    0.02589    0.01362  47.4%  1388 1498s  2665   525    0.01477   81   11    0.02589    0.01362  47.4%  1318 1538s  2825   529    0.01477  133    1    0.02589    0.01362  47.4%  1277 1598s* 2826   306             134       0.0147692    0.01362  7.81%  1277 1598s  2938   330     cutoff   64         0.01477    0.01362  7.81%  1258 1685s  3036   342     cutoff   63         0.01477    0.01362  7.81%  1273 1720s  3083   426    0.01365   34   10    0.01477    0.01362  7.81%  1300 1812s  3231   434 infeasible   44         0.01477    0.01362  7.81%  1284 1858s  3277   436     cutoff   34         0.01477    0.01362  7.81%  1312 1910s  3349   449    0.01362   24   18    0.01477    0.01362  7.81%  1331 1972s  3410   446 infeasible   29         0.01477    0.01362  7.81%  1355 2031s  3447   449     cutoff   29         0.01477    0.01362  7.81%  1374 2095s  3506   455    0.01362   22   45    0.01477    0.01362  7.81%  1413 2182s  3576   457    0.01362   35   46    0.01477    0.01362  7.81%  1443 2257s  3606   467     cutoff   35         0.01477    0.01362  7.81%  1469 2314s  3654   474    0.01419   46   18    0.01477    0.01362  7.81%  1481 2400s  3691   466 infeasible   32         0.01477    0.01362  7.81%  1525 2444s  3739   486    0.01366   31   14    0.01477    0.01362  7.80%  1595 2594s  3841   504    0.01362   32  489    0.01477    0.01362  7.80%  1593 2664s  3895   516    0.01400   33   19    0.01477    0.01362  7.80%  1609 2707s  3989   522    0.01416   34   20    0.01477    0.01362  7.80%  1601 2840s  4051   519     cutoff   37         0.01477    0.01366  7.54%  1633 2915s  4084   526 infeasible   37         0.01477    0.01372  7.09%  1660 2964s  4183   545    0.01396   29   14    0.01477    0.01372  7.09%  1666 3009s  4300   571    0.01372   27   17    0.01477    0.01372  7.09%  1656 3046s  4422   591    0.01406   23   12    0.01477    0.01372  7.09%  1654 3084s  4518   596     cutoff   27         0.01477    0.01372  7.09%  1664 3153s  4593   602    0.01419   51   18    0.01477    0.01372  7.09%  1675 3212s  4655   626    0.01374   27   14    0.01477    0.01374  6.95%  1695 3289s  4745   653    0.01374   30   15    0.01477    0.01374  6.95%  1709 3364s  4918   660 infeasible   39         0.01477    0.01377  6.80%  1704 3448s  5029   663    0.01377   42   23    0.01477    0.01377  6.80%  1734 3535s  5130   667     cutoff   35         0.01477    0.01377  6.79%  1767 3605s  5218   676    0.01379   24   13    0.01477    0.01379  6.61%  1796 3675s  5365   672     cutoff   34         0.01477    0.01381  6.47%  1793 3766s  5459   669    0.01477   22   13    0.01477    0.01384  6.32%  1836 3975s  5540   668     cutoff   23         0.01477    0.01384  6.32%  1871 4051s  5589   658     cutoff   31         0.01477    0.01388  5.99%  1917 4116s  5625   658    0.01434   27   12    0.01477    0.01389  5.98%  1952 4199s  5731   658 infeasible   41         0.01477    0.01395  5.57%  1967 4277s  5795   661    0.01419   55   18    0.01477    0.01396  5.46%  1995 4338s  5914   667     cutoff   30         0.01477    0.01396  5.46%  2012 4420s  6050   665 infeasible   39         0.01477    0.01397  5.44%  2029 4482s  6174   664    0.01409   27   12    0.01477    0.01397  5.43%  2051 4553s  6281   670    0.01397   29   11    0.01477    0.01397  5.42%  2082 4775s  6403   703    0.01410   43   10    0.01477    0.01397  5.42%  2097 4879s  6586   725    0.01419   55   18    0.01477    0.01399  5.27%  2103 4956s  6824   752 infeasible   32         0.01477    0.01399  5.27%  2103 5048s  7143   844    0.01399   34   12    0.01477    0.01399  5.27%  2075 5235s  7847  1134    0.01401   34   13    0.01477    0.01401  5.11%  1958 5303s* 8144   506             113       0.0084364    0.00844  0.00%  1911 5303sExplored 9203 nodes (15696379 simplex iterations) in 5303.68 secondsThread count was 8 (of 8 available processors)Solution count 4: 0.00843635 0.0147692 0.0258932 0.0259407Optimal solution found (tolerance 1.00e-04)Warning: max constraint violation (5.4127e-09) exceeds toleranceBest objective 8.436353171704e-03, best bound 8.436353171704e-03, gap 0.0000%Obj:  0.008436353171704492q:  {0: <gurobi.Var q[0] (value 267.562984392664)>, 1: <gurobi.Var q[1] (value 55.53033496620341)>}Solution quality statistics for model MyModel2 :  Maximum violation:    Bound       : 0.00000000e+00    Constraint  : 5.41270806e-09 (QC250)    Integrality : 0.00000000e+00
I want to know if my problem reformulation is correct, how to solve this problem, and how to determine if the output solution is correct.

All the best,

Yali

• Gurobi Staff

Your formulation of using auxiliary variables to translate your model into a form which can be implemented using gurobipy looks fine to me.

how to determine if the output solution is correct.

• To some extent, you need to trust that Gurobi provides the optimal solution. You can of course do a number of different sanity checks:
• For a smaller size instance of your model, you can find the optimal solution via enumerating all possible solutions and check that the optimal solution found matches the solution returned by Gurobi.
• If Gurobi reports $$z^*$$ as the optimal objective value for a minimization problem, you can try to solve the model again with the additional constraint $$z \leq z^* - \epsilon$$. We expect for Gurobi to terminate with an infeasibility status because Gurobi has already concluded that it is not possible to have a solution better than $$z^*$$.

Best regards,

Maliheh

Dear Maliheh,

The optimization problem is as follows, and optimization variables are underlined in blue.

Currently, our code and existing issues are as follows.

from gurobipy import *import gurobipy as gpfrom numpy import *import numpy as npimport mathimport random#define global variablesglobal KK = 10global MM = 5global ATtheta1ATtheta1=0.42global ATtheta2ATtheta2=2.00global CPtheta1CPtheta1=0.11global CPtheta2CPtheta2=0.68global CStheta1CStheta1=0.05global CStheta2CStheta2=0.62global HH=100global PmaxPmax=10**(-2)  global N0w#N0w=10**(-11)N0w=0global beta0beta0=10**(-4)  global thresholdthreshold=10**(-3)global sumcombinationsumcombination=int(math.factorial(K)/math.factorial(M)/math.factorial(K-M))#The Positions of K users are given randomlyglobal uposupos = []for i in range(0,K):    upos.append([400*random.random(),400*random.random()])#define distancedef 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 distanceglobal DSsetDSset= []global DSmaxfor i in range(0,K):    for j in range(0,K):            DSset.append(DS(i,j))DSmax=max(DSset)  #define functiondef cov(i,j):    return math.exp(-((DS(i,j)/DSmax/ATtheta1)**ATtheta2))global deltazsum0=0for k in range(0,K):    for k1 in range(0,K):            sum0=sum0+cov(k,k1)deltaz=sum0/K/K#initialized variablesb0=math.sqrt(Pmax)*np.array(ones((K)))q0=np.array([250,250])#create modelmodel = gp.Model("MyModel2")#declare variablesq=model.addVars(2, lb=0, ub=400, vtype=GRB.CONTINUOUS, name="q")a=model.addVars(K, lb=0, ub=(400*math.sqrt(2))**2, name="a")v=model.addVars(K, lb=beta0/(H**2+(400*math.sqrt(2))**2), ub=beta0/H/H, name="v")r=model.addVars(K, lb=math.sqrt(beta0/(H**2+(400*math.sqrt(2))**2)), ub=math.sqrt(beta0/H/H), name="r")t=model.addVars(K, K, lb=beta0/(H**2+(400*math.sqrt(2))**2), ub=beta0/H/H, name="t")nu=model.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=0for i in range(0, sumcombination):      sumvalue=sumvalue+deltaz-nu[i]obj=math.factorial(M)*math.factorial(K-M)/math.factorial(K)*sumvaluemodel.setObjective(obj, GRB.MINIMIZE)model.addConstrs(a[k]==(q[0]-upos[k][0])**2+(q[1]-upos[k][1])**2 for k in range(0,K))model.addConstrs(v[k]*(H**2+a[k])/beta0 == 1 for k in range(0,K))#for k in range(0,K):#    model.addGenConstrPow(v[k], r[k], 0.5, "FuncPieces=-2 FuncPieceError=0.0001")model.addConstrs(v[k] == (r[k])**2 for k in range(0,K))for k in range(0,K):    for k1 in range(0,K):        model.addConstr(t[k,k1] == r[k]*r[k1])combination=[]index=0def 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]))          model.addConstr(sum1**2 == nu[index]*K*K*(sum2+N0w))        index=index+1        return    for idx in range(last_idx+1,K):        combination.append(idx)        dfs1(idx, cnt+1)        combination.pop()dfs1(-1,0)model.params.NonConvex=2model.params.NumericFocus=3model.params.FeasibilityTol=1e-9model.params.OptimalityTol=1e-9model.params.ScaleFlag=2model.params.Aggregate=0model.optimize()print("Obj: ", model.objVal)print("q: ", q)opt=model.objValmodel.printQuality()

For the constraint (c), I would like to ask a question. Is it appropriate to write as

for k in range(0,K):    model.addGenConstrPow(v[k], r[k], 0.5, "FuncPieces=-2 FuncPieceError=0.0001")
or as
model.addConstrs(v[k] == (r[k])**2 for k in range(0,K))
(It should be noted that r[k] is already a positive value in the variable definition)
The problems with the former writing method are as follows.
Changed value of parameter NonConvex to 2   Prev: -1  Min: -1  Max: 2  Default: -1Changed value of parameter NumericFocus to 3   Prev: 0  Min: 0  Max: 3  Default: 0Changed value of parameter FeasibilityTol to 1e-09   Prev: 1e-06  Min: 1e-09  Max: 0.01  Default: 1e-06Changed value of parameter OptimalityTol to 1e-09   Prev: 1e-06  Min: 1e-09  Max: 0.01  Default: 1e-06Changed value of parameter ScaleFlag to 2   Prev: -1  Min: -1  Max: 3  Default: -1Changed value of parameter Aggregate to 0   Prev: 1  Min: 0  Max: 1  Default: 1Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)Thread count: 4 physical cores, 8 logical processors, using up to 8 threadsOptimize a model with 0 rows, 384 columns and 0 nonzerosModel fingerprint: 0x640229d0Model has 372 quadratic constraintsModel has 10 general constraintsVariable types: 384 continuous, 0 integer (0 binary)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     [3e-10, 3e+05]  RHS range        [0e+00, 0e+00]  QRHS range       [1e+00, 2e+05]Presolve removed 0 rows and 22 columnsPresolve time: 0.02sPresolved: 25842 rows, 6662 columns, 61620 nonzerosPresolved model has 6400 bilinear constraint(s)Variable types: 6662 continuous, 0 integer (0 binary)Root relaxation: objective 1.248260e-02, 13165 iterations, 1.81 seconds    Nodes    |    Current Node    |     Objective Bounds      |     Work Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time     0     0    0.01248    0   15          -    0.01248      -     -    1sH    0     0                       0.0251561    0.01248  50.4%     -    2s     0     2    0.01248    0   15    0.02516    0.01248  50.4%     -    2s*   14    14               3       0.0130558    0.01304  0.14%   422    4sExplored 23 nodes (19937 simplex iterations) in 4.89 secondsThread count was 8 (of 8 available processors)Solution count 2: 0.0130558 0.0251561Optimal solution found (tolerance 1.00e-04)Warning: max constraint violation (5.1011e-09) exceeds toleranceWarning: max general constraint violation (3.7379e-09) exceeds toleranceBest objective 1.305580428846e-02, best bound 1.305580428846e-02, gap 0.0000%Obj:  0.013055804288459483q:  {0: <gurobi.Var q[0] (value 0.0)>, 1: <gurobi.Var q[1] (value 0.0)>}Solution quality statistics for model MyModel2 :  Maximum violation:    Bound       : 0.00000000e+00    Constraint  : 5.10105021e-09 (QC226)    Integrality : 0.00000000e+00
I would like to ask you how to solve the two warning issues.

Then, the problems with the latter writing method are as follows.
Changed value of parameter NonConvex to 2   Prev: -1  Min: -1  Max: 2  Default: -1Changed value of parameter NumericFocus to 3   Prev: 0  Min: 0  Max: 3  Default: 0Changed value of parameter FeasibilityTol to 1e-09   Prev: 1e-06  Min: 1e-09  Max: 0.01  Default: 1e-06Changed value of parameter OptimalityTol to 1e-09   Prev: 1e-06  Min: 1e-09  Max: 0.01  Default: 1e-06Changed value of parameter ScaleFlag to 2   Prev: -1  Min: -1  Max: 3  Default: -1Changed value of parameter Aggregate to 0   Prev: 1  Min: 0  Max: 1  Default: 1Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)Thread count: 4 physical cores, 8 logical processors, using up to 8 threadsOptimize a model with 0 rows, 384 columns and 0 nonzerosModel fingerprint: 0x04e81293Model has 382 quadratic constraintsCoefficient 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     [3e-10, 3e+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.02sPresolved: 25928 rows, 6686 columns, 61820 nonzerosPresolved model has 6422 bilinear constraint(s)Variable types: 6686 continuous, 0 integer (0 binary)Root simplex log...Iteration    Objective       Primal Inf.    Dual Inf.      Time   17237   -1.8045301e-02   2.306531e-04   3.255809e+11      5s   21918    1.8956010e-02   3.334482e-04   1.765671e+12     10s   29181    8.9006723e-04   9.741333e-05   0.000000e+00     15s   34863    1.1348627e-02   0.000000e+00   0.000000e+00     18sRoot relaxation: objective 1.134863e-02, 34863 iterations, 17.89 seconds    Nodes    |    Current Node    |     Objective Bounds      |     Work Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time     0     0    0.01135    0   10          -    0.01135      -     -   17s     0     2    0.01135    0   10          -    0.01135      -     -   19s     3     6 infeasible    2               -    0.01135      -  86.7   20s    15    14 infeasible    4               -    0.01135      -   132   37s    19    16    0.01165    5   11          -    0.01135      -  2239   44s    23    18 infeasible    5               -    0.01135      -  2305   72s    29    21 infeasible    6               -    0.01135      -  3371   92s    34    25    0.01135    6   20          -    0.01135      -  5219  125s    40    33    0.01135    7   20          -    0.01135      -  5483  207s    48    35    0.01135    8   27          -    0.01135      -  5608  217s    66    54    0.01135   10   20          -    0.01135      -  4301  228s    79    58    0.01135   11   20          -    0.01135      -  3727  234s    93    62    0.01135   12   20          -    0.01135      -  3249  246s   107    83    0.01135   13   10          -    0.01135      -  3080  268s   130    99    0.01155   14   10          -    0.01135      -  2851  272s   152   118    0.01135   15   10          -    0.01135      -  2443  276s   177   142    0.01189   18   16          -    0.01135      -  2174  281s   209   166    0.01135   20   21          -    0.01135      -  1923  301s   241   180    0.01135   22   21          -    0.01135      -  1825  323s   269   207    0.01337   26   12          -    0.01135      -  1828  336s   300   244    0.01424   33   10          -    0.01135      -  1709  358s   337   296    0.01424   41   10          -    0.01135      -  1578  360s   405   363    0.01424   46   10          -    0.01135      -  1319  371s   488   429    0.01424   56   10          -    0.01135      -  1131  380s   574   493    0.01424   69   10          -    0.01135      -   991  394s*  665   302              60       0.0121132    0.01135  6.31%   888  394s   666   234     cutoff   84         0.01211    0.01135  6.31%   887  409s   705   258 infeasible   11         0.01211    0.01135  6.31%   905  428s   786   270    0.01137   36   10    0.01211    0.01135  6.31%   849  440s   965   334 infeasible   36         0.01211    0.01135  6.31%   727  456s* 1048   330              56       0.0121093    0.01135  6.28%   673  456s  1076   363    0.01198   63    6    0.01211    0.01135  6.28%   669  464s  1213   446    0.01198   87    6    0.01211    0.01135  6.28%   615  468s  1440   584    0.01204   15   14    0.01211    0.01135  6.28%   523  471s* 1579   352              60       0.0115466    0.01135  1.71%   478  471s  1660   199     cutoff   16         0.01155    0.01135  1.71%   457  489s  1743   208    0.01135   15   23    0.01155    0.01135  1.71%   460  504s  1782   205 infeasible   19         0.01155    0.01135  1.71%   474  511s  1831   206     cutoff   22         0.01155    0.01135  1.71%   477  518s  1874   220 infeasible   15         0.01155    0.01135  1.71%   480  524s  1910   290     cutoff   25         0.01155    0.01135  1.71%   479  570s  2015   224     cutoff  106         0.01155    0.01135  1.71%   490  600s  2059   225    0.01149   17    0    0.01155    0.01135  1.71%   528  613s  2061   226    0.01149   43    0    0.01155    0.01149  0.47%   528  650sExplored 2065 nodes (1133934 simplex iterations) in 652.39 secondsThread count was 8 (of 8 available processors)Solution count 3: 0.0115466 0.0121093 0.0121132Optimal solution found (tolerance 1.00e-04)Warning: max constraint violation (4.7465e-09) exceeds toleranceBest objective 1.154663448745e-02, best bound 1.154663448745e-02, gap 0.0000%Obj:  0.01154663448744575q:  {0: <gurobi.Var q[0] (value 0.01176034045085275)>, 1: <gurobi.Var q[1] (value 238.42679656337583)>}Solution quality statistics for model MyModel2 :  Maximum violation:    Bound       : 0.00000000e+00    Constraint  : 4.74646360e-09 (QC316)    Integrality : 0.00000000e+00

I would like to ask you how to solve this warning issue.

In summary, I would like to ask you three questions.

(1) How to write the constraint (c)? Both methods should be correct, but the warning issues are different.

(2) For the former writing method, how to deal with "Warning: max constraint violation (5.1011e-09) exceeds tolerance" and "Warning: max general constraint violation (3.7379e-09) exceeds tolerance"

(3) For the latter writing method, how to deal with "Warning: max constraint violation (4.7465e-09) exceeds tolerance"

Best regards,

Yali

Dear Maliheh,

About the questions listed above, we really need and look forward to your comments. Thank you very much！

Best regards,

Yali

• Gurobi Staff

(1) How to write the constraint (c)? Both methods should be correct, but the warning issues are different.

First approach: when a general function constraint such as power function is used to model a nonlinear relationship (square root), Gurobi will automatically add a piecewise-linear approximation of the function to the model.

Second approach: a quadratic equality constraint would be handled directly by Gurobi the as a non-convex quadratic constraint.

(2) For the former writing method, how to deal with "Warning: max constraint violation (5.1011e-09) exceeds tolerance" and "Warning: max general constraint violation (3.7379e-09) exceeds tolerance"

(3) For the latter writing method, how to deal with "Warning: max constraint violation (4.7465e-09) exceeds tolerance

Gurobi solves the model in the scaled presolved space. Both scaling and presolve have implicit effects on the feasibility tolerance value. This can lead to constraint and bound violations beyond the feasibility tolerance value 1e-9 when the solution is mapped back to the unscaled and original problem space. Are violations as small as 5e-9 unacceptable for your problem?

• The log shows that there are variables with bounds as small as 3e-10 in your problem. Does such a small value have any physical meaning in your problem? If not, I would suggest zeroing it out. If it does, you would need to consider scaling your variables.
• You can experiment with lower values for the ScaleFlag parameter but this can increase the runtime.
• For question 2, you can also experiment with parameters/attributes such as FuncPieces, FuncPieceLength, FuncPieceError, FuncPieceRatio which control the accuracy of the piecewise linear approximation. Of course, there is a speed-versus-accuracy tradeoff here. For example, the bigger the FuncPieces parameter is, the more accurate the approximation is at the cost of longer runtime. (See section Function Constraints in General Constraints documentation)

Best regards,

Maliheh

Hi, Maliheh,

(1) Regarding the "Warning: max constraint violation (4.7465e-09) exceeds tolerance", as the violation is relatively small, if I can accept it, can the output result be guaranteed to be correct? Warnings sometimes do not need to be addressed?

(2) Why did I explicitly add constraints to make the objective function positive, but the output result was negative and there was no error and no warnings reported.

Solution count 4: -7.77156e-16 0.000672095 0.00068015 0.000683127Optimal solution found (tolerance 1.00e-04)Warning: max constraint violation (6.4220e-07) exceeds toleranceBest objective -7.771561172376e-16, best bound -7.771561172376e-16, gap 0.0000%Obj:  -7.771561172376096e-16

Best regards,

Yali

Hi, Maliheh,

In addition to the two questions I asked you above, I also want to ask you new questions.

The optimization problem is shown as follows. Optimization variables are underlined in red.

To facilitate the solution, we introduce auxiliary variables. Then, the optimizaion problem is reformulatd as folllows. Optimization variables are underlined in blue.

My code is as follows，

from gurobipy import *import gurobipy as gpfrom numpy import *import numpy as npimport mathimport random#define global variablesglobal KK = 10global MM = 5global ATtheta1ATtheta1=0.42global ATtheta2ATtheta2=2.00global CPtheta1CPtheta1=0.11global CPtheta2CPtheta2=0.68global CStheta1CStheta1=0.05global CStheta2CStheta2=0.62global HH=100global PmaxPmax=10**(-2)  global N0w#N0w=10**(-11)N0w=0global beta0beta0=10**(-4)  global thresholdthreshold=10**(-3)global sumcombinationsumcombination=int(math.factorial(K)/math.factorial(M)/math.factorial(K-M))#The Positions of K users are given randomlyglobal uposupos = []for i in range(0,K):    upos.append([400*random.random(),400*random.random()])#define distancedef 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 distanceglobal DSsetDSset= []global DSmaxfor i in range(0,K):    for j in range(0,K):            DSset.append(DS(i,j))DSmax=max(DSset)  #define functiondef cov(i,j):    return math.exp(-((DS(i,j)/DSmax/ATtheta1)**ATtheta2))global deltazsum0=0for k in range(0,K):    for k1 in range(0,K):            sum0=sum0+cov(k,k1)deltaz=sum0/K/K    q=np.array([250,250])#create modelmodel2 = gp.Model("MyModel2")#declare variablesb=model2.addVars(K, lb=0, ub=math.sqrt(Pmax), name="b")t1=model2.addVars(K, K, lb=0, ub=Pmax, name="t1")nu1=model2.addVars(sumcombination, lb=0, ub=deltaz, name="nu1")#set objectiveobj = LinExpr(0)for i in range(sumcombination):    obj.addTerms(-math.factorial(M)*math.factorial(K-M)/math.factorial(K), nu1[i])model2.setObjective(obj+deltaz, GRB.MINIMIZE)for k in range(0,K):    for k1 in range(0,K):        model2.addConstr(t1[k,k1] == b[k]*b[k1])combination=[]index=0def 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+(b[combination[km]]*math.sqrt(beta0/(H**2+(q[0]-upos[combination[km]][0])**2+(q[1]-upos[combination[km]][1])**2))*cov(k,combination[km]))                          sum2=0        for km in range(0,M):            for km1 in range(0,M):                sum2=sum2+(t1[combination[km],combination[km1]]*math.sqrt(beta0/(H**2+(q[0]-upos[combination[km]][0])**2+(q[1]-upos[combination[km]][1])**2))*math.sqrt(beta0/(H**2+(q[0]-upos[combination[km1]][0])**2+(q[1]-upos[combination[km1]][1])**2))*cov(combination[km],combination[km1]))            model2.addConstr(sum1**2 == nu1[index]*K*K*(sum2+N0w))        index=index+1        return    for idx in range(last_idx+1,K):        combination.append(idx)        dfs1(idx, cnt+1)        combination.pop()dfs1(-1,0)model2.params.NonConvex=2model2.params.NumericFocus=2model2.params.FeasibilityTol=1e-9model2.params.OptimalityTol=1e-9model2.params.ScaleFlag=0model2.params.Aggregate=0model2.optimize()print("Obj: ", model2.objVal)print("b: ", b)model2.printQuality()

The output results are as follows，

Set parameter NonConvex to value 2Set parameter NumericFocus to value 2Set parameter FeasibilityTol to value 1e-09Set parameter OptimalityTol to value 1e-09Set parameter ScaleFlag to value 0Set parameter Aggregate to value 0Gurobi 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 threadsOptimize a model with 0 rows, 362 columns and 0 nonzerosModel fingerprint: 0x48473f75Model has 352 quadratic constraintsCoefficient statistics:  Matrix range     [0e+00, 0e+00]  QMatrix range    [7e-10, 1e+00]  QLMatrix range   [1e+00, 1e+00]  Objective range  [4e-03, 4e-03]  Bounds range     [1e-02, 4e-01]  RHS range        [0e+00, 0e+00]Continuous model is non-convex -- solving as a MIPFound heuristic solution: objective 0.3558289Presolve time: 0.00sPresolved: 25842 rows, 6662 columns, 61250 nonzerosPresolved model has 6400 bilinear constraint(s)Variable types: 6662 continuous, 0 integer (0 binary)Root relaxation: objective 5.551115e-17, 0 iterations, 0.01 seconds (0.00 work units)    Nodes    |    Current Node    |     Objective Bounds      |     Work Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time*    0     0               0      -0.0000000   -0.00000  0.00%     -    0sExplored 1 nodes (0 simplex iterations) in 0.13 seconds (0.06 work units)Thread count was 8 (of 8 available processors)Solution count 2: -1.11022e-16 0.355829Optimal solution found (tolerance 1.00e-04)Best objective -1.665334536938e-16, best bound -1.665334536938e-16, gap 0.0000%Obj:  -1.6653345369377348e-16b:  {0: <gurobi.Var b[0] (value 0.0)>, 1: <gurobi.Var b[1] (value 0.0)>, 2: <gurobi.Var b[2] (value 0.0)>, 3: <gurobi.Var b[3] (value 0.0)>, 4: <gurobi.Var b[4] (value 0.0)>, 5: <gurobi.Var b[5] (value 0.0)>, 6: <gurobi.Var b[6] (value 0.0)>, 7: <gurobi.Var b[7] (value 0.0)>, 8: <gurobi.Var b[8] (value 0.0)>, 9: <gurobi.Var b[9] (value 0.0)>}Solution quality statistics for model MyModel2 :  Maximum violation:    Bound       : 0.00000000e+00    Constraint  : 0.00000000e+00    Integrality : 0.00000000e+00

Firstly, the variable “b” should be greater than 0, not equal to 0. The result is incorrect.

Secondly, according to the current output, when b is equal to 0, the result should be 0.355829 in "Solution count 2: -1.11022e-16 0.355829", rather than the current optimal solution.

Thirdly, for the optimization variable, as it is greater than 0, when we set "lb" to 1e-9 rather than 0, the output result is still similar, and the problem still exists.

Best regards,

Yali

Hi, Maliheh,

In addition to the above questions, I have accumulated additional questions to ask you. How to scaling the variables? For example, when I multiply a variable "b" with a very small lower bound 3e-10 by 1e5, i.e., the lower bound of "b" is 3e-5. Then, I need to divide the "b" by 1e5 in the code to ensure that the output result does not change? Looking forward to your comments to the questions in the three posts mentioned above.

Best regards,

Yali

• Gurobi Staff

(1) Regarding the "Warning: max constraint violation (4.7465e-09) exceeds tolerance", as the violation is relatively small, if I can accept it, can the output result be guaranteed to be correct? Warnings sometimes do not need to be addressed?

Gurobi is a floating-point based MIP solver because computers represent real numbers within finite precision. Imprecise arithmetic results in round-off errors. To deal with this reality, the tolerance values need to be introduced to allow for potential acceptance of a solution with a violation of at most $$\epsilon$$ being equal to the feasibility tolerance. Whether a violation as small as 4.7465e-09 is acceptable or not is a question that needs to be answered by the user of the model.

(2) Why did I explicitly add constraints to make the objective function positive, but the output result was negative and there was no error and no warnings reported.

Firstly, the variable “b” should be greater than 0, not equal to 0. The result is incorrect.

• It is not possible to have strict inequality constraints in Gurobi (see "Why doesn't Gurobi support strict inequality constraints like x < a?" for more information).
• The variable $$b$$ is defined with lower bound 0 in your model. So, it can take value 0.
• I do not see any constraint that enforces the objective to be strictly positive in your model. Furthermore, any constraint in the form $$Ax \geq b$$ is considered satisfied if $$Ax - b \geq -\mbox{FeasibilityTol}$$.

Secondly, according to the current output, when b is equal to 0, the result should be 0.355829 in "Solution count 2: -1.11022e-16 0.355829", rather than the current optimal solution.

• If all $$b$$ decision variables are zero, according to constraints $$(b)$$, all decision variables $$t$$ are also zero.
• According to constraints $$(c)$$, you then have $$0 = 0 + K^2 \sigma_w^2 \nu_{\mathcal{I}}$$. Since the $$\sigma_w^2$$ parameter is set to 0, you then have $$0 = 0$$. This means that the variable $$\nu_{\mathcal{I}}$$ can take any value within its bounds without violating any constraints. To minimize the objective function, the optimal value for the $$\nu_{\mathcal{I}}$$ variable equals to its upper bound $$\mbox{deltaz}$$, resulting in an optimal objective value of 0. If you would like to force $$\nu_{\mathcal{I}}$$ to equal 0 if $$b$$ and $$t$$ variables are 0, you need to set the $$\sigma_w^2$$ parameter to a value strictly greater than 0 (few orders of magnitude larger than the feasibility tolerance).

Thirdly, for the optimization variable, as it is greater than 0, when we set "lb" to 1e-9 rather than 0, the output result is still similar, and the problem still exists.

Defining a lower bound for a variable essentially means having $$x \geq l$$. The lower bound of a variable is considered to be respected as long as $$x - l \geq - \mbox{FeasibilityTol}$$.

For example, when I multiply a variable "b" with a very small lower bound 3e-10 by 1e5, i.e., the lower bound of "b" is 3e-5. Then, I need to divide the "b" by 1e5 in the code to ensure that the output result does not change?

Let us assume that you would like to scale the variable $$x$$ as $$x^{\prime} = 10^5 x$$. You need to replace the variable $$x$$ by $$10^{-5} x^{\prime}$$. The lower and upper bounds of the variable $$x^{\prime}$$ should be defined by multiplying the lower and upper bounds of the variable $$x$$ by $$10^5$$.

Best regards,

Maliheh

Hi, Maliheh,

(1) It is not possible to have strict inequality constraints in Gurobi, if the variable "b>0" not "b>=0", how should we express it? Is it represented in the lower bound or in the constraint?

(2) "If you would like to force νI to equal 0 if b and t variables are 0, you need to set the σw2 parameter to a value strictly greater than 0 (few orders of magnitude larger than the feasibility tolerance)." the σw2 parameter is noise, It should have been taken as 1e-11, but due to numerical issues caused by being too small, it was changed to 0. Can it be arbitrarily changed to a number which is few orders of magnitude larger than the feasibility tolerance, such as 1e-8?

(3) How to scaling the variables is effctive?

(4) How to solve the waring "Warning: 1 variables dropped from basis

Best regards,

Yali

• Gurobi Staff

(1) It is not possible to have strict inequality constraints in Gurobi, if the variable "b>0" not "b>=0", how should we express it? Is it represented in the lower bound or in the constraint?

To model $$x > 0$$, you need to model $$x \geq \epsilon$$ either as a constraint or set the lower bound of the variable $$x$$ to $$\epsilon$$. The value of $$\epsilon$$ depends on your application and it should be at least one order of magnitude larger than the value of the FeasibilityTol parameter.

(2) "If you would like to force νI to equal 0 if b and t variables are 0, you need to set the σw2 parameter to a value strictly greater than 0 (few orders of magnitude larger than the feasibility tolerance)." the σw2 parameter is noise, It should have been taken as 1e-11, but due to numerical issues caused by being too small, it was changed to 0. Can it be arbitrarily changed to a number which is few orders of magnitude larger than the feasibility tolerance, such as 1e-8?

I was just explaining why the optimal objective function was zero and not the value you expected to see.

• Ideally, the model should not include any parameter/coefficient smaller than the feasibility tolerance value. The rule of thumb is for the smallest coefficient to be at least one order of magnitude bigger.

(3) How to scaling the variables is effctive?

One source of numerical issues in the model is bad numbers such as large matrix coefficient range, large objective coefficient values, and large RHS values. Scaling the variables/constraints helps to get reasonable coefficients and address this source of numerical issues in the model.

(4) How to solve the waring "Warning: 1 variables dropped from basis

I did not see such warnings in your log. These warning messages indicate that the basis matrix got near singular. The Gurobi Optimizer then replaced one of the structural variables by a slack or switched to a higher precision to address the singularity of the basis matrix.

Hope this concludes all your questions.

Best regards,

Maliheh

Hi, Maliheh,
My optimization problem is as follows, and optimization variables are underlined in red.

First, because the lower bound of my variables v and t can reach 1e-10, which is smaller than the "model1.params.FeasibilityTol=1e-9", 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")
Second, 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.

My complete code is as follows,

from gurobipy import *import gurobipy as gpfrom numpy import *import numpy as npimport mathimport random#define global variablesglobal KK = 10global MM = 5global ATtheta1ATtheta1=0.42global ATtheta2ATtheta2=2.00global HH=100global PmaxPmax=10**(-2)  global N0w#N0w=10**(-11)  #to avoid large coefficients numerical issueN0w=0global beta0beta0=10**(-4)  global sumcombinationsumcombination=int(math.factorial(K)/math.factorial(M)/math.factorial(K-M))#The Positions of K users are given randomlyglobal uposupos=[[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 distancedef 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 distanceglobal DSsetDSset= []global DSmaxfor i in range(0,K):    for j in range(0,K):            DSset.append(DS(i,j))DSmax=max(DSset)  #define functiondef cov(i,j):    return math.exp(-((DS(i,j)/DSmax/ATtheta1)**ATtheta2))global deltazsum0=0for k in range(0,K):    for k1 in range(0,K):            sum0=sum0+cov(k,k1)deltaz=sum0/K/K#initialized variablesb0=math.sqrt(Pmax)*np.array(ones((K)))#Obtain the maximum distance to tighten the range of variable aglobal DSset_aDSset_a= []global DSmax_afor 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 modelmodel1 = gp.Model("MyModel1")#declare variablesq=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=(10**3)*beta0/(H**2+DSmax_a), ub=(10**3)*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=(10**3)*beta0/(H**2+DSmax_a), ub=(10**3)*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=0for i in range(0, sumcombination):      sumvalue=sumvalue+deltaz-nu[i]obj=math.factorial(M)*math.factorial(K-M)/math.factorial(K)*sumvaluemodel1.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*(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")combination=[]index=0def 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]]*(10**(-3))*t[combination[km],combination[km1]]*cov(combination[km],combination[km1]))          model1.addConstr((sum1**2)*(10**3) == nu[index]*K*K*(sum2+N0w)*(10**3), name="c4")        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=2model1.params.NumericFocus=2model1.params.FeasibilityTol=1e-9model1.params.OptimalityTol=1e-9model1.params.ScaleFlag=2model1.params.Aggregate=0model1.optimize()print("Obj: ", model1.objVal)print("q: ", q)model1.printQuality()
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 threadsOptimize a model with 0 rows, 384 columns and 0 nonzerosModel fingerprint: 0x1643e273Model has 382 quadratic constraintsCoefficient 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-07, 2e+05]  RHS range        [0e+00, 0e+00]  QRHS range       [1e+03, 2e+05]Continuous model is non-convex -- solving as a MIPPresolve time: 0.01sPresolved: 26138 rows, 6741 columns, 87620 nonzerosPresolved model has 6477 bilinear constraint(s)Variable types: 6741 continuous, 0 integer (0 binary)Root relaxation: objective 9.714181e-04, 24208 iterations, 2.44 seconds (5.82 work units)
However, the code runs continuously, and the simulation time is very long and it is difficult to output results. I want to know what the reason is and how to solve it.

Best regards,

Yali

Hi, Maliheh,

Connect to the previous post,  Although scaling variables and "multiply both sides of the equation by 1e3" can avoid some warnings, they greatly prolong the simulation time, resulting in no results at all.

Best regards,

Yali

• Gurobi Staff

The difficulty of the new-version of your model is not related to scaling. It is because some of the input data to the model has changed, for example, the values of $$\texttt{upos}$$ parameters have changed. Even if you remove the scaling, it still takes time to find a feasible solution.

Please note that tightening the FeasibilityTol/OptimalityTol might come at the cost of a large number of iterations. I would suggest experimenting with longer TimeLimit and parameters such as MIPFocus=1, Heuristics, and NoRelHeurTime. You might also want to experiment with the Model.feasRelaxS() to figure out relaxing which variable bounds results in a feasible solution.

If you need any further help with the performance improvement of your model, please open a new post.

Best regards,

Maliheh