Gurobi Python: request for checking the construction of variables

Hi there, I am using the Gurobi Python and trying to solve a MILP optimization problem as shown below:

where x_{s,n} is the decision variable, and R is a function of x_{s,n}. The code can normally run, however, the results are not logical. We want to make sure that the syntax and structure for our code are right (so we can focus on modifying the mathematical logic). We provide our fully-reproducible codes and results.

The code is a bit long, we divided it into 3 parts to make it trackable:

part 1: define the fixed values and matrix shapes

part 2: define all the intermediate variables related to x (we use x_vars to denote x)

part 3: define the constraints of the MILP

import gurobipy as gp
from gurobipy import GRB
import numpy as np
# Create a new model
m = gp.Model("mip1")
# ---------------------------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------------------------
# Part 1, define the fixed values and matrix shapes
S = 2
N = 8
shape = (S, N)
C_exp_c = 8960 * 10 ** 9
C_exp_e = 4480 * 10 ** 9
RB_s = 9
f_CPU_c = 2.5 * 10 ** 9
f_CPU_e = 2.2 * 10 ** 9
alpha_DL = np.array([32.583, 1.072, 0.03])
alpha_UL = np.array([35.545, 1.623, 0.086])
i_DL = 27
i_UL = 16
lambda_n_c = 0
lambda_n_e = 0
for k in range(3):
    lambda_n_c = lambda_n_c + (C_exp_c * RB_s / f_CPU_c) * (alpha_DL[k] * i_DL + alpha_UL[k] * i_UL)
    lambda_n_e = lambda_n_e + (C_exp_e * RB_s / f_CPU_e) * (alpha_DL[k] * i_DL + alpha_UL[k] * i_UL)

lambda_sn_c = np.ones(shape) * lambda_n_c
lambda_sn_e = np.ones(shape) * lambda_n_e
b = np.array([1, 3, 3, 3, 200, 500, 10 ** 4, 2 * 10 ** 3])
b_sn = np.array([b] * S)
f = np.array([1, 3, 3, 3, 200, 500, 10 ** 4, 2 * 10 ** 3])
f_sn = np.array([f] * S)
C_sn_c = np.zeros(shape)
C_sn_e = np.zeros(shape)
for s in range(S):
    for n in range(1, N - 1):  # don't consider head and tail
        C_sn_c[s, n] = lambda_sn_c[s, n] / min(b_sn[s, n], f_sn[s, n + 1])  # central
        C_sn_e[s, n] = lambda_sn_e[s, n] / min(b_sn[s, n], f_sn[s, n + 1])  # edge

d_ec = 100 * 10 ** 3
d_e = 100 * 10 ** 3
d_c = 150 * 10 ** 3
v = 200 / (1 * 10 ** (-6))

C_sn_c_plus = np.zeros(shape)
C_sn_e_plus = np.zeros(shape)
C_sn_c_minus = np.zeros(shape)
C_sn_e_minus = np.zeros(shape)
for s in range(S):
    for n in range(1, N - 1):
        C_sn_c_plus[s, n] = max(lambda_sn_c[s, n] / (f_sn[s, n + 1] - d_ec / v), C_sn_c[s, n])
        C_sn_e_plus[s, n] = max(lambda_sn_e[s, n] / (f_sn[s, n + 1] - d_ec / v), C_sn_e[s, n])
        C_sn_c_minus[s, n] = max(lambda_sn_c[s, n] / (b_sn[s, n] - d_ec / v), C_sn_c[s, n])
        C_sn_e_minus[s, n] = max(lambda_sn_e[s, n] / (b_sn[s, n] - d_ec / v), C_sn_e[s, n])

delta_C_sn_c_minus = np.zeros(shape)
delta_C_sn_e_minus = np.zeros(shape)
delta_C_sn_c_plus = np.zeros(shape)
delta_C_sn_e_plus = np.zeros(shape)
# end Part 1, defining the fixed values and matrix shapes
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# Part 2, define all the constraints related to x, which is denoted by x_vars
x_vars = m.addVars(S, N, vtype=GRB.BINARY, name='X')
aux1 = m.addVars(S, N, vtype=GRB.BINARY)
aux2 = m.addVars(S, N, vtype=GRB.BINARY)
aux3 = m.addVars(S, N, vtype=GRB.BINARY)
aux4 = m.addVars(S, N, vtype=GRB.BINARY)
delta_x_sn_c_minus = m.addVars(S, N, vtype=GRB.BINARY)
delta_x_sn_c_plus = m.addVars(S, N, vtype=GRB.BINARY)
delta_x_sn_e_minus = m.addVars(S, N, vtype=GRB.BINARY)
delta_x_sn_e_plus = m.addVars(S, N, vtype=GRB.BINARY)

R_c_sn = m.addVars(S, N, vtype=GRB.CONTINUOUS)
R_e_sn = m.addVars(S, N, vtype=GRB.CONTINUOUS)
aux5 = m.addVars(S, N, vtype=GRB.CONTINUOUS)
aux6 = m.addVars(S, N, vtype=GRB.CONTINUOUS)
aux7 = m.addVars(S, N, vtype=GRB.CONTINUOUS)
aux8 = m.addVars(S, N, vtype=GRB.CONTINUOUS)
aux9 = m.addVars(S, N, vtype=GRB.CONTINUOUS)
aux10 = m.addVars(S, N, vtype=GRB.CONTINUOUS)


