KeyError: 0. sum1=sum1+b0[combination[km]][0]*math.sqrt(beta0/(H**2+a[combination[km]][0]))*cov(k,combination[km])

Answered

Yali Chen

August 25, 2023 14:32

We use Gurobi to sove the following problem.

Optimization variables are underlined in green.

#initialized variables

b0=math.sqrt(Pmax)*np.mat(ones((K,1)))

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

#create model

model = Model("MyModel")

#declare variables

a=model.addVars(K,1,lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name="a")

q=model.addVars(1,2,lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name="q")

#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

combination=[]

def dfs(last_idx, cnt):

if cnt == M:

print(combination)

sum1=0

for k in range(0,K):

for km in range(0,M):

sum1=sum1+b0[combination[km]][0]*math.sqrt(beta0/(H**2+a[combination[km]][0]))*cov(k,combination[km])

sum2=0

for km in range(0,M):

for km1 in range(0,M):

sum2=sum2+b0[combination[km]][0]*math.sqrt(beta0/(H**2+((q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2)))*b0[combination[km1]][0]*math.sqrt(beta0/(H**2+((q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2)))*cov(combination[km],combination[km1])

sum3=0

for km in range(0,M):

sum4=0

for km1 in range(0,M):

sum4=sum4+b0[combination[km1]][0]*math.sqrt(beta0/(H**2+((q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2)))*cov(combination[km],combination[km1])

sum3=sum3+(a[combination[km]][0]-((q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2))*(-b0[combination[km]][0]*math.sqrt(beta0))/2/(H**2+(q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2)**(3/2)*sum4/(sum2+N0w)

sum5=0

for km1 in range(0,M):

sum6=0

for km in range(0,M):

sum6=sum6+b0[combination[km]][0]*math.sqrt(beta0/(H**2+((q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2)))*cov(combination[km],combination[km1])

sum5=sum5+(a[combination[km1]][0]-((q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2))*(-b0[combination[km1]][0])*math.sqrt(beta0)/2/(H**2+((q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2))**(3/2)*sum6/(sum2+N0w)

sum7=2*math.log(sum1)-math.log(K**2)-(math.log(sum2+N0w)+sum3+sum5)

sumvalue=sumvalue+sum7

return

for idx in range(last_idx+1,K):

combination.append(idx)

dfs(idx, cnt+1)

combination.pop()

dfs(-1,0)

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

model.setObjective(obj, GRB.MAXIMIZE)

model.addConstr(a[k][0]==(q[0]-upos[k][0])**2+(q[1]-upos[k][1])**2 for k in range(0,K))

#optimize

model.optimize()

print("Obj: ", model.objVal)

opt=model.objVal

Comments

13 comments

Jaromił Najman
Gurobi Staff

August 28, 2023 07:25

Edited
Your code is not executable.

You should edit your post to use a code block and make sure that your code can be copy-pasted and executed resulting in the error you are trying to solve. Please refer to minimal reproducible example for more details.
0

Yali Chen

August 28, 2023 14:47

from gurobipy import *

from numpy import *

import numpy as np

import math

import random

#define global variables

global K

K = 5

global M

M = 2

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**(-8)*10**(-3)  

global beta0

beta0=10**(-4)  

global threshold

threshold=10**(-3)

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

#initialized variables

b0=math.sqrt(Pmax)*np.mat(ones((K,1)))

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

#create model

model = Model("MyModel")

#declare variables

a=model.addVars(K,1,lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name="a")

q=model.addVars(1,2,lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name="q")

#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  

combination=[]

def dfs(last_idx, cnt):

    if cnt == M:

        print(combination)

        sum1=0

        for k in range(0,K):

            for km in range(0,M):

                sum1=sum1+b0[combination[km]][0]*math.sqrt(beta0/(H**2+a[combination[km]][0]))*cov(k,combination[km])

        sum2=0

        for km in range(0,M):

            for km1 in range(0,M):

                sum2=sum2+b0[combination[km]][0]*math.sqrt(beta0/(H**2+((q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2)))*b0[combination[km1]][0]*math.sqrt(beta0/(H**2+((q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2)))*cov(combination[km],combination[km1])  

        sum3=0

        for km in range(0,M):

            sum4=0

            for km1 in range(0,M):

                sum4=sum4+b0[combination[km1]][0]*math.sqrt(beta0/(H**2+((q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2)))*cov(combination[km],combination[km1])

            sum3=sum3+(a[combination[km]][0]-((q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2))*(-b0[combination[km]][0]*math.sqrt(beta0))/2/(H**2+(q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2)**(3/2)*sum4/(sum2+N0w)

        sum5=0

        for km1 in range(0,M):

            sum6=0

            for km in range(0,M):

                sum6=sum6+b0[combination[km]][0]*math.sqrt(beta0/(H**2+((q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2)))*cov(combination[km],combination[km1])

            sum5=sum5+(a[combination[km1]][0]-((q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2))*(-b0[combination[km1]][0])*math.sqrt(beta0)/2/(H**2+((q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2))**(3/2)*sum6/(sum2+N0w)      

        sum7=2*math.log(sum1)-math.log(K**2)-(math.log(sum2+N0w)+sum3+sum5)

        sumvalue=sumvalue+sum7

        return

    for idx in range(last_idx+1,K):

        combination.append(idx)

        dfs(idx, cnt+1)

        combination.pop()

dfs(-1,0)

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

model.setObjective(obj, GRB.MAXIMIZE)

model.addConstr(a[k][0]==(q[0]-upos[k][0])**2+(q[1]-upos[k][1])**2 for k in range(0,K))

#optimize

model.optimize()

print("Obj: ", model.objVal)

opt=model.objVal

Yali Chen

August 28, 2023 14:54
Hi, thank you very much for your comments. For the optimization problem given by the above picture, optimization variables are underlined in green, and we use Gurobi to sove the problem. According to your comments, we have given the complete code, looking forward to your reply, many thanks!

All the best,

Yali
0
Jaromił Najman
Gurobi Staff

August 29, 2023 07:05
Hi Yali,

\(\texttt{a}\) is a tupledict of Var objects returned by the addVars method and cannot be accessed via two indices but rather by a tuple, i.e., you have to access the entries of \(\texttt{a}\) via \(\texttt{a[i,j]}\) instead of \(\texttt{a[i][j]}\). The same applies to your variables \(\texttt{q}\).

You will also run into the error
```
  File "test3.py", line 72, in dfs
    sum1=sum1+b0[combination[km]][0]*math.sqrt(beta0/(H**2+a[combination[km],0]))*cov(k,combination[km])
  File "src/gurobipy/linexpr.pxi", line 570, in gurobipy.LinExpr.__truediv__
gurobipy.GurobiError: Divisor must be a constant
```
Gurobi currently does not support division by a variable. However, we discuss a workaround in the Knowledge Base article How do I divide by a variable in Gurobi?

Best regards,
Jaromił
0
Yali Chen

August 29, 2023 12:32
Hi Jaromił, thanks for your help!

All the best,

Yali
0

Yali Chen

August 30, 2023 11:40

Hi Jaromił, according to your comments, we have changed the code, but there are still some problems. Looking forward to your reply! Many thanks!

We use Gurobi to sove the following problem. Optimization variables are underlined in green.

Then, I will attach my complete code below. sum1=sum1+b0[combination[km]]*math.sqrt(v[combination[km]])*cov(k,combination[km])
TypeError: must be real number, not Var

from gurobipy import *

from numpy import *

import numpy as np

import math

import random

#define global variables

global K

K = 5

global M

M = 2

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**(-8)*10**(-3)  

global beta0

beta0=10**(-4)  

global threshold

threshold=10**(-3)

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

#initialized variables

b0=math.sqrt(Pmax)*np.mat(ones((K)))

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

#create model

model = Model("MyModel")

#declare variables

q=model.addVars(2, lb=0, ub=400, vtype=GRB.CONTINUOUS, name="q")

a=model.addVars(K, name="a")

v=model.addVars(K, name="v")

sumvariable=model.addVars(1, name="sumvariable")

#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  

combination=[]

def dfs(last_idx, cnt):

    if cnt == M:

        print(combination)

        sum1=0

        for k in range(0,K):

            for km in range(0,M):

                sum1=sum1+b0[combination[km]]*math.sqrt(v[combination[km]])*cov(k,combination[km])

        #Since lg(sum1) is a component of the objective function, we assign sum1 to a defined variable, that is sumvariable

        sumvariable=sum1        

        sum2=0

        for km in range(0,M):

            for km1 in range(0,M):

                sum2=sum2+b0[combination[km]]*math.sqrt(beta0/(H**2+(q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2))*b0[combination[km1]]*math.sqrt(beta0/(H**2+(q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2))*cov(combination[km],combination[km1])  

        sum3=0

        for km in range(0,M):

            sum4=0

            for km1 in range(0,M):

                sum4=sum4+b0[combination[km1]]*math.sqrt(beta0/(H**2+(q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2))*cov(combination[km],combination[km1])

            sum3=sum3+(v[combination[km]]-(beta0/(H**2+(q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2)))*1/2*b0[combination[km]]*(beta0/(H**2+(q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2))**(-1/2)*sum4/(sum2+N0w)        

        sum5=0

        for km1 in range(0,M):

            sum6=0

            for km in range(0,M):

                sum6=sum6+b0[combination[km]]*math.sqrt(beta0/(H**2+(q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2))*cov(combination[km],combination[km1])

            sum5=sum5+(v[combination[km1]]-(beta0/(H**2+(q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2)))*1/2*b0[combination[km1]]*(beta0/(H**2+(q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2))**(-1/2)*sum6/(sum2+N0w)      

        sum7=2*math.log(sumvariable)-math.log(K**2)-(math.log(sum2+N0w)+sum3+sum5)

        sumvalue=sumvalue+sum7

        return

    for idx in range(last_idx+1,K):

        combination.append(idx)

        dfs(idx, cnt+1)

        combination.pop()

dfs(-1,0)

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

model.setObjective(obj, GRB.MAXIMIZE)

model.addConstr(a[k]==(q[0]-upos[k][0])**2+(q[1]-upos[k][1])**2 for k in range(0,K))

model.addConstr(v[k]*(H**2+a[k])/beta0 == 1 for k in range(0,K))

model.params.NonConvex=2

#optimize

model.optimize()

print("Obj: ", model.objVal)

opt=model.objVal

Jaromił Najman
Gurobi Staff

August 30, 2023 12:31
Hi Yali,

you are calling the \(\sqrt{}\) function from the \(\texttt{math}\) package on a Gurobi Var. This is not possible. You have to use Gurobi's general constraints feature to implement nonlinear functions which have optimization variables or expressions holding optimization variables as input. Note that Gurobi currently uses a piecewise-linear approximation of the nonlinear functions such that the final solution will very likely suffer from some (small) approximation error.

In order to model the particular function \(\sqrt{f(x)}\) where \(x\) is an optimization variable, you have to introduce an auxiliary variable for the argument variable and for the result. This means that you have to model

\[\begin{align*}
w_1 & = f(x)\\
w_2 & = \sqrt{w_1}
\end{align*}\]

with appropriate bounds for the auxiliary variables \(w_1, w_2\). You can then use \(w_2\) in the rest of your model. In the particular case of the \(\sqrt{}\) function, you can use the addGenConstrPow method with \(a=0.5\). You can find the list of all supported general functions here.

Best regards,
Jaromił
0

Yali Chen

August 31, 2023 06:09

Hi Jaromił, according to your comments, we have changed the code, but there are still some problems. TypeError: unsupported operand type(s) for *: 'int' and 'tupledict'：sum7=2*nu-math.log(K**2)-(math.log(sum2+N0w)+sum3+sum5). Looking forward to your reply! Many thanks!

We use Gurobi to sove the following problem. Optimization variables are underlined in green.

from gurobipy import *

from numpy import *

import numpy as np

import math

import random

#define global variables

global K

K = 5

global M

M = 2

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**(-8)*10**(-3)  

global beta0

beta0=10**(-4)  

global threshold

threshold=10**(-3)

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

#initialized variables

b0=math.sqrt(Pmax)*np.array(ones((K)))

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

#create model

model = Model("MyModel")

#declare variables

q=model.addVars(2, lb=0, ub=400, vtype=GRB.CONTINUOUS, name="q")

a=model.addVars(K, name="a")

v=model.addVars(K, name="v")

r=model.addVars(K, name="r")

varpi=model.addVars(1, name="varpi")

nu=model.addVars(1, 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  

combination=[]

def dfs(last_idx, cnt):

    if cnt == M:

        print(combination)    

        sum2=0

        for km in range(0,M):

            for km1 in range(0,M):

                sum2=sum2+b0[combination[km]]*math.sqrt(beta0/(H**2+(q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2))*b0[combination[km1]]*math.sqrt(beta0/(H**2+(q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2))*cov(combination[km],combination[km1])  

        sum3=0

        for km in range(0,M):

            sum4=0

            for km1 in range(0,M):

                sum4=sum4+b0[combination[km1]]*math.sqrt(beta0/(H**2+(q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2))*cov(combination[km],combination[km1])

            sum3=sum3+(v[combination[km]]-(beta0/(H**2+(q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2)))*1/2*b0[combination[km]]*(beta0/(H**2+(q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2))**(-1/2)*sum4/(sum2+N0w)        

        sum5=0

        for km1 in range(0,M):

            sum6=0

            for km in range(0,M):

                sum6=sum6+b0[combination[km]]*math.sqrt(beta0/(H**2+(q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2))*cov(combination[km],combination[km1])

            sum5=sum5+(v[combination[km1]]-(beta0/(H**2+(q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2)))*1/2*b0[combination[km1]]*(beta0/(H**2+(q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2))**(-1/2)*sum6/(sum2+N0w)      

        sum7=2*nu-math.log(K**2)-(math.log(sum2+N0w)+sum3+sum5)

        sumvalue=sumvalue+sum7

        return

    for idx in range(last_idx+1,K):

        combination.append(idx)

        dfs(idx, cnt+1)

        combination.pop()

dfs(-1,0)

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

model.setObjective(obj, GRB.MAXIMIZE)

model.addConstr(a[k]==(q[0]-upos[k][0])**2+(q[1]-upos[k][1])**2 for k in range(0,K))

model.addConstr(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)

varpi=LinExpr(0)

for k in range(0,K):

    for km in range(0,M):

        varpi.addTerms(1, b0[combination[km]]*r[combination[km]]*cov(k,combination[km]))

model.addGenConstrLogA(varpi, nu, 10.0)

model.params.NonConvex=2

#optimize

model.optimize()

print("Obj: ", model.objVal)

opt=model.objVal

Jaromił Najman
Gurobi Staff

August 31, 2023 07:22
The errror is caused because the \(\texttt{nu}\) is a tupledict of Var objects returned by the addVars method and you have to access its indices. In this case, this would be \(\texttt{nu[0]}\). It holds only 1 variable so it is probably better to just use the addVar method, then you don't have to change the code for \(\texttt{sum7}\). The same holds for the optimization variable \(\texttt{varpi}\).
0

Yali Chen

September 01, 2023 13:58

Hi Jaromił, sorry to bother you again. IndexError: list index out of range: constr_1.addTerms(1, b0[combination[km]]*r[combination[km]]*cov(k,combination[km])).

from gurobipy import *

from numpy import *

import numpy as np

import math

import random

#define global variables

global K

K = 5

global M

M = 2

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**(-8)*10**(-3)  

global beta0

beta0=10**(-4)  

global threshold

threshold=10**(-3)

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

#initialized variables

b0=math.sqrt(Pmax)*np.array(ones((K)))

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

#create model

model = Model("MyModel")

#declare variables

q=model.addVars(2, lb=0, ub=400, vtype=GRB.CONTINUOUS, name="q")

a=model.addVars(K, name="a")

v=model.addVars(K, name="v")

r=model.addVars(K, name="r")

varpi=model.addVars(1, name="varpi")

nu=model.addVars(1, 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  

combination=[]

def dfs(last_idx, cnt):

    global sumvalue

    if cnt == M:

        print(combination)    

        sum2=0

        for km in range(0,M):

            for km1 in range(0,M):

                sum2=sum2+b0[combination[km]]*math.sqrt(beta0/(H**2+(q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2))*b0[combination[km1]]*math.sqrt(beta0/(H**2+(q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2))*cov(combination[km],combination[km1])  

        sum3=0

        for km in range(0,M):

            sum4=0

            for km1 in range(0,M):

                sum4=sum4+b0[combination[km1]]*math.sqrt(beta0/(H**2+(q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2))*cov(combination[km],combination[km1])

            sum3=sum3+(v[combination[km]]-(beta0/(H**2+(q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2)))*1/2*b0[combination[km]]*(beta0/(H**2+(q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2))**(-1/2)*sum4/(sum2+N0w)        

        sum5=0

        for km1 in range(0,M):

            sum6=0

            for km in range(0,M):

                sum6=sum6+b0[combination[km]]*math.sqrt(beta0/(H**2+(q0[0]-upos[combination[km]][0])**2+(q0[1]-upos[combination[km]][1])**2))*cov(combination[km],combination[km1])

            sum5=sum5+(v[combination[km1]]-(beta0/(H**2+(q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2)))*1/2*b0[combination[km1]]*(beta0/(H**2+(q0[0]-upos[combination[km1]][0])**2+(q0[1]-upos[combination[km1]][1])**2))**(-1/2)*sum6/(sum2+N0w)      

        sum7=2*nu[0]-math.log(K**2)-(math.log(sum2+N0w)+sum3+sum5)

        sumvalue=sumvalue+sum7

        return

    for idx in range(last_idx+1,K):

        combination.append(idx)

        dfs(idx, cnt+1)

        combination.pop()

dfs(-1,0)

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

model.setObjective(obj, GRB.MAXIMIZE)

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)

constr_1 = LinExpr(0)

for k in range(0,K):

    for km in range(0,M):

        constr_1.addTerms(1, b0[combination[km]]*r[combination[km]]*cov(k,combination[km]))  

model.addConstr(varpi[0] == constr_1)

model.addGenConstrLogA(varpi[0], nu[0], 10.0)

model.params.NonConvex=2

#optimize

model.optimize()

print("Obj: ", model.objVal)

opt=model.objVal

Thank you very much for your help and look forward to your reply!

All the best,

Yali

Eli Towle
Gurobi Staff

September 01, 2023 17:14
At that line of code, \(\texttt{combination}\) is an empty list, so you cannot access the list at index \(\texttt{km}\) (here, \(\texttt{km}\) is \(\texttt{0}\)). If you're expecting the \(\texttt{combination}\) list to be nonempty at this point, there must be an issue in how you construct \(\texttt{combination}\).

Your \(\texttt{dfs}\) function:
```
def dfs(last_idx, cnt):
    ...
    for idx in range(last_idx+1,K):
        combination.append(idx)
        dfs(idx, cnt+1)
        combination.pop()
```
Here, elements are being added to \(\texttt{combination}\), the function is called recursively, and then the last element is removed from \(\texttt{combination}\).
0
Yali Chen

September 02, 2023 15:14
Hi Eli, sorry to bother you again. According to your comments about the code I showed above，the error is in the following part，
constr_1 = LinExpr(0)

for k in range(0,K):

for km in range(0,M):

constr_1.addTerms(1, b0[combination[km]]*r[combination[km]]*cov(k,combination[km]))

model.addConstr(varpi[0] == constr_1)

First, let me introduce you to the details. We 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. Now，if the constraint corresponding to the error code is changed as follows:
Optimization variables are underlined in red. For each combination, we need to write this constraint, so please tell me how should I add it to the my dfs function? Where it should be added is before the optimization objective
```
model.setObjective(obj, GRB.MAXIMIZE)
```
appears, but it is a constraint so I don't understand how to add it.
In addition, you also know that M numbers are selected from 1 to K. If K=50, the number of combinations will be very large, and the computational complexity of this code will be very high, which is why we use depth-first search. But in this constraint, for each combination, the equation must be satisfied, which also involves complexity issues, we cannot directly list the constraints for each combination. I really want to ask you these questions, looking forward to your reply！Many thanks！

All the best，

Yali
0
Eli Towle
Gurobi Staff

September 05, 2023 19:31
I see you created a follow-up post The nonlinear function in the optimization objective contains optimization variables with this same question. We can continue the discussion there.
0

Please sign in to leave a comment.

Comments

Didn't find what you were looking for?