對於Vrptw問題來說，數學模型主要由以下部分組成。首先我們定義一些相關引數，一個圖可以表示為 G ( V , A ) G(V,A)

接下來構建數學模型。
目標函式為最小化運輸成本：
m i n ∑ k ∈ K ∑ ( i , j ) ∈ A c i j x i j (1) min \sum_{k\in K} \sum_{(i,j)\in A}c_{ij}x_{ij} \tag{1} mink∈K∑​(i,j)∈A∑​cij​xij​(1)
約束一讓車輛駛出倉庫(depot)：
∑ ( 0 , j ) ∈ A x 0 j k = 1 , ∀ k ∈ K (2) \sum_{(0,j)\in A}x_{0j}^k =1,\ \forall k\in K \tag{2} (0,j)∈A∑​x0jk​=1,∀k∈K(2)
約束二為流平衡(除去depot之外的點)：
∑ ( i , j ) ∈ A x i j k − ∑ ( j , i ) ∈ A x j i k = 0 , ∀ k ∈ K , i ∈ V ∖ { s , t } (3) \sum_{(i,j)\in A}x_{ij}^k - \sum_{(j,i)\in A}x_{ji}^k = 0,\ \forall k\in K ,i \in V\setminus \{s,t\} \tag{3} (i,j)∈A∑​xijk​−(j,i)∈A∑​xjik​=0,∀k∈K,i∈V∖{s,t}(3)
約束三讓車輛駛回depot：
∑ ( i , n + 1 ) ∈ A x i , n + 1 k = 1 , ∀ k ∈ K (4) \sum_{(i,n+1)\in A}x_{i,n+1}^k =1,\ \forall k\in K \tag{4} (i,n+1)∈A∑​xi,n+1k​=1,∀k∈K(4)
約束四保證每個客戶點都被服務：
∑ k ∈ K ∑ ( i , j ) ∈ A x i , j k = 1 , ∀ i ∈ V ∖ { s , t } (5) \sum_{k \in K}\sum_{(i,j)\in A}x_{i,j}^k =1,\ \forall i \in V\setminus \{s,t\} \tag{5} k∈K∑​(i,j)∈A∑​xi,jk​=1,∀i∈V∖{s,t}(5)
約束五保證被服務的相鄰節點開始服務時間的大小關係(去迴路)：
s i + t i j + s e r v i − M ( 1 − x i j ) ≤ s j , ∀ ( i , j ) ∈ A (6) s_i+t_{ij}+serv_i-M(1-x_{ij})\le s_j,\ \forall (i,j) \in A \tag{6} si​+tij​+servi​−M(1−xij​)≤sj​,∀(i,j)∈A(6)
約束六保證不違反客戶的時間窗：
e i ≤ s i ≤ l i , ∀ i ∈ V ∖ { s , t } (7) e_i \le s_i \le l_i ,\ \forall i \in V\setminus \{s,t\} \tag{7} ei​≤si​≤li​,∀i∈V∖{s,t}(7)
約束七保證不違反車輛的載重約束：
∑ ( i , j ) ∈ A x i j k q i ≤ Q k , ∀ k ∈ K (8) \sum_{(i,j)\in A}x_{ij}^kq_i \le Q_k ,\ \forall k \in K \tag{8} (i,j)∈A∑​xijk​qi​≤Qk​,∀k∈K(8)
最後是決策變數的取值約束：
x i j k ∈ { 0 , 1 } ∀ ( i , j ) ∈ A , k ∈ K s i ≥ 0 , ∀ i ∈ V ∖ { s , t } (9) x_{ij}^k\in \{0,1\} \ \forall (i,j) \in A,k \in K \\ s_i \ge 0,\ \forall i \in V\setminus \{s,t\} \tag{9} xijk​∈{0,1}∀(i,j)∈A,k∈Ksi​≥0,∀i∈V∖{s,t}(9)
那麼(1)~(9)式就組成了Vrptw問題的數學模型。

class Data:
    customerNum = 0
    nodeNum     = 0
    vehicleNum  = 0
    capacity    = 0
    cor_X       = []
    cor_Y       = []
    demand      = []
    serviceTime = []
    readyTime   = []
    dueTime     = []
    disMatrix   = [[]]   

# function to read data from .txt files
def readData(data, path, customerNum):
    data.customerNum = customerNum
    data.nodeNum = customerNum + 2
    f = open(path, 'r')
    lines = f.readlines()
    count = 0
    # read the info
    for line in lines:
        count = count + 1
        if (count == 5):
            line = line[:-1].strip()
            str = re.split(r" +", line)
            data.vehicleNum = int(str[0])
            data.capacity = float(str[1])
        elif (count >= 10 and count <= 10 + customerNum):
            line = line[:-1]
            str = re.split(r" +", line)
            data.cor_X.append(float(str[2]))
            data.cor_Y.append(float(str[3]))
            data.demand.append(float(str[4]))
            data.readyTime.append(float(str[5]))
            data.dueTime.append(float(str[6]))
            data.serviceTime.append(float(str[7]))

    data.cor_X.append(data.cor_X[0])
    data.cor_Y.append(data.cor_Y[0])
    data.demand.append(data.demand[0])
    data.readyTime.append(data.readyTime[0])
    data.dueTime.append(data.dueTime[0])
    data.serviceTime.append(data.serviceTime[0])

    # compute the distance matrix
    data.disMatrix = [([0] * data.nodeNum) for p in range(data.nodeNum)]  # 初始化距離矩陣的維度,防止淺拷貝
    # data.disMatrix = [[0] * nodeNum] * nodeNum]; 這個是淺拷貝，容易重複
    for i in range(0, data.nodeNum):
        for j in range(0, data.nodeNum):
            temp = (data.cor_X[i] - data.cor_X[j]) ** 2 + (data.cor_Y[i] - data.cor_Y[j]) ** 2
            data.disMatrix[i][j] = math.sqrt(temp)
            temp = 0
            
    return data

