Interpret the formulation of a pricing model in crowdshipping

Question

I am trying to run the pricing model from the paper "Designing pricing and compensation schemes by integrating matching and routing models for crowd-shipping systems" on python with Gurobi, but the constraints given in this model do not preserve the network flow as expected. Maybe there should be more constraints to secure the networkflow from the origin nodes to the pick up nodes. I would like to get some feedback about the model and my 'translation' to python, as I don't get why it does not work as expected. First the scenario: We have m drivers h in H, H=[1,..,m], each having an origin node tau[h] and a destination node tau_[h]. We have also n requests i in A, A=[1,..,n] (pickup nodes), which have to be delivered in the delivery nodes A_=[n+1,...,2n]. Resulting from this we have a Graph G(V,E) with V=A+A_+tau+tau_ nodes and E=V*V lines. A driver should drive from his origin passing a pickup and the appropriate delivery to his destination.(or directly from his origin to his destination).
This is the model from the paper:

The package size is not relevant and is neglected in the paper. That's why K is also not considered in the programm.

And here is my code:

import gurobipy as gp
from gurobipy import GRB, quicksum
def pricing(n, m, A, A_, H, WTP, ETP, R, d2, V, V1, E2, d, t, s, l, a, b, M):

    global mod
    mod = gp.Model("Pricing")

    # Decision variables
    x = mod.addVars(E, H, vtype=GRB.BINARY, name="x")
    S = mod.addVars(V, H, vtype=GRB.CONTINUOUS,lb=0, name="S")
    L = mod.addVars(V, H, vtype=GRB.CONTINUOUS,lb=0, name="L")
    z = mod.addVars(A, vtype=GRB.BINARY, name="z")
    p = mod.addVars(E, H, vtype=GRB.CONTINUOUS, lb=0, name="p")
    c = mod.addVars(E, H, vtype=GRB.CONTINUOUS, lb=0, name="c")
    mod.setObjective(quicksum(p[i,j,h]*d2[i]-c[i, j, h]*d[i, j] for i in A for j in A_ 
    for h in H if j!=i ), gp.GRB.MAXIMIZE) # (40)
    mod.addConstrs((p[i, j, h] <= WTP[i]/d2[i] + (1-x[i, j, h])*M for i in A for j in A_ 
    for h in H if i!=j), name='Price_not_greater_than_WTP (02)/(28)')
    mod.addConstrs((c[i, j, h] >= ETP[h]*x[i, j, h] for j in A_ for h in H for i in A if 
    i!=j), name='Compensation_not_less_than_ETP (03)/(29)')
    mod.addConstrs((p[i, j, h] >= c[i, j, h] for i in A for j in A_ for h in H), 
    name="constraint4")
    mod.addConstrs((quicksum(x[i, j, h] for j in A_ for h in H) + z[i] == 1 for i in A), 
    name="constraint5")
    mod.addConstrs((quicksum(x[i, j, h] for j in V if i!=j ) - quicksum(x[j, n + i, h] 
    for j in V if j!=n+i) == 0 for i in A for h in H), name="constraint6")
    mod.addConstrs((quicksum(x[tau[h], j, h] for j in pickup_destination[h] ) == 1 for h 
    in H), name="constraint7")
    mod.addConstrs((quicksum(x[i, tau_[h], h] for i in delivery_origin[h] ) == 1 for h 
    in H), name="constraint8")
    mod.addConstrs((quicksum(x[i, j, h] for i in V if i!=j) - quicksum(x[j, i, h] for i 
    in V if i!=j) == 0 for j in A_ for h in H), name="constraint9")
    mod.addConstrs((S[i, h] + s[i] + t[i, j] <= S[j, h] + (1-x[i, j, h])*M for (i, j) in 
    E for h in H), name='courier_follows_matched_paths (10)/(34)')
    mod.addConstrs(( S[i, h] == [a[i], b[i]] for i in V for h in H),name="constraint11")
    mod.addConstrs((S[i, h] <= S[n + i, h] for i in V[0:-n] for h in H), 
    name="constraint12")
    mod.addConstrs((L[i, h] + l[j] <= L[j, h]+ (1-x[i , j, h])*M for (i,j) in E for h in 
    H), name='capacity_for_next_loading (13)/(35)')
    mod.addConstrs((L[i, h] <= R[h] for i in V for h in H ), name="constraint14")
    mod.addConstrs((L[tau[h], h] == 0 for h in H), name="constraint15.1")
    mod.addConstrs((L[tau_[h], h] == 0 for h in H), name="constraint15.2")
    mod.addConstrs((p[i, j, h] <= M*x[i, j, h] for i, j in E for h in H if i!=j), 
    name="constraint23")
    mod.addConstrs((c[i, j, h] <= M*x[i, j, h] for i, j in E for h in H), 
    name="constraint24")
    mod.addConstrs((p[i, j, h]*x[i, j, h]==p[i, j, h] for i, j in E for h in H if i!=j), 
    name='constraint25')
    mod.addConstrs((c[i, j, h]*x[i, j, h]==c[i, j, h] for i, j in E for h in H 
    ),name='constraint26')
    mod.optimize()

    active_arcs=[(i, j, h) for i, j in E for h in H if x[i, j, h].x>0.99] 
    active_prices=[(i, j, h) for i, j in E for h in H if p[i, j, h].x>0.99]
    active_costs=[(i, j, h) for i, j in E for h in H if c[i, j, h].x>0.99]   
    active_load=[(i, h) for i in V for h in H  if L[i, h].x>0.99] 
    active_z=[i for i in A if z[i].x>0.99] 
    active_Start=[(i, h) for i in V for h in H if S[i, h].x>0.99]

    if mod.status == GRB.OPTIMAL:
          
        optimal_x = mod.getAttr("x", x)
        optimal_S = mod.getAttr("x", S)
        optimal_L = mod.getAttr("x", L)
        optimal_z = mod.getAttr("x", z)
        optimal_p = mod.getAttr("x", p)
        optimal_c = mod.getAttr("x", c)

         
        return active_arcs, active_prices, active_costs, active_load, 
        active_Start, active_z, optimal_L, optimal_S, optimal_c, optimal_p, optimal_x, 
        optimal_z
    else:
        return None
