Help with constraints
ユーザーの入力を待っています。In the given code, T1[i]== (TL[i] * (1-D[i])+ TR[i]*D[i]) constraint set the decision variable D[i] such that minimize T1
when i run the code it produces D[7] 1 D[8] 1 whereas TL[7] 3.71598 TL[8] 3.54706 and TR[7] 7.46835 TR[8] 7.1529
In my understanding, it should have produced T1[7] = 3.71598 , T1[8] 3.54706 and D[7] 0 D[8] 0 but it produces the T1[7] 7.46835 T1[8] 7.1529 and D[7] 1 D[8] 1 , here i am unable to find where is the mistake in the code, please help thank you. I have given working code below
import gurobipy as gp
from gurobipy import GRB
import numpy.random as nmp
import numpy as nmp1
import random
import math
P = 20
M = 6
R = 2
r_d = nmp1.array(
[
0,
3096.4327405932236,
1683.420450973956,
2781.2908051446766,
2740.9173654214164,
1789.9154088453154,
2504.217255624966,
2464.5070229415232,
1972.5534679018601,
2816.3383083746858,
2958.861446074776,
1833.8277583380036,
3305.358630760327,
2636.231647029585,
3312.1566452892807,
2346.8253565506584,
1994.681937809064,
3223.8572814704094,
2068.2525962614286,
1702.6147488906336,
1873.5679354356862,
]
)
x = nmp.uniform(0, 0, size=(P + 1))
a_d = nmp1.array(x)
for i in range(P + 1):
a_d[i] = r_d[i] / 5.0
# Computational intensity for each UE for player i
d = nmp1.array(
[
0,
0.9712271049819241,
0.8307421603424311,
0.8676264169227274,
0.5142155139168529,
0.9284516701363555,
0.6764654313048653,
0.7666229899262247,
0.8384578221643,
0.914441391384507,
0.698365237398614,
0.7711195295204585,
0.9706405098202289,
0.6749256027760206,
0.9226048032962242,
0.7897752688520887,
0.650164188957424,
0.8943314332418693,
0.6771326886174707,
0.5884800517308626,
0.9533348226987488,
]
)
# Computational intensity for each MEC host
d1 = nmp1.array(
[
0,
0.5165766471693557,
0.5737616270564242,
0.8939235367344246,
0.6579113772307663,
0.6561324790580341,
0.7590596765918923,
]
)
# Service rate for player i
CU = nmp1.array(
[
0,
697.0872334338458,
672.2448227405582,
947.9580302911515,
685.1722454187703,
621.8412692364695,
838.2480083345274,
844.4865709431585,
778.4426915632525,
537.8263447320076,
731.566226751524,
503.75277854718513,
529.3135692647196,
795.5989805996362,
946.0251485083222,
559.4746724838026,
580.2454108559145,
785.0229864868584,
659.0084532594252,
622.9189439555269,
922.8649075642064,
]
)
# Service rate of MEC host k
CM = nmp1.array([0, 25000.0, 25000.0, 25000.0, 25000.0, 25000.0, 25000.0])
# it keeps the information about regions of MEC host for e.g., CR[k]=j says mec host k belongs to regions j
CR = nmp1.array([0, 2, 2, 2, 1, 2, 1])
# Data Rate between player i and MEC host k
bw_initial = nmp1.array(
[
[0, 0, 0, 0, 0, 0, 0],
[0, 235.702260, 34.401046, 20.721217, 19.532889, 138.675049, 17.966327],
[0, 18.550354, 26.948403, 29.000740, 19.822394, 21.417647, 38.749213],
[0, 43.113061, 29.297937, 35.158128, 16.998054, 58.520574, 18.837771],
[0, 46.574643, 116.247639, 22.530295, 30.042088, 61.429512, 30.700278],
[0, 36.810509, 447.213595, 18.477377, 43.810795, 41.099747, 34.964180],
[0, 18.908356, 37.635496, 16.286559, 44.721360, 20.615270, 164.398987],
[0, 17.652863, 32.376195, 16.755332, 34.921515, 19.356450, 333.333333],
[0, 277.350098, 35.623525, 20.761370, 19.912577, 156.173762, 18.278756],
[0, 57.928445, 108.465229, 19.885000, 31.295877, 68.680282, 26.938623],
[0, 30.289127, 43.314808, 34.001020, 22.310537, 39.043440, 31.465839],
[0, 22.727273, 53.528773, 15.135887, 117.041147, 24.042352, 47.891314],
[0, 19.048483, 21.693046, 71.247050, 15.205718, 22.260731, 22.054055],
[0, 24.325213, 21.632053, 89.442719, 14.284257, 28.964222, 17.541160],
[0, 22.518867, 22.815832, 124.034735, 15.101330, 26.938623, 19.952172],
[0, 14.278431, 20.104818, 20.370021, 17.811762, 15.841118, 35.691531],
[0, 15.357378, 21.452116, 23.749017, 17.831574, 17.246507, 34.921515],
[0, 18.027531, 32.108065, 18.784519, 29.361011, 20.092640, 164.398987],
[0, 19.532889, 21.081851, 110.431526, 14.648968, 22.863490, 20.323692],
[0, 71.247050, 35.136418, 26.938623, 18.938850, 124.034735, 19.473541],
[0, 20.395425, 44.324221, 16.228549, 55.470020, 22.173004, 98.058068],
]
)
# Latency between MEC hosts
latency_mec_host = nmp1.array(
[
[0, 0, 0, 0, 0, 0, 0],
[0, 0.000000, 2.849748, 3.956402, 3.804876, 2.000000, 4.083115],
[0, 2.849748, 0.000000, 3.994486, 2.769632, 2.697032, 2.931996],
[0, 3.956402, 3.994486, 0.000000, 5.000000, 3.639770, 4.231276],
[0, 3.804876, 2.769632, 5.000000, 0.000000, 3.751287, 2.925613],
[0, 2.000000, 2.697032, 3.639770, 3.751287, 0.000000, 3.853929],
[0, 4.083115, 2.931996, 4.231276, 2.925613, 3.853929, 0.000000],
]
)
# Bandwidth between MEC hosts
bw_mec_host = nmp1.array(
[
[0, 0, 0, 0, 0, 0, 0],
[
0,
0.0,
3472.565599444328,
1704.7221178954644,
4104.194795577694,
3603.923997986016,
6215.600083555815,
],
[
0,
2329.8664567380893,
0.0,
9341.183260142236,
8403.785612711014,
7139.612692187351,
9486.888684590038,
],
[
0,
5239.108150849817,
6494.508770952193,
0.0,
5568.454323146592,
7238.0277643526515,
1422.7675824017113,
],
[
0,
1201.9849626813082,
7122.958215026663,
2195.564995976808,
0.0,
1082.6887982980793,
3030.6003011988014,
],
[
0,
3318.444020573778,
7787.940892379947,
6392.915740494214,
3724.520614146774,
0.0,
5796.001194494097,
],
[
0,
7400.059939491781,
8872.406342836044,
2295.2523882268415,
1873.9010866628894,
9256.691445957442,
0.0,
],
]
)
# intial region coordinators region 1 has mec host 2, and region 2 has mec host 1
E_initial = nmp1.array([0, 6, 2])
# The MEC host to which player i is conne
C_initial = nmp1.array([0, 1, 6, 5, 2, 2, 6, 6, 1, 2, 2, 4, 3, 4, 3, 6, 6, 6, 3, 5, 6])
m = gp.Model()
# create optimize and auxiliary variables
E = m.addVars(range(1, R + 1), range(1, M + 1), lb=0, ub=1, vtype=GRB.BINARY, name="E")
TL = m.addVars(range(1, P + 1), vtype=GRB.CONTINUOUS, name="TL")
TL_temp = m.addVars(range(1, P + 1), vtype=GRB.CONTINUOUS, name="TL_temp")
C = m.addVars(range(1, P + 1), range(1, M + 1), lb=0, ub=1, vtype=GRB.BINARY, name="C")
r_t = m.addVar(vtype=GRB.CONTINUOUS, name="r_t")
bw = m.addVars(range(1, P + 1), range(1, M + 1), vtype=GRB.CONTINUOUS, name="bw")
TR = m.addVars(range(1, P + 1), vtype=GRB.CONTINUOUS, name="TR")
D = m.addVars(
range(1, P + 1), lb=0, ub=1, vtype=GRB.BINARY, name="D"
) # NO NEED TO DEFINE
T1 = m.addVars(range(1, P + 1), vtype=GRB.CONTINUOUS, name="T1")
T1_temp = m.addVars(range(1, P + 1), vtype=GRB.CONTINUOUS, name="T1_temp")
latency = m.addVars(
range(1, M + 1), range(1, M + 1), vtype=GRB.CONTINUOUS, name="latency"
)
b = m.addVars(range(1, P + 1), range(1, M + 1), vtype=GRB.BINARY, name="b")
for i in range(1, P + 1):
for j in range(1, M + 1):
m.addConstr((bw[i, j] == (1 / bw_initial[i][j])))
for i in range(1, M + 1):
for j in range(1, M + 1):
m.addConstr(latency[i, j] == latency_mec_host[i][j])
for i in range(1, P + 1):
m.addConstr(
TL_temp[i] == (sum(bw[i, k] * C[i, k] for k in range(1, M + 1)))
) # optimize3 C[i,k] decision variable,
for i in range(1, P + 1):
for k in range(1, M + 1):
m.addConstr(
(C[i, k] == 1) >> (TL[i] == (r_d[i] * d[i] / CU[i] + a_d[i] * TL_temp[i]))
) # optimize
for i in range(1, P + 1):
for k in range(1, M + 1):
m.addConstr(
(
(C[i, k] == 1)
>> (TR[i] == (r_d[i] * TL_temp[i] + r_d[i] * d1[k] / CM[k]))
)
)
for i in range(1, P + 1):
m.addConstr((T1[i] == (TL[i] * (1 - D[i]) + TR[i] * D[i]))) # Optimize model,
for i in range(1, P + 1):
m.addConstr((sum(C[i, k] for k in range(1, M + 1)) == 1)) # optimize model 2
for j in range(1, R + 1):
temp_list = []
for k in range(1, M + 1):
if CR[k] == j:
temp_list.append(k)
for i in range(1, P + 1):
for k in range(1, M + 1):
if CR[k] == j:
m.addConstr(
(C[i, k] == 1)
>> (
T1_temp[i]
== T1[i]
+ sum(latency_mec_host[k][k1] * E[j, k1] for k1 in temp_list)
)
)
for j in range(1, R + 1):
temp_list = []
for k in range(1, M + 1):
if CR[k] == j:
temp_list.append(k)
m.addConstr(sum(E[j, k] for k in temp_list) == 1)
m.addConstr((r_t == (gp.max_(T1_temp[j] for j in range(1, P + 1)))), name="r_tc")
m.setObjective(r_t, GRB.MINIMIZE)
m.optimize()
for v in m.getVars():
if v.X != 0:
print("%s %g" % (v.VarName, v.X)) # printing the gurobi variables
print("Obj: %g" % m.ObjVal) # printing the value of object,
optimal_value = m.ObjVal
-
-
Gurobi Staff
The following constraints look suspicious to me:
R163: C[8,1] + C[8,2] + C[8,3] + C[8,4] + C[8,5] + C[8,6] = 1
GC42: C[8,1] = 1 -> TL[8] - 394.510693580372 TL_temp[8]
= 2.124630242822753
GC43: C[8,2] = 1 -> TL[8] - 394.510693580372 TL_temp[8]
= 2.124630242822753
GC44: C[8,3] = 1 -> TL[8] - 394.510693580372 TL_temp[8]
= 2.124630242822753
GC45: C[8,4] = 1 -> TL[8] - 394.510693580372 TL_temp[8]
= 2.124630242822753
GC46: C[8,5] = 1 -> TL[8] - 394.510693580372 TL_temp[8]
= 2.124630242822753
GC47: C[8,6] = 1 -> TL[8] - 394.510693580372 TL_temp[8]
= 2.124630242822753By \(\texttt{R163}\), exactly one of the binary variables \(\texttt{C[8,1]}\), \(\texttt{C[8,2]}\), \(\texttt{C[8,3]}\), \(\texttt{C[8,4]}\), \(\texttt{C[8,5]}\), and \(\texttt{C[8,6]}\) must equal \( 1 \). However, no matter which of these binary variables equals one, constraints \(\texttt{GC42}\) through \(\texttt{GC47}\) enforce the exact same constraint:
TL[8] - 394.510693580372 TL_temp[8] = 2.124630242822753
In your Python code, these indicator constraints are built as follows:
for i in range(1,P+1):
for k in range(1,M+1):
m.addConstr((C[i,k]==1)>> (TL[i]==(r_d[i]*d[i]/CU[i] + a_d[i]*TL_temp[i])))Should the constraints enforced by these indicator constraints have some dependency on \(\texttt{k}\)?
1 -
-
Thank you!
m.addConstr((C[i,k]==1)>> (TL[i]==(r_d[i]*d[i]/CU[i] + a_d[i]*TL_temp[i])))this constraints have indirect dependency on k, through TL_temp[i] (m.addConstr(TL_temp[i]==(sum(bw[i,k]*C[i,k] for k in range(1,M+1))))
to make direct dependency this TL_temp[i] can be replaced with bw[i,k]. will this solve the problem? shall i also replace TL_temp[i] in TR constraints?, will it not affect the TL_temp constraints? Previous constraints m.addConstr(TL_temp[i]==(sum(bw[i,k]*C[i,k] for k in range(1,M+1)))) this assign the value of C[i,k]=1 for which k it has highest bandwidth (bw).
0 -
-
-
I can't say if those changes are valid, as I'm not familiar with your application. That said, maybe it would help to understand why there exists no feasible solution with \(\texttt{T1[8]}=3.54706\) and \(\texttt{D[8]}=0\). To do this, let's fix \(\texttt{T1[8]}\) to \(3.54706\) and \(\texttt{D[8]}\) to \(0\) and see where the infeasibility comes from.
First, consider the constraint \(\texttt{qc1}\):
qc1: - TL[8] + T1[8] + [ TL[8] * D[8] - TR[8] * D[8] ] = 0
With \(\texttt{T1[8]}\) and \(\texttt{D[8]}\) fixed to \(3.54706\) and \(0\) respectively, this constraint yields \(\texttt{TL[8]} = 3.54706\).
Next, recall constraints \(\texttt{R163}\) and \(\texttt{GC42}\) through \(\texttt{GC47}\) that we discussed previously. Together, these constraints enforce the following relationship:
TL[8] = 2.124630242822753 + 394.510693580372 TL_temp[8]
Finally, let's examine the following set of constraints:
R42: bw[8,1] = 0.00360555127692798
R43: bw[8,2] = 0.0280713376904728
R44: bw[8,3] = 0.0481663782303384
R45: bw[8,4] = 0.0502195170419178
R46: bw[8,5] = 0.0064031242328657
R47: bw[8,6] = 0.0547083182247194
R163: C[8,1] + C[8,2] + C[8,3] + C[8,4] + C[8,5] + C[8,6] = 1
qc0: TL_temp[8] + [ - C[8,1] * bw[8,1] - C[8,2] * bw[8,2]
- C[8,3] * bw[8,3] - C[8,4] * bw[8,4] - C[8,5] * bw[8,5]
- C[8,6] * bw[8,6] ] = 0Together, these constraints tell us that \(\texttt{TL_temp[8]}\) must equal one of the following values:
0.00360555127692798
0.0280713376904728
0.0481663782303384
0.0502195170419178
0.0064031242328657
0.0547083182247194Given the relationship \(\texttt{TL[8]} = 2.124630242822753 + 394.510693580372 \cdot \texttt{TL_temp[8]}\), it follows that \(\texttt{TL[8]}\) must equal one of the following values:
3.547058777823206
13.199073144820016
21.126781525728084
21.936766742301057
4.650731225011888
23.70764681027251This contradicts our conclusion that \(\texttt{TL[8]}\) equals \(3.54706\). The first value is very close, but it differs by ~1.22e-6, which is enough for Gurobi to reasonably claim that the model is infeasible.
Looking at your data, my best guess for the source of the infeasibility is the \(\texttt{bw_initial}\) and \(\texttt{latency_mec_host}\) arrays. All of the values included in those arrays are truncated to the sixth decimal point, which can easily cause calculation errors on the order of 1e-6 or higher.
0 -
-
okay thank you so much for detaild explainations
0 -
-
if i removed TL_temp constraints
and update the
TL and TR as
for i in range(1,P+1):m.addConstr(TL[i]==sum((raw_data[i]*delta[i]/CU[i] + action_data[i]*bw[i,k])*C[i,k] for k in range(1,M+1))) #optimizefor i in range(1,P+1):m.addConstr(TR[i]==sum((raw_data[i]*bw[i,k]+ raw_data[i] *delta1[k]/ CM[k])*C[i,k] for k in range(1,M+1)))Further, i truncate bw_initial and latency_bw_host upto 12 decimal points. still having same problem...0 -
-
-
I would try computing an irreducible infeasible subsystem (IIS) by calling Model.computeIIS() after optimizing. Then, write out an ILP model file and inspect it:
m.computeIIS()
m.write("model.ilp")The ILP file will contain a subset of constraints and bounds that together form an infeasible system. However, this infeasible system can be made feasible by removing any single constraint or bound. This should give you a good idea of why the model is still infeasible.
See the article How do I determine why my model is infeasible? for more information.
0 -
サインインしてコメントを残してください。
コメント
6件のコメント