data = Data()
path = 'c101.txt' //讀取Solomon資料集
customerNum = 20  //設定客戶數量
readData(data, path, customerNum)
data.vehicleNum = 10  //設定車輛數
printData(data, customerNum)
BigM = 100000 //定義一個極大值

x = {}
s = {}  //定義字典用來存放決策變數

for i in range(data.nodeNum):
    for k in range(data.vehicleNum):
        name = 's_' + str(i) + '_' + str(k)
        s[i,k] = model.addVar(0
                              , 1500
                              , vtype= GRB.CONTINUOUS
                              , name= name)  //定義訪問時間為連續變數
        for j in range(data.nodeNum):
            if(i != j):
                name = 'x_' + str(i) + '_' + str(j) + '_' + str(k)
                x[i,j,k] = model.addVar(0
                                        , 1
                                        , vtype= GRB.BINARY
                                        , name= name)  //定義是否服務為0-1變數

//首先定義一個線性表示式
obj = LinExpr(0) 
for i in range(data.nodeNum):
    for k in range(data.vehicleNum):
        for j in range(data.nodeNum):
            if(i != j):
           		 //將目標函式係數與決策變數相乘，並進行連加
                obj.addTerms(data.disMatrix[i][j], x[i,j,k]) 
//將表示目標函式的線性表示式加入模型，並定義為求解最小化問題
model.setObjective(obj, GRB.MINIMIZE)  

約束一（vehicle_depart）：

for k in range(data.vehicleNum):
	//同樣先定義一個線性表示式
    lhs = LinExpr(0) 
    for j in range(data.nodeNum):
        if(j != 0):
       		//約束係數與決策變數相乘	
            lhs.addTerms(1, x[0,j,k])  
    //將約束加入模型
    model.addConstr(lhs == 1, name= 'vehicle_depart_' + str(k)) 

約束二（flow_conservation）：

for k in range(data.vehicleNum):
    for h in range(1, data.nodeNum - 1):
        expr1 = LinExpr(0)
        expr2 = LinExpr(0)
        for i in range(data.nodeNum):
            if (h != i):
                expr1.addTerms(1, x[i,h,k])

        for j in range(data.nodeNum):
            if (h != j):
                expr2.addTerms(1, x[h,j,k])

        model.addConstr(expr1 == expr2, name= 'flow_conservation_' + str(i))
        expr1.clear()
        expr2.clear()

約束三（vehicle_enter）：

for k in range(data.vehicleNum):
    lhs = LinExpr(0)
    for j in range(data.nodeNum - 1):
        if(j != 0):
            lhs.addTerms(1, x[j, data.nodeNum-1, k])
    model.addConstr(lhs == 1, name= 'vehicle_enter_' + str(k))

約束四（customer_visit）：

for i in range(1, data.nodeNum - 1):
    lhs = LinExpr(0) 
    for k in range(data.vehicleNum):
        for j in range(1, data.nodeNum):
            if(i != j):
                lhs.addTerms(1, x[i,j,k])  
    model.addConstr(lhs == 1, name= 'customer_visit_' + str(i)) 

約束五（time_windows）：

for k in range(data.vehicleNum):
    for i in range(data.nodeNum):
        for j in range(data.nodeNum):
            if(i != j):
                model.addConstr(s[i,k] + data.disMatrix[i][j] + data.serviceTime[i] - s[j,k]- BigM + BigM * x[i,j,k] <= 0 , name= 'time_windows_')

約束六（ready_time and due_time）：

for i in range(1,data.nodeNum-1):
    for k in range(data.vehicleNum):
        model.addConstr(data.readyTime[i] <= s[i,k], name= 'ready_time')
        model.addConstr(s[i,k] <= data.dueTime[i], name= 'due_time')

約束七（capacity_vehicle）：

for k in range(data.vehicleNum):
    lhs = LinExpr(0)
    for i in range(1, data.nodeNum - 1):
        for j in range(data.nodeNum):
            if(i != j):
                lhs.addTerms(data.demand[i], x[i,j,k])
    model.addConstr(lhs <= data.capacity, name= 'capacity_vehicle' + str(k))

求解模型（solve）並輸出解：

model.optimize()
print("\n\n-----optimal value-----")
print(model.ObjVal)

for key in x.keys():
    if(x[key].x > 0 ):
        print(x[key].VarName + ' = ', x[key].x)

匯出模型：

model.write('VRPTW.lp')

求解結果（部分）：

這樣就完成了，感謝大家的閱讀。

作者：夏暘，清華大學，工業工程系/深圳國際研究生院 (碩士在讀)
劉興祿，清華大學，清華伯克利深圳學院（博士在讀）

完整Project檔案請關注運小籌公眾號，回覆【Vrptw】獲取。