Skip to main content

Quadratic constraints contain large coefficient range on linear part/max constraint violation (1.0000e+00) exceeds tolerance

Answered

Comments

15 comments

  • Maliheh Aramon
    Gurobi Staff 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

    0
  • Yali Chen
    Curious
    Collaborator

    Hi, Maliheh

    Thank you very much for your detailed reply! 

    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 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 CPtheta1

    CPtheta1=0.11

    global CPtheta2

    CPtheta2=0.68

    global CStheta1

    CStheta1=0.05

    global CStheta2

    CStheta2=0.62

    global H

    H=100

    global Pmax

    Pmax=10**(-2)  

    global N0w

    #N0w=10**(-11)

    N0w=0

    global beta0

    beta0=10**(-4)  

    global threshold

    threshold=10**(-3)

    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 = []

    for i in range(0,K):

        upos.append([400*random.random(),400*random.random()])

    #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)))

    q0=np.array([250,250])

    #create model

    model = gp.Model("MyModel2")

    #declare variables

    q=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=0

    for i in range(0, sumcombination):  

        sumvalue=sumvalue+deltaz-nu[i]

    obj=math.factorial(M)*math.factorial(K-M)/math.factorial(K)*sumvalue

    model.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=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]))  

            #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=0

    model.params.NumericFocus=2

    model.params.FeasibilityTol=10**(-9)

    model.params.OptimalityTol=10**(-9)

    model.params.ScaleFlag=2

    model.params.Aggregate=0

    #

    model.optimize()

    print("Obj: ", model.objVal)

    print("q: ", q)

    opt=model.objVal

    model.printQuality()

    Changed value of parameter NonConvex to 2
       Prev: -1  Min: -1  Max: 2  Default: -1
    Changed value of parameter NumericFocus to 2
       Prev: 0  Min: 0  Max: 3  Default: 0
    Changed value of parameter FeasibilityTol to 1e-09
       Prev: 1e-06  Min: 1e-09  Max: 0.01  Default: 1e-06
    Changed value of parameter OptimalityTol to 1e-09
       Prev: 1e-06  Min: 1e-09  Max: 0.01  Default: 1e-06
    Changed value of parameter ScaleFlag to 2
       Prev: -1  Min: -1  Max: 3  Default: -1
    Changed value of parameter Aggregate to 0
       Prev: 1  Min: 0  Max: 1  Default: 1
    Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)
    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: 0x7abb6cb6
    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     [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.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.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 5303s

    Explored 9203 nodes (15696379 simplex iterations) in 5303.68 seconds
    Thread count was 8 (of 8 available processors)

    Solution count 4: 0.00843635 0.0147692 0.0258932 0.0259407

    Optimal solution found (tolerance 1.00e-04)
    Warning: max constraint violation (5.4127e-09) exceeds tolerance
    Best objective 8.436353171704e-03, best bound 8.436353171704e-03, gap 0.0000%
    Obj:  0.008436353171704492
    q:  {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.

    Looking forward to your reply! Many thanks!

    All the best,

    Yali

    0
  • Maliheh Aramon
    Gurobi Staff 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

    0
  • Yali Chen
    Curious
    Collaborator

    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 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 CPtheta1

    CPtheta1=0.11

    global CPtheta2

    CPtheta2=0.68

    global CStheta1

    CStheta1=0.05

    global CStheta2

    CStheta2=0.62

    global H

    H=100

    global Pmax

    Pmax=10**(-2)  

    global N0w

    #N0w=10**(-11)

    N0w=0

    global beta0

    beta0=10**(-4)  

    global threshold

    threshold=10**(-3)

    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 = []

    for i in range(0,K):

        upos.append([400*random.random(),400*random.random()])

    #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)))

    q0=np.array([250,250])

    #create model

    model = gp.Model("MyModel2")

    #declare variables

    q=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=0

    for i in range(0, sumcombination):  

        sumvalue=sumvalue+deltaz-nu[i]

    obj=math.factorial(M)*math.factorial(K-M)/math.factorial(K)*sumvalue

    model.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=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]))  

            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.NumericFocus=3

    model.params.FeasibilityTol=1e-9

    model.params.OptimalityTol=1e-9

    model.params.ScaleFlag=2

    model.params.Aggregate=0

    model.optimize()

    print("Obj: ", model.objVal)

    print("q: ", q)

    opt=model.objVal

    model.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: -1
    Changed value of parameter NumericFocus to 3
       Prev: 0  Min: 0  Max: 3  Default: 0
    Changed value of parameter FeasibilityTol to 1e-09
       Prev: 1e-06  Min: 1e-09  Max: 0.01  Default: 1e-06
    Changed value of parameter OptimalityTol to 1e-09
       Prev: 1e-06  Min: 1e-09  Max: 0.01  Default: 1e-06
    Changed value of parameter ScaleFlag to 2
       Prev: -1  Min: -1  Max: 3  Default: -1
    Changed value of parameter Aggregate to 0
       Prev: 1  Min: 0  Max: 1  Default: 1
    Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)
    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: 0x640229d0
    Model has 372 quadratic constraints
    Model has 10 general constraints
    Variable 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 columns
    Presolve time: 0.02s
    Presolved: 25842 rows, 6662 columns, 61620 nonzeros
    Presolved 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      -     -    1s
    H    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    4s

    Explored 23 nodes (19937 simplex iterations) in 4.89 seconds
    Thread count was 8 (of 8 available processors)

    Solution count 2: 0.0130558 0.0251561

    Optimal solution found (tolerance 1.00e-04)
    Warning: max constraint violation (5.1011e-09) exceeds tolerance
    Warning: max general constraint violation (3.7379e-09) exceeds tolerance
    Best objective 1.305580428846e-02, best bound 1.305580428846e-02, gap 0.0000%
    Obj:  0.013055804288459483
    q:  {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: -1
    Changed value of parameter NumericFocus to 3
       Prev: 0  Min: 0  Max: 3  Default: 0
    Changed value of parameter FeasibilityTol to 1e-09
       Prev: 1e-06  Min: 1e-09  Max: 0.01  Default: 1e-06
    Changed value of parameter OptimalityTol to 1e-09
       Prev: 1e-06  Min: 1e-09  Max: 0.01  Default: 1e-06
    Changed value of parameter ScaleFlag to 2
       Prev: -1  Min: -1  Max: 3  Default: -1
    Changed value of parameter Aggregate to 0
       Prev: 1  Min: 0  Max: 1  Default: 1
    Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)
    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: 0x04e81293
    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     [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.02s
    Presolved: 25928 rows, 6686 columns, 61820 nonzeros
    Presolved 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     18s

    Root 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  650s

    Explored 2065 nodes (1133934 simplex iterations) in 652.39 seconds
    Thread count was 8 (of 8 available processors)

    Solution count 3: 0.0115466 0.0121093 0.0121132

    Optimal solution found (tolerance 1.00e-04)
    Warning: max constraint violation (4.7465e-09) exceeds tolerance
    Best objective 1.154663448745e-02, best bound 1.154663448745e-02, gap 0.0000%
    Obj:  0.01154663448744575
    q:  {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"

    Looking forward to your reply! Many thanks!

    Best regards,

    Yali

     
     
    0
  • Yali Chen
    Curious
    Collaborator

    Dear Maliheh,

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

     

    Best regards,

    Yali

    0
  • Maliheh Aramon
    Gurobi Staff 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

    0
  • Yali Chen
    Curious
    Collaborator

    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.000683127

    Optimal solution found (tolerance 1.00e-04)
    Warning: max constraint violation (6.4220e-07) exceeds tolerance
    Best objective -7.771561172376e-16, best bound -7.771561172376e-16, gap 0.0000%
    Obj:  -7.771561172376096e-16

    Looking forward to your reply! Many thanks!

    Best regards,

    Yali

     

    0
  • Yali Chen
    Curious
    Collaborator

    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 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 CPtheta1

    CPtheta1=0.11

    global CPtheta2

    CPtheta2=0.68

    global CStheta1

    CStheta1=0.05

    global CStheta2

    CStheta2=0.62

    global H

    H=100

    global Pmax

    Pmax=10**(-2)  

    global N0w

    #N0w=10**(-11)

    N0w=0

    global beta0

    beta0=10**(-4)  

    global threshold

    threshold=10**(-3)

    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 = []

    for i in range(0,K):

        upos.append([400*random.random(),400*random.random()])

    #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    

    q=np.array([250,250])

    #create model

    model2 = gp.Model("MyModel2")

    #declare variables

    b=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 objective

    obj = 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=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+(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=2

    model2.params.NumericFocus=2

    model2.params.FeasibilityTol=1e-9

    model2.params.OptimalityTol=1e-9

    model2.params.ScaleFlag=0

    model2.params.Aggregate=0

    model2.optimize()

    print("Obj: ", model2.objVal)

    print("b: ", b)

    model2.printQuality()

    The output results are as follows,

    Set parameter NonConvex to value 2
    Set parameter NumericFocus to value 2
    Set parameter FeasibilityTol to value 1e-09
    Set parameter OptimalityTol to value 1e-09
    Set parameter ScaleFlag to value 0
    Set parameter Aggregate to value 0
    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, 362 columns and 0 nonzeros
    Model fingerprint: 0x48473f75
    Model has 352 quadratic constraints
    Coefficient 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 MIP

    Found heuristic solution: objective 0.3558289
    Presolve time: 0.00s
    Presolved: 25842 rows, 6662 columns, 61250 nonzeros
    Presolved 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%     -    0s

    Explored 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.355829

    Optimal solution found (tolerance 1.00e-04)
    Best objective -1.665334536938e-16, best bound -1.665334536938e-16, gap 0.0000%
    Obj:  -1.6653345369377348e-16
    b:  {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.

    Therefore, I would like to ask you two questions about the above post and three questions about this post. Looking forward to your reply! Many thanks!

    Best regards,

    Yali

     

    0
  • Yali Chen
    Curious
    Collaborator

    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

    0
  • Maliheh Aramon
    Gurobi Staff 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

    0
  • Yali Chen
    Curious
    Collaborator

    Hi, Maliheh,

    Thank you for your reply! Regarding your comments, I would like to ask a few questions.

    (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
    Warning: switch to quad precision"?

    Looking forward to your comments.

    Best regards,

    Yali

     

    0
  • Maliheh Aramon
    Gurobi Staff 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
    Warning: switch to quad precision"?

    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

    0
  • Yali Chen
    Curious
    Collaborator
    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 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=(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=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*(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=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]]*(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=2

    model1.params.NumericFocus=2

    model1.params.FeasibilityTol=1e-9

    model1.params.OptimalityTol=1e-9

    model1.params.ScaleFlag=2

    model1.params.Aggregate=0

    model1.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 threads

    Optimize a model with 0 rows, 384 columns and 0 nonzeros
    Model fingerprint: 0x1643e273
    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-07, 2e+05]
      RHS range        [0e+00, 0e+00]
      QRHS range       [1e+03, 2e+05]

    Continuous model is non-convex -- solving as a MIP

    Presolve time: 0.01s
    Presolved: 26138 rows, 6741 columns, 87620 nonzeros
    Presolved 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.
     

    Looking forward to your comments.

    Best regards,

    Yali

     
    0
  • Yali Chen
    Curious
    Collaborator

    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.

    Looking forward to your comments.

    Best regards,

    Yali

     
    0
  • Maliheh Aramon
    Gurobi Staff 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

    0

Please sign in to leave a comment.