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!