容量制約付き配送計画問題（CVRP）のMTZ制約式のpythonプログラムの書き方がわかりません

Question

容量付き制約問題(CVRP)に関して、以下のURLを参考にして定式化およびpythonでのシミュレーションを試みました。

参考文献：配送計画問題をpythonで最適化する
URL:https://www.letsopt.com/entry/2020/08/30/180258
作成日：2020/8/30
作成者：cresselia2012

しかし、試みる過程で、以下のMTZ制約に関する制約式とpythonのプログラムがなぜ一致しているのか理解できませんでした。
参考にさせていただいたURL内のpythonのプログラム、特にMTZ制約が正しい理由をご教授いただけますと大変ありがたいです。

理解できなかった理由としては2点あります。

1. 本来のMTZ制約の制約式ではu[j]とdemand[j]に入力される顧客jは同様のはずですが、上記のpythonのプログラムではu[j-1]とdemand[j]と記述されており、一見異なる顧客を入力しているように思います。

### MTZ制約
for i in range(1,num_nodes):
    for j in range(1,num_nodes):
        if i != j and demand[i] + demand[j] <= capacity:
            problem += u[i-1] - u[j-1] + capacity * x[i][j] <= capacity - demand[j]

1. 変数u[i]の入力iは0を含まないと定義されています。しかし、MTZ制約に関するpythonのプログラムではi=1のとき、u[0]が生み出されてしまい、矛盾が起きてしまうのかなと疑問に思います。

## 変数を定義
u = [ pulp.LpVariable( 'u_{}'.format( i ), demand[i], capacity, cat="Integer" ) \
        for i in range(1,num_nodes) ]

### MTZ制約
for i in range(1,num_nodes):
    for j in range(1,num_nodes):
        if i != j and demand[i] + demand[j] <= capacity:
            problem += u[i-1] - u[j-1] + capacity * x[i][j] <= capacity - demand[j]

ちなみに、以下のように本来のMTZ制約の制約式のままでpythonのプログラムを書いてシミュレーションをしたところ、

for i in range(1,num_nodes):
    for j in range(1,num_nodes):
        if i != j and demand[i] + demand[j] <= capacity:
            problem += u[i] - u[j] + capacity * x[i][j] <= capacity - demand[j]

結果としては以下のようなエラーが起き、プログラムが正しくないことがわかりました。
u[i-1] - u[j-1]をu[i] - u[j]に変更したことでエラーが起こったのは明らかなのですが、具体的な原因はわかりませんでした。

Traceback (most recent call last):
  File "C:\Users\Control-PC\Documents\easy_vrp5_compver.py", line 62, in <module>
    problem += u[i] - u[j] + capacity * x[i][j] <= capacity - demand[j]
IndexError: list index out of range

最後に私がURLを参考にして記述したpythonのプログラムを記載させていただきます。
以下のプログラムのMTZ制約は参考文献と同様のプログラムを使用しており、シミュレーションを行うとエラーは起きません。

import math

def makeCVRP():
    num_nodes = 7
    capacity = 30
    demand = [ 0, 9, 11, 13, 7,\
               19, 17 ]
    coordinate = []
    coordinate.append( ( 190, 3 ) )
    coordinate.append( ( 98, 184 ) )
    coordinate.append( ( 5, 42 )  )
    coordinate.append( ( 117, 89 )  )
    coordinate.append( ( 61, 162 ) )
    coordinate.append( ( 9, 97 )  )
    coordinate.append( ( 80, 15 ) )
    def computeDistance( c1, c2 ):
        return math.sqrt( pow( c2[0] - c1[0], 2 ) + pow( c2[1] - c1[1], 2 ) )
    distance = [ [ round(computeDistance( c1, c2 )) for c1 in coordinate ] \
                    for c2 in coordinate ]
    return num_nodes, capacity, demand, distance, coordinate

# make problem
num_nodes, capacity, demand, distance, coordinate = makeCVRP()

import pulp


# 最適化問題を定義
problem = pulp.LpProblem( "CVRP", pulp.LpMinimize )

## 変数を定義
x = [ [ pulp.LpVariable( 'x_{}_{}'.format( i, j ), cat="Binary" ) \
        if i != j else None for j in range(num_nodes) ] \
        for i in range(num_nodes) ]
u = [ pulp.LpVariable( 'u_{}'.format( i ), demand[i], capacity, cat="Integer" ) \
        for i in range(1,num_nodes) ]
num_vehicle= [ pulp.LpVariable( 'num_vehicle', 1, num_nodes, cat="Integer" ) ]