n = 2  # Number of packages
m = 2  # Number of couriers
# Sets/nodes
A = list(range(m+1, n+m + 1))  # Set of package pickup nodes
A_ = list(range(n+m + 1, 2 * n +m+ 1))  # Set of package delivery nodes
H = list(range(1, m + 1))  # Set of couriers
tau=[h for h in H]
tau.insert(0,0)
tau_=[2*n + m +h for h in H]
tau_.insert(0,0)
V= tau[1:] + A + A_ + tau_[1:]
destination={(h): tau_[h] for h in H}
pickup_destination= {h: A+[destination[h]] for h in H}
origin={(h): tau[h] for h in H}
delivery_origin={h: A_+[origin[h]] for h in H}
Vh={(h): [origin[h]]+A+A_+[destination[h]] for h in H }
# Edges
E=[(i, j) for i in V for j in V if i!=j ]
M = 10000  # Large positive number
# Other parameters
WTP = {i: 10 for i in A}  # Maximum price that a Sender is willing to pay
ETP = {h: 2 for h in H}  # Minimum compensation that courier h expects
R = {(h): 5 for h in H }  # Capacity of courier h for package size k
d = {(i, j): 1 for i, j in E}  # Travel distance for link (i, j)
for key,value in d.items():  
    if key[0]==key[1]:
        value= 0     
    d[key]=value
d1= {i:1 for i in V} #Travel distance for price relevant arcs 
d2= d1.copy()
for i in d2:
    if i not in A:
        d2[i]=0
t = {(i, j): 1 for i, j in E}  # Travel time for link (i, j)
s = {i: 1 for i in V}  # Service time at node i
for i in s:
    if i in tau:
        s[i]=0
    if i in tau_:
        s[i]=0
l = {i: 1 for i in V}  # Amount of packages that need to be loaded at node i
for key, value in l.items():
    global newvalue
    newvalue=value
    if key in A_:
        newvalue=value*(-1)
    if key in tau :
        newvalue=0
    if key in tau_:
        newvalue=0
    l[key]=newvalue #positive load in pickup and negative load at dropoff
a = {i:0 for i in V}  # Time window lower bound
b = {i: 10 for i in V}  # Time window upper bound


result = pricing(n, m, A, A_, H, WTP, ETP, R, d2, V, V1, E2, d, t, s, l, a, b, M)

When I run it like this, the plotted result looks like this: So there is no constraint (in my programm) forcing the driver actually to begin at his origin. It would be very helpful to get some feedback where the mistake could be.

Answer 1

For source/dest node you may add constraints

$ \sum_j x_{\tau_h,j}^{h}- \sum_i x_{i,\tau_h}^h = 1$

$ \sum_j x_{\tau'_h,j}^{h}- \sum_i x_{i,\tau'_h}^h = -1 \ \ \forall h$