for s in range(S):
    for n in range(1, N - 1):
        delta_C_sn_c_minus[s, n] = C_sn_c_minus[s, n] - C_sn_c[s, n]
        delta_C_sn_e_minus[s, n] = C_sn_e_minus[s, n] - C_sn_e[s, n]
        delta_C_sn_c_plus[s, n] = C_sn_c_plus[s, n] - C_sn_c[s, n]
        delta_C_sn_e_plus[s, n] = C_sn_e_plus[s, n] - C_sn_e[s, n]

        m.addConstr(aux1[s, n] == x_vars[s, n] - x_vars[s, n - 1])
        m.addConstr(delta_x_sn_c_minus[s, n] == gp.max_(aux1[s, n],constant=0))
        m.addConstr(aux2[s, n] == x_vars[s, n] - x_vars[s, n + 1])
        m.addConstr(delta_x_sn_c_plus[s, n] == gp.max_(aux2[s, n], constant=0))
        m.addConstr(aux3[s, n] == -x_vars[s, n] - x_vars[s, n - 1])
        m.addConstr((delta_x_sn_e_minus[s, n] == gp.max_(aux3[s, n], 0)))
        m.addConstr(aux4[s, n] == -x_vars[s, n] + x_vars[s, n + 1])
        m.addConstr((delta_x_sn_e_plus[s, n] == gp.max_(aux4[s, n], 0)))

        m.addConstr(aux5[s, n] == C_sn_c[s, n] * x_vars[s, n])
        m.addConstr(aux6[s, n] == gp.max_(delta_x_sn_c_minus[s, n], delta_x_sn_c_plus[s, n]))
        m.addConstr(aux7[s, n] == gp.max_(aux6[s, n]))
        m.addConstr(R_c_sn[s, n] == aux5[s, n] + aux7[s, n])

        m.addConstr(aux8[s, n] == C_sn_e[s, n] * (1 - x_vars[s, n]))
        m.addConstr(aux9[s, n] == gp.max_(delta_x_sn_e_minus[s, n],
                                          delta_x_sn_e_plus[s, n], constant=delta_C_sn_e_minus[s, n]))
        m.addConstr(aux10[s, n] == gp.max_(aux9[s, n], constant=delta_C_sn_e_plus[s, n]))
        m.addConstr(R_e_sn[s, n] == aux8[s, n] + aux10[s, n])
# end Part 2, defining all the constraints related to x
# -----------------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------------
# Part 3, Defining the constraints of MILP

# <= Constraint: Equation (8b)
m.addConstr((gp.quicksum(R_c_sn[s, n] for n in range(N) for s in range(S))) <= C_exp_c)

# <= Constraint: Equation (8c)
m.addConstr((gp.quicksum(R_e_sn[s, n] for n in range(N) for s in range(S))) <= C_exp_e)

# <= Constraint: Equation (8d)
aux11 = m.addVars(S, N, vtype=GRB.CONTINUOUS)
aux12 = m.addVars(S, N, vtype=GRB.CONTINUOUS)
for s in range(S):
    for n in range(N - 1):
        m.addConstr(aux11[s, n] == x_vars[s, n] - x_vars[s, n + 1])
        m.addConstr(aux12[s, n] == gp.abs_(aux11[s, n]))
        m.addConstr((d_ec * aux12[s, n]) <= v * f_sn[s, n])

# <= Constraint: Equation (8e)
aux13 = m.addVars(S, N, vtype=GRB.CONTINUOUS)
aux14 = m.addVars(S, N, vtype=GRB.CONTINUOUS)
for s in range(S):
    for n in range(1, N):
        m.addConstr(aux13[s, n] == x_vars[s, n] - x_vars[s, n - 1])
        m.addConstr(aux14[s, n] == gp.abs_(aux13[s, n]))
        m.addConstr((d_ec * aux14[s, n]) <= v * b_sn[s, n])

# <= Constraint: Equation (8f)
aux15 = m.addVars(S, vtype=GRB.CONTINUOUS)
aux16 = m.addVars(S, vtype=GRB.CONTINUOUS)
for s in range(S):
    m.addConstr(aux15[s] == d_c * x_vars[s, 0])
    m.addConstr(aux16[s] == d_e * (1 - x_vars[s, 0]))
    m.addConstr(aux15[s] + aux16[s] <= v * b_sn[s, 0])

# End Part 3, Defining the constraints of MILP
# -----------------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------------
# Part 4, optimization procedure
objective = gp.quicksum((R_c_sn[s, n] + R_e_sn[s, n]) for n in range(N) for s in range(S))
m.ModelSense = GRB.MINIMIZE
m.setObjective(objective)
m.update()
m.optimize()

for v in x_vars.values():
    print("{}: {}".format(v.varName, v.X))

the results are shown as:

since x is a binary variable, and we are sure that under this parameter setting, it is impossible for all variables x to be equal to 0. If you can help us check whether the syntax and structure for our code are right, then we know that there must be something wrong with our mathematical logic.

Thank you so much for your support!

正式なコメント

Simranjit Kaur

Gurobi Staff

2025年11月16日 23:40

This post is more than three years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum, or try Gurobot, our chatbot interface offering instant, expert-level support.

Jaromił Najman

2022年04月11日 09:02 編集済み

You model seems to be correct. Note that the default lower bound for continuous variables is \(0\). Thus, it might be necessary to set the lower bound of variables \(\texttt{aux11, aux13, aux15}\) to \(-1\).

You could write the model to an LP file, insert \(x_{s,n}=0\) by hand and check why every optimization variable can indeed be \(0\).

In part 2, you have many \(\max\) expressions not shown in your equation system. Maybe there is some issue hidden in there.

For example you have

C129 = MAX ( C65 , 0 )
C241 = MAX ( C113 , C129 , 1533.1883313898 )

with \(\texttt{C113,C129}\) being binary. Thus, you force \(\texttt{C241}\) to \(1533.188\).

Best regards,
Jaromił

Gurobi Python: request for checking the construction of variables

コメント

お探しのものを見つけられませんでしたか？