## 最小化したい関数を定義
problem += pulp.lpSum( distance[i][j] * x[i][j] for i in range(num_nodes) \
                        for j in range(num_nodes) if i != j )

## 制約条件を定義
### \sum x_ij = 1
for j in range(1,num_nodes):
    problem += pulp.lpSum( x[i][j] for i in range(num_nodes) if i != j ) == 1

for i in range(1,num_nodes):
    problem += pulp.lpSum( x[i][j] for j in range(num_nodes) if i != j ) == 1

### \sum x_0j = |K|
problem += pulp.lpSum(x[0][j] for j in range(1,num_nodes)) == num_vehicle

### MTZ制約
for i in range(1,num_nodes):
    for j in range(1,num_nodes):
        if i != j:
            problem += u[i-1] - u[j-1] + capacity * x[i][j] <= capacity - demand[j]

# solve
result = problem.solve(pulp.CPLEX_CMD())

print("objective value = {}".format(pulp.value(problem.objective)))


# pulp（CBC）の結果から辺を定義
edges = [ ( i, j ) for i in range(num_nodes) for j in range(num_nodes)
            if i != j and pulp.value(x[i][j]) == 1 ]
            

# edgesから経路を取得
paths = []
for i,j in edges:
    if i == 0:
        path = [ i, j ]
        while path[-1] != 0:
            for v, u in edges:
                if v == path[-1]:
                    path.append(u)
                    break
        paths.append(path)

import networkx as nx
import matplotlib.pyplot as plt

# 有向グラフの作成
G = nx.DiGraph()
G.add_edges_from( edges )

color = [ "r", "g", "y", "m", "c" ]
edge_color = []

# 経路毎に色をセット
for i,j in G.edges:
    for t,path in enumerate(paths):
        if i in path and j in path:
            edge_color.append( color[t] )
            break
assert len(edges) == len(edge_color)

# グラフの描画
pos = { i : coordinate[i] for i in range(num_nodes) }

fig = plt.figure()
nx.draw_networkx( G, pos, edge_color=edge_color, alpha=0.5)

# 画像を保存
plt.axis("off")
fig.savefig("test.png")

print(plt.show())

初めての質問であり、質問の仕方が間違っていたら申し訳ありません。
何卒よろしくお願い申し上げます。

Answer 1

詳しいわけではありませんが、

u = [ pulp.LpVariable( 'u_{}'.format( i ), demand[i], capacity, cat="Integer" ) \
        for i in range(1,num_nodes) ]

この式で生成されるリストuは通常のpythonのリストですから0始まりですが、内部で計算に使われているiは1始まりです。
なので、uの添字0の場所にあるのは、i=1のときの値ということになります。
また、内部に含まれるdemand[1]と対応するのもu[0]となります。

Comments

[robocat]

夜分遅くに失礼いたします。

お忙しい中にもかかわらずご丁寧にご回答下さり誠にありがとうございます。
pythonの基礎的な知識が不足しており、恐れ入ります。

僭越ながら、回答を拝読させていただいた私の認識を段階的に、以下に記載させていただきます。
基本的に回答者様の回答をもとに記述させていただきましたが、私が誤った理解をしている可能性がございます。
そのため、もし認識に誤りがございましたらご指摘いただけますと幸いに存じます。
また、python初心者でして、用語の使い方が誤っておりましたらご指摘いただけますと幸いに存じます。

step1:
uとdのpython内の定義より、uの添字0の場所(uの配列内の0番目)にu[1]の値が入っており、またdの添字1の場所(dの配列内の1番目)にd[1]の値が入っている。
step2:
よって、u[1]の値とd[1]の値を対応させるためには、uの添字0の場所(uの配列内の0番目)からu[1]の値を、dの添字1の場所からd[1]の値を取り出す必要がある。
step3:
そのため、事前に定式化したMTZ制約を記述した際はu[i] - u[j] + capacity * x[i][j] <= capacity - demand[j]であったが、
今回の場合(uとdのpython内の定義的に)、事前に定式化したMTZ制約と同様の意味を持つMTZ制約をpythonで記述するには、uの添字内を-1し、u[i-1] - u[j-1] + capacity * x[i][j] <= capacity - demand[j]と記述する必要がある。

（例えばj=1をpythonのプログラムに入れてみると、uの添字0の場所から変数u[1]を、dの添字1の場所から変数d[1]を取り出すことができ、u[j-1]とdemand[j]が対応していることがわかる。）

[TakaiY]

はい。 僕の思ったとおりです。

[robocat]

お忙しい中にもかかわらずご確認下さりありがとうございます。

承知いたしました。
この度はご教授下さりありがとうございました。