python-constraint

python-constraint

python-constraintは、制約プログラミングを簡単に行うためのPythonライブラリです。

このライブラリを使うことで、制約条件を満たす組み合わせを効率的に探索することができます。

以下に、python-constraintの基本的な使い方と便利な機能について説明します。

インストール

まず、ライブラリをインストールします。

1
pip install python-constraint

基本的な使い方

以下は、基本的な制約プログラムの例です。

この例では、変数に対していくつかの制約を追加し、制約を満たす解を求めます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from constraint import *

# 問題を定義
problem = Problem()

# 変数の追加
problem.addVariable('x', range(10)) # xは0から9までの値をとる
problem.addVariable('y', range(10)) # yは0から9までの値をとる

# 制約条件の追加
def constraint_func(x, y):
return x + y == 5

problem.addConstraint(constraint_func, ['x', 'y'])

# 解を取得
solutions = problem.getSolutions()

# 解の表示
for solution in solutions:
print(solution)

便利な機能

1. 変数の範囲設定

変数の範囲をリストやジェネレータで指定できます。

1
problem.addVariable('z', [1, 2, 4, 8, 16])

2. 複数変数への同じ制約の追加

同じ制約を複数の変数に対して追加することができます。

1
2
# 全ての変数の値が異なるように制約を追加
problem.addConstraint(AllDifferentConstraint(), ['x', 'y', 'z'])

3. 事前定義された制約

python-constraintにはいくつかの事前定義された制約が用意されています。

  • AllDifferentConstraint(): すべての変数が異なる値を持つように制約する。
  • ExactSumConstraint(total): 変数の合計がtotalに等しくなるように制約する。
1
2
3
4
5
# すべての変数が異なる値を持つ
problem.addConstraint(AllDifferentConstraint(), ['x', 'y', 'z'])

# 変数の合計が10になるように制約する
problem.addConstraint(ExactSumConstraint(10), ['x', 'y', 'z'])

4. カスタム制約の追加

カスタム制約関数を定義して追加することができます。

1
2
3
4
def custom_constraint(a, b):
return a * b <= 20

problem.addConstraint(custom_constraint, ['x', 'y'])

5. 特定の値を除外

特定の値を変数から除外する制約も可能です。

1
problem.addConstraint(lambda x: x != 5, ['x'])  # xは5を取らない

実践例

以下に、より複雑な制約問題の例を示します。

この例では、簡単なナンプレ(数独)パズルの解を求めます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
from constraint import *

# 問題を定義
problem = Problem()

# 変数の追加
rows = range(1, 10)
cols = range(1, 10)

# 各セルに1から9の値を割り当てる
for row in rows:
for col in cols:
problem.addVariable((row, col), range(1, 10))

# 各行に対して全ての数が異なるように制約
for row in rows:
problem.addConstraint(AllDifferentConstraint(), [(row, col) for col in cols])

# 各列に対して全ての数が異なるように制約
for col in cols:
problem.addConstraint(AllDifferentConstraint(), [(row, col) for row in rows])

# 各3x3のボックスに対して全ての数が異なるように制約
for box_row in range(3):
for box_col in range(3):
cells = [(row, col) for row in range(box_row*3+1, box_row*3+4) for col in range(box_col*3+1, box_col*3+4)]
problem.addConstraint(AllDifferentConstraint(), cells)

# 解を取得
solutions = problem.getSolutions()

# 解の表示(最初の解のみ)
if solutions:
solution = solutions[0]
for row in rows:
print([solution[(row, col)] for col in cols])
else:
print("No solution found.")

この例では、ナンプレパズルの各セルに$1$から$9$の値を割り当て、行、列、および$3\times3$のボックス内で全ての数が異なるように制約しています。

python-constraintを使用することで、様々な制約付き問題を効率的に解くことができます。

制約を定義し、それに従って解を探索するこのライブラリの使い方を習得することで、複雑な問題にも対応できるようになります。

Google OR-Tools 4. ジョブショップスケジューリング問題(Job Shop Scheduling)

Google OR-Tools

Google OR-Toolsは、組合せ最適化問題を解決するための強力なツールキットです。

PythonC++, JavaC#などで利用でき、様々な最適化問題(線形最適化、整数最適化、車両ルーティング問題、スケジューリング問題など)に対応しています。

4. ジョブショップスケジューリング問題(Job Shop Scheduling)

ジョブショップスケジューリングは、複数のジョブとマシンがある状況で、ジョブを最適にスケジュールする問題です。

以下のソースコードは、ジョブショップスケジューリング問題をGoogle OR-ToolsCP-SATソルバーを使用して解く例です。

ジョブショップスケジューリング問題は、いくつかのジョブ(作業)があり、それぞれのジョブは特定の順序でマシンを必要とする複数のタスクで構成されています。

この問題は、すべてのジョブが完了するまでの最短時間(メイクスパン)を求めるものです。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
from ortools.sat.python import cp_model

def job_shop_scheduling():
# Create the model.
model = cp_model.CpModel()

jobs_data = [ # (machine_id, processing_time).
[(0, 3), (1, 2), (2, 2)],
[(0, 2), (2, 1), (1, 4)],
[(1, 4), (2, 3)]
]

machines_count = 3
all_tasks = {}
machine_to_intervals = [[] for _ in range(machines_count)]
horizon = sum(task[1] for job in jobs_data for task in job)

# Create job variables and constraints.
for job_id, job in enumerate(jobs_data):
for task_id, task in enumerate(job):
machine = task[0]
duration = task[1]
suffix = '_%i_%i' % (job_id, task_id)
start_var = model.NewIntVar(0, horizon, 'start' + suffix)
end_var = model.NewIntVar(0, horizon, 'end' + suffix)
interval_var = model.NewIntervalVar(start_var, duration, end_var, 'interval' + suffix)
all_tasks[job_id, task_id] = (start_var, end_var, interval_var)
machine_to_intervals[machine].append(interval_var)

# Create and add disjunctive constraints.
for machine_intervals in machine_to_intervals:
model.AddNoOverlap(machine_intervals)

# Precedences inside a job.
for job_id, job in enumerate(jobs_data):
for task_id in range(len(job) - 1):
model.Add(all_tasks[job_id, task_id + 1][0] >= all_tasks[job_id, task_id][1])

# Makespan objective.
obj_var = model.NewIntVar(0, horizon, 'makespan')
model.AddMaxEquality(obj_var, [all_tasks[job_id, len(job) - 1][1] for job_id, job in enumerate(jobs_data)])
model.Minimize(obj_var)

# Solve the model.
solver = cp_model.CpSolver()
status = solver.Solve(model)

# Print solution.
if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
print('Optimal Schedule Length: %i' % solver.ObjectiveValue())
print()
for job_id, job in enumerate(jobs_data):
print('Job %i:' % job_id)
for task_id, task in enumerate(job):
start_var, end_var, _ = all_tasks[job_id, task_id]
print('Task %i starts at %i and ends at %i' % (
task_id, solver.Value(start_var), solver.Value(end_var)))
else:
print('No solution found.')

job_shop_scheduling()

ソースコード解説

ソースコードの各部分を説明します。

1. インポートと関数定義

1
2
3
4
5
from ortools.sat.python import cp_model

def job_shop_scheduling():
# Create the model.
model = cp_model.CpModel()

まず、Google OR-ToolsCP-SATソルバーからcp_modelモジュールをインポートし、スケジューリングを行う関数job_shop_schedulingを定義します。

その後、制約プログラミングモデルを作成します。

2. ジョブデータの定義

1
2
3
4
5
6
7
8
9
10
jobs_data = [  # (machine_id, processing_time).
[(0, 3), (1, 2), (2, 2)],
[(0, 2), (2, 1), (1, 4)],
[(1, 4), (2, 3)]
]

machines_count = 3
all_tasks = {}
machine_to_intervals = [[] for _ in range(machines_count)]
horizon = sum(task[1] for job in jobs_data for task in job)

ジョブデータを定義します。

ここでは、各ジョブはタスクのリストで構成されており、各タスクは使用するマシンIDと処理時間を持っています。

また、マシンの数、すべてのタスクを保持する辞書、各マシンに対するタスクのリスト、およびスケジュールの全期間(ホライズン)を計算します。

3. 変数と制約の作成

1
2
3
4
5
6
7
8
9
10
11
# Create job variables and constraints.
for job_id, job in enumerate(jobs_data):
for task_id, task in enumerate(job):
machine = task[0]
duration = task[1]
suffix = '_%i_%i' % (job_id, task_id)
start_var = model.NewIntVar(0, horizon, 'start' + suffix)
end_var = model.NewIntVar(0, horizon, 'end' + suffix)
interval_var = model.NewIntervalVar(start_var, duration, end_var, 'interval' + suffix)
all_tasks[job_id, task_id] = (start_var, end_var, interval_var)
machine_to_intervals[machine].append(interval_var)

各タスクに対して、開始時間、終了時間、およびそれらを結びつけるインターバル変数を作成します。

これらの変数をall_tasks辞書に格納し、対応するマシンのインターバルリストに追加します。

4. 排他制約の追加

1
2
3
# Create and add disjunctive constraints.
for machine_intervals in machine_to_intervals:
model.AddNoOverlap(machine_intervals)

同じマシンで同時に複数のタスクが実行されないようにするため、AddNoOverlap制約を追加します。

この制約は、同じマシンに割り当てられたすべてのインターバル変数に対して適用されます。

5. ジョブ内の順序制約の追加

1
2
3
4
# Precedences inside a job.
for job_id, job in enumerate(jobs_data):
for task_id in range(len(job) - 1):
model.Add(all_tasks[job_id, task_id + 1][0] >= all_tasks[job_id, task_id][1])

各ジョブ内のタスクが定められた順序で実行されるようにするための制約を追加します。

具体的には、あるタスクの終了時間が次のタスクの開始時間よりも早くなるようにします。

6. メイクスパンの目標設定

1
2
3
4
# Makespan objective.
obj_var = model.NewIntVar(0, horizon, 'makespan')
model.AddMaxEquality(obj_var, [all_tasks[job_id, len(job) - 1][1] for job_id, job in enumerate(jobs_data)])
model.Minimize(obj_var)

メイクスパン(全ジョブが完了するまでの最短時間)を最小化することを目標として設定します。

AddMaxEquality制約を使用して、メイクスパン変数がすべてのジョブの最終タスクの終了時間の最大値に等しくなるようにします。

7. モデルの解決

1
2
3
# Solve the model.
solver = cp_model.CpSolver()
status = solver.Solve(model)

CP-SATソルバーを作成し、モデルを解決します。

8. 解の表示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
    # Print solution.
if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
print('Optimal Schedule Length: %i' % solver.ObjectiveValue())
print()
for job_id, job in enumerate(jobs_data):
print('Job %i:' % job_id)
for task_id, task in enumerate(job):
start_var, end_var, _ = all_tasks[job_id, task_id]
print('Task %i starts at %i and ends at %i' % (
task_id, solver.Value(start_var), solver.Value(end_var)))
else:
print('No solution found.')

job_shop_scheduling()

解が見つかった場合、そのメイクスパンと各タスクの開始時間と終了時間を表示します。

解が見つからなかった場合は、その旨を表示します。


このコード全体で、ジョブショップスケジューリング問題を解決するための制約プログラミングモデルを定義し、最適なスケジュールを求めています。

結果解説

[実行結果]

Optimal Schedule Length: 11

Job 0:
Task 0 starts at 2 and ends at 5
Task 1 starts at 5 and ends at 7
Task 2 starts at 7 and ends at 9
Job 1:
Task 0 starts at 0 and ends at 2
Task 1 starts at 2 and ends at 3
Task 2 starts at 7 and ends at 11
Job 2:
Task 0 starts at 0 and ends at 4
Task 1 starts at 4 and ends at 7

この結果は、ジョブショップスケジューリング問題に対して最適なスケジュールが見つかり、そのメイクスパン(全ジョブが完了するまでの最短時間)が$11$であることを示しています。

以下、解説します。

1. メイクスパンの解説

Optimal Schedule Length: 11

  • メイクスパンは$11$です。
    これは、すべてのジョブが完了するまでに最短で$11$単位時間かかることを意味します。

2. ジョブ0のスケジュール

Job 0:

  • Task 0 starts at 2 and ends at 5: このタスクはマシン$0$で処理され、処理時間は$3$単位時間です。
  • Task 1 starts at 5 and ends at 7: このタスクはマシン$1$で処理され、処理時間は$2$単位時間です。
  • Task 2 starts at 7 and ends at 9: このタスクはマシン$2$で処理され、処理時間は$2$単位時間です。

3. ジョブ1のスケジュール

Job 1:

  • Task 0 starts at 0 and ends at 2: このタスクはマシン$0$で処理され、処理時間は$2$単位時間です。
  • Task 1 starts at 2 and ends at 3: このタスクはマシン$2$で処理され、処理時間は$1$単位時間です。
  • Task 2 starts at 7 and ends at 11: このタスクはマシン$1$で処理され、処理時間は$4$単位時間です。

4. ジョブ2のスケジュール

Job 2:

  • Task 0 starts at 0 and ends at 4: このタスクはマシン$1$で処理され、処理時間は$4$単位時間です。
  • Task 1 starts at 4 and ends at 7: このタスクはマシン$2$で処理され、処理時間は$3$単位時間です。

5. スケジュール全体の概要

このスケジュールにより、以下のようなタスクの配置が見られます:

  • マシン0:

    • Job 1のTask 0: $0$から$2$
    • Job 0のTask 0: $2$から$5$
  • マシン1:

    • Job 2のTask 0: $0$から$4$
    • Job 0のTask 1: $5$から$7$
    • Job 1のTask 2: $7$から$11$
  • マシン2:

    • Job 1のTask 1: $2$から$3$
    • Job 2のTask 1: $4$から$7$
    • Job 0のTask 2: $7$から$9$

このように、各マシンでタスクが重ならないようにスケジュールされています。

また、ジョブ内のタスクが指定された順序で実行されています。

結果として、このスケジュールによりすべてのジョブが最短時間($11$単位時間)で完了することが確認できます。

Google OR-Tools 3. 車両ルーティング問題(Vehicle Routing Problem, VRP)

Google OR-Tools

Google OR-Toolsは、組合せ最適化問題を解決するための強力なツールキットです。

PythonC++, JavaC#などで利用でき、様々な最適化問題(線形最適化、整数最適化、車両ルーティング問題、スケジューリング問題など)に対応しています。

3. 車両ルーティング問題(Vehicle Routing Problem, VRP)

VRPは、指定された地点を巡回する車両の最適なルートを見つける問題です。

サンプルコード:VRP

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
from ortools.constraint_solver import pywrapcp
from ortools.constraint_solver import routing_enums_pb2

def create_data_model():
data = {}
data['distance_matrix'] = [
[0, 2, 9, 10],
[1, 0, 6, 4],
[15, 7, 0, 8],
[6, 3, 12, 0],
]
data['num_vehicles'] = 1
data['depot'] = 0
return data

def print_solution(manager, routing, solution):
print('Objective: {}'.format(solution.ObjectiveValue()))
index = routing.Start(0)
plan_output = 'Route for vehicle 0:\n'
route_distance = 0
while not routing.IsEnd(index):
plan_output += ' {} ->'.format(manager.IndexToNode(index))
previous_index = index
index = solution.Value(routing.NextVar(index))
route_distance += routing.GetArcCostForVehicle(previous_index, index, 0)
plan_output += ' {}\n'.format(manager.IndexToNode(index))
plan_output += 'Distance of the route: {}m\n'.format(route_distance)
print(plan_output)

def vehicle_routing_problem():
data = create_data_model()
manager = pywrapcp.RoutingIndexManager(len(data['distance_matrix']), data['num_vehicles'], data['depot'])
routing = pywrapcp.RoutingModel(manager)
def distance_callback(from_index, to_index):
from_node = manager.IndexToNode(from_index)
to_node = manager.IndexToNode(to_index)
return data['distance_matrix'][from_node][to_node]
transit_callback_index = routing.RegisterTransitCallback(distance_callback)
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)
search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)
solution = routing.SolveWithParameters(search_parameters)
if solution:
print_solution(manager, routing, solution)

vehicle_routing_problem()

ソースコード解説

ソースコードの各部分を説明します。

1. インポート

1
2
from ortools.constraint_solver import pywrapcp
from ortools.constraint_solver import routing_enums_pb2
  • pywrapcpOR-Tools制約ソルバーを提供します。
  • routing_enums_pb2 はルーティング用の定数を提供します。

2. データモデルの作成

1
2
3
4
5
6
7
8
9
10
11
def create_data_model():
data = {}
data['distance_matrix'] = [
[0, 2, 9, 10],
[1, 0, 6, 4],
[15, 7, 0, 8],
[6, 3, 12, 0],
]
data['num_vehicles'] = 1
data['depot'] = 0
return data
  • create_data_model 関数は、距離行列や車両の数、出発点(デポ)などのデータを定義します。
  • distance_matrix は各地点間の距離を示します。
  • num_vehicles は使用する車両の数を示します。
  • depot は出発点のインデックスを示します。

3. 解の表示

1
2
3
4
5
6
7
8
9
10
11
12
13
def print_solution(manager, routing, solution):
print('Objective: {}'.format(solution.ObjectiveValue()))
index = routing.Start(0)
plan_output = 'Route for vehicle 0:\n'
route_distance = 0
while not routing.IsEnd(index):
plan_output += ' {} ->'.format(manager.IndexToNode(index))
previous_index = index
index = solution.Value(routing.NextVar(index))
route_distance += routing.GetArcCostForVehicle(previous_index, index, 0)
plan_output += ' {}\n'.format(manager.IndexToNode(index))
plan_output += 'Distance of the route: {}m\n'.format(route_distance)
print(plan_output)
  • print_solution 関数は、解(ルートと距離)を出力します。
  • solution.ObjectiveValue()目的関数の値を取得します。
  • routing.Start(0) で車両$ 0 $のルートを開始します。
  • ループ内で、各地点のインデックスを取得し、次の地点へ進みます。
  • ルート全体の距離を計算し、最終的にルート距離を出力します。

4. 車両経路問題の定義と解の計算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def vehicle_routing_problem():
data = create_data_model()
manager = pywrapcp.RoutingIndexManager(len(data['distance_matrix']), data['num_vehicles'], data['depot'])
routing = pywrapcp.RoutingModel(manager)

def distance_callback(from_index, to_index):
from_node = manager.IndexToNode(from_index)
to_node = manager.IndexToNode(to_index)
return data['distance_matrix'][from_node][to_node]

transit_callback_index = routing.RegisterTransitCallback(distance_callback)
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)

solution = routing.SolveWithParameters(search_parameters)
if solution:
print_solution(manager, routing, solution)

vehicle_routing_problem()
  • vehicle_routing_problem 関数は、車両経路問題を定義し、解を計算します。
  • create_data_model を呼び出してデータを取得します。
  • RoutingIndexManager を使用してルーティングインデックスマネージャを作成します。
  • RoutingModel を使用してルーティングモデルを作成します。
  • distance_callback 関数を定義し、各地点間の距離を返します。
  • RegisterTransitCallback を使用してコールバックを登録し、距離行列を評価します。
  • SetArcCostEvaluatorOfAllVehicles を使用してすべての車両のアークコスト評価器を設定します。
  • DefaultRoutingSearchParameters を使用して検索パラメータを設定し、最初の解法戦略を指定します。
  • SolveWithParameters を使用して問題を解き、解が見つかった場合は print_solution を呼び出して結果を表示します。

5. 全体の流れ

  1. データモデルを作成し、距離行列車両情報を定義します。
  2. ルーティングインデックスマネージャとルーティングモデルを作成します。
  3. 距離評価用のコールバック関数を定義し、登録します。
  4. 検索パラメータを設定し、最適なルートを求めます。
  5. 解が見つかった場合、ルート距離を出力します。

以上がこのコードの詳細な説明です。

結果解説

[実行結果]

Objective: 21
Route for vehicle 0:
 0 -> 2 -> 3 -> 1 -> 0
Distance of the route: 21m

詳細な解説

  1. Objective: 21

    • 目的関数の値を示しています。
      この場合、車両が訪問するルートの総距離が$21$メートルであることを意味します。
      ここでは最小化された距離が$21$メートルという結果を得ています。
  2. Route for vehicle 0: 0 -> 2 -> 3 -> 1 -> 0

    • ルートの詳細を示しています。
      車両$0$が訪問する順番を示しており、出発点から始まり、各ノードを訪問し、最後に出発点に戻ります。
      • 0 は出発点(デポ)。
      • 2 は2番目の地点。
      • 3 は3番目の地点。
      • 1 は1番目の地点。
      • 0 は再び出発点に戻ります。
  3. Distance of the route: 21m

    • ルート全体の距離が$21$メートルであることを示しています。
      ルートの距離は、各ステップ間の距離の合計です。

Google OR-Tools 2. 整数最適化(Integer Optimization)

Google OR-Tools

Google OR-Toolsは、組合せ最適化問題を解決するための強力なツールキットです。

PythonC++, JavaC#などで利用でき、様々な最適化問題(線形最適化、整数最適化、車両ルーティング問題、スケジューリング問題など)に対応しています。

2. 整数最適化(Integer Optimization)

整数最適化では、変数が整数値を取るような最適化問題を解きます。

サンプルコード:整数最適化

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
from ortools.linear_solver import pywraplp

def integer_optimization_example():
# Solverの作成
solver = pywraplp.Solver.CreateSolver('SCIP')
if not solver:
return

# 変数の作成
x = solver.IntVar(0, 10, 'x')
y = solver.IntVar(0, 10, 'y')

# 制約条件の追加
solver.Add(2 * x + 3 * y <= 12)

# 目的関数の設定
solver.Maximize(3 * x + 4 * y)

# 解を求める
status = solver.Solve()

# 結果を表示
if status == pywraplp.Solver.OPTIMAL:
print('Objective value =', solver.Objective().Value())
print('x =', x.solution_value())
print('y =', y.solution_value())
else:
print('The problem does not have an optimal solution.')

integer_optimization_example()

ソースコード解説

1. pywraplp モジュールのインポート

最初に、ortools.linear_solver モジュールから pywraplp をインポートします。
このモジュールは、Google OR-Toolsの線形および整数ソルバーを使用するためのインターフェースを提供します。

1
from ortools.linear_solver import pywraplp

2. ソルバーの作成

次に、整数最適化問題を解くためのソルバーを作成します。
CreateSolver('SCIP') メソッドは、SCIP(Solving Constraint Integer Programs)ソルバーを作成します。

1
2
3
solver = pywraplp.Solver.CreateSolver('SCIP')
if not solver:
return
  • solver が作成されない場合、None が返され、関数を終了します。

3. 変数の作成

ソルバーで使用する整数変数を作成します。
IntVar メソッドを使って変数を定義します。
ここでは、変数 xy がそれぞれ0から10の範囲に制限されています。

1
2
x = solver.IntVar(0, 10, 'x')
y = solver.IntVar(0, 10, 'y')
  • xy は整数変数で、それぞれ$0$から$10$の範囲の整数値を取ります。

4. 制約条件の追加

問題の制約条件を追加します。

ここでは、2 * x + 3 * y <= 12 という制約を追加しています。

1
solver.Add(2 * x + 3 * y <= 12)
  • この制約は、2倍の x と3倍の y の合計が$12$以下であることを要求します。

5. 目的関数の設定

最適化する目的関数を設定します。ここでは、3 * x + 4 * y を最大化するように設定しています。

1
solver.Maximize(3 * x + 4 * y)
  • 目的関数は 3 * x + 4 * y で、これを最大化します。

6. 解を求める

ソルバーを呼び出して、最適解を求めます。

Solve メソッドは、ソルバーを実行して最適化問題を解きます。

1
status = solver.Solve()
  • Solve() メソッドは、ソルバーを実行し、解が見つかったかどうかを示すステータスコードを返します。

7. 結果の表示

最適解が見つかったかどうかを status で確認し、結果を表示します。

pywraplp.Solver.OPTIMAL は、問題に最適解が見つかったことを示します。

1
2
3
4
5
6
if status == pywraplp.Solver.OPTIMAL:
print('Objective value =', solver.Objective().Value())
print('x =', x.solution_value())
print('y =', y.solution_value())
else:
print('The problem does not have an optimal solution.')
  • statuspywraplp.Solver.OPTIMAL であれば、最適解が見つかっています。
  • solver.Objective().Value()目的関数の値を取得します。
  • x.solution_value() および y.solution_value() で各変数の最適解を取得します。

8. 関数の実行

最後に、定義した関数を実行します。

1
integer_optimization_example()

このソースコード全体は、SCIPソルバーを使用して、整数最適化問題を解き、その結果を表示するプロセスを示しています。

結果解説

[実行結果]

Objective value = 18.0
x = 6.0
y = 0.0

この結果は、与えられた整数最適化問題の最適解を示しています。

結果の詳細

1
2
3
Objective value = 18.0
x = 6.0
y = 0.0

問題の設定の再確認

まず、ソースコードで設定した問題の内容を再確認します。

  • 変数の範囲:

    • x は$0$から$10$の範囲の整数
    • y は$0$から$10$の範囲の整数
  • 制約条件:

    • 2 * x + 3 * y <= 12
  • 目的関数:

    • 3 * x + 4 * y を最大化

結果の意味

  1. Objective value = 18.0

    • 目的関数 3 * x + 4 * y の最大値が$18$であることを示しています。
  2. x = 6.0

    • 変数 x の最適値が$6$であることを示しています。
  3. y = 0.0

    • 変数 y の最適値が$0$であることを示しています。

制約条件の確認

最適解 (x = 6, y = 0)制約条件を満たしていることを確認します。

  • 制約条件: 2 * x + 3 * y <= 12
    • 2 * 6 + 3 * 0 = 12 なので、この制約を満たしています。

目的関数の値の計算

最適解 (x = 6, y = 0) に対する目的関数の値を計算します。

  • 目的関数: 3 * x + 4 * y
    • 3 * 6 + 4 * 0 = 18

これにより、目的関数の値が$18$であることが確認できます。

結果の解釈

この結果から、以下のことがわかります:

  • 制約条件 2 * x + 3 * y <= 12 を満たしながら、目的関数 3 * x + 4 * y最大化するためには、x = 6y = 0 に設定するのが最適であり、そのときの目的関数の値は$18$である。
  • つまり、x を6に、y を$0$にすることで、与えられた制約条件の下で最大の利益(目的関数の値)を得ることができます。

このようにして、整数最適化問題の最適解を求め、その結果を評価することができます。

Google OR-Tools 1. 線形最適化(Linear Optimization)

Google OR-Tools

Google OR-Toolsは、組合せ最適化問題を解決するための強力なツールキットです。

PythonC++, JavaC#などで利用でき、様々な最適化問題(線形最適化、整数最適化、車両ルーティング問題、スケジューリング問題など)に対応しています。

1. 線形最適化(Linear Optimization)

線形最適化では、線形の目的関数最大化または最小化する問題を解きます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
!pip install ortools
from ortools.linear_solver import pywraplp

def linear_optimization_example():
# Solverの作成
solver = pywraplp.Solver.CreateSolver('GLOP')
if not solver:
return

# 変数の作成
x = solver.NumVar(0, 1, 'x')
y = solver.NumVar(0, 2, 'y')

# 制約条件の追加
solver.Add(x + y <= 2)

# 目的関数の設定
solver.Maximize(x + y)

# 解を求める
status = solver.Solve()

# 結果を表示
if status == pywraplp.Solver.OPTIMAL:
print('Objective value =', solver.Objective().Value())
print('x =', x.solution_value())
print('y =', y.solution_value())
else:
print('The problem does not have an optimal solution.')

linear_optimization_example()

コード解説

ソースコードを説明します。

1. インポート

1
from ortools.linear_solver import pywraplp

Google OR-Toolsの線形最適化ライブラリであるpywraplpをインポートします。

これにより、線形最適化問題を解くための機能を利用できます。

2. 関数定義

1
def linear_optimization_example():

線形最適化問題の例を示す関数linear_optimization_exampleを定義します。

3. ソルバーの作成

1
2
3
solver = pywraplp.Solver.CreateSolver('GLOP')
if not solver:
return

ここでは、GLOP(Google Linear Optimization Package)という線形最適化ソルバーを作成します。

CreateSolver('GLOP')は新しいソルバーオブジェクトを作成します。

ソルバーが作成できなかった場合、関数を終了します。

4. 変数の作成

1
2
x = solver.NumVar(0, 1, 'x')
y = solver.NumVar(0, 2, 'y')

線形最適化問題の変数を作成します。ここでは、変数xyを定義します。

  • xは$0$から$1$までの値を取ります。
  • yは$0$から$2$までの値を取ります。

5. 制約条件の追加

1
solver.Add(x + y <= 2)

変数に対する制約を追加します。ここでは、x + y <= 2という制約条件を設定しています。

これは、変数xyの和が2以下でなければならないことを意味します。

6. 目的関数の設定

1
solver.Maximize(x + y)

最適化の目的を設定します。ここでは、x + yを最大化するように設定しています。

7. 解を求める

1
status = solver.Solve()

設定した制約条件目的関数に基づいて、最適解を求めます。

このメソッドは、解が見つかったかどうかを示すステータスコードを返します。

8. 結果の表示

1
2
3
4
5
6
if status == pywraplp.Solver.OPTIMAL:
print('Objective value =', solver.Objective().Value())
print('x =', x.solution_value())
print('y =', y.solution_value())
else:
print('The problem does not have an optimal solution.')

最適解が見つかった場合、目的関数の値と変数xおよびyの値を表示します。

解が見つからなかった場合は、最適解が存在しないことを表示します。

9. 関数の呼び出し

1
linear_optimization_example()

定義した関数を呼び出し、線形最適化問題を実行します。

まとめ

このコードは、OR-Toolsを使用して簡単な線形最適化問題を解決する方法を示しています。

具体的には、2つの変数xyを定義し、制約条件x + y <= 2のもとでx + yを最大化する問題を解いています。

最適解が見つかった場合、その解を表示します。

結果解説

[実行結果]

Objective value = 2.0
x = 0.0
y = 2.0

この結果は、線形最適化問題の解として得られた最適解を示しています。
各部分の意味を以下に説明します。

結果の詳細

Objective value = 2.0

目的関数の値が$2.0$であることを示しています。
つまり、x + yを最大化する目的関数に対して得られた最大値が2.0です。

x = 0.0

変数xの最適解が$0.0$であることを示しています。

y = 2.0

変数yの最適解が$2.0$であることを示しています。

解の詳細

この結果に基づいて、以下の点を解説します。

  1. 目的関数の最大化:

    • 目的関数はx + yを最大化することを目指しています。
  2. 制約条件の確認:

    • 制約条件はx + y <= 2です。
      この条件により、xyの和が$2$以下でなければなりません。
  3. 最適解の妥当性:

    • 得られた解はx = 0.0y = 2.0です。
    • この解は制約条件x + y <= 2を満たしており、0.0 + 2.0 = 2.0となっています。
      制約条件の上限である$2$に一致しています。
  4. 最大化の達成:

    • 目的関数x + yを最大化した結果が$2.0$であることは、制約条件の下で可能な最大値であることを示しています。
    • 他の値の組み合わせ(例えば、x = 1y = 1など)では、x + yの値が2を超えることはありません。

結論

この結果は、与えられた線形最適化問題の最適解として、変数xが$0.0$、変数yが$2.0$であり、そのときの目的関数x + yの値が最大で$2.0$であることを示しています。

制約条件x + y <= 2を満たしながら、目的関数の値を最大化するための最適な値の組み合わせがx = 0.0y = 2.0であることが確認されました。

DEAP(Distributed Evolutionary Algorithms)

DEAP(Distributed Evolutionary Algorithms)

DEAP(Distributed Evolutionary Algorithms)は、遺伝的アルゴリズム(GA)進化的計算を実装するための強力で柔軟なライブラリです。

以下に、DEAPの便利な使い方をいくつか紹介します。

インストール

まず、DEAPをインストールします。

1
pip install deap

基本的な使い方

1. 必要なモジュールのインポート

1
2
import random
from deap import base, creator, tools, algorithms

2. 問題の設定

最大化問題の適応度を定義します。

1
2
creator.create("FitnessMax", base.Fitness, weights=(1.0,))
creator.create("Individual", list, fitness=creator.FitnessMax)

3. 個体と集団の初期化

1
2
3
4
toolbox = base.Toolbox()
toolbox.register("attr_float", random.uniform, -10, 10)
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, n=10)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)

4. 適応度関数の定義

1
2
3
4
def eval_function(individual):
return sum(individual),

toolbox.register("evaluate", eval_function)

5. 遺伝的操作の定義

1
2
3
toolbox.register("mate", tools.cxTwoPoint)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.2)
toolbox.register("select", tools.selTournament, tournsize=3)

6. 進化の実行

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
def main():
random.seed(42)
population = toolbox.population(n=300)
ngen = 40
cxpb = 0.5
mutpb = 0.2

# 適応度の計算
fitnesses = list(map(toolbox.evaluate, population))
for ind, fit in zip(population, fitnesses):
ind.fitness.values = fit

# 進化のループ
for g in range(ngen):
print(f"-- Generation {g} --")

# 選択
offspring = toolbox.select(population, len(population))
offspring = list(map(toolbox.clone, offspring))

# 交叉と突然変異
for child1, child2 in zip(offspring[::2], offspring[1::2]):
if random.random() < cxpb:
toolbox.mate(child1, child2)
del child1.fitness.values
del child2.fitness.values

for mutant in offspring:
if random.random() < mutpb:
toolbox.mutate(mutant)
del mutant.fitness.values

# 適応度の再計算
invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
fitnesses = map(toolbox.evaluate, invalid_ind)
for ind, fit in zip(invalid_ind, fitnesses):
ind.fitness.values = fit

population[:] = offspring

# 統計の表示
fits = [ind.fitness.values[0] for ind in population]
length = len(population)
mean = sum(fits) / length
sum2 = sum(x*x for x in fits)
std = abs(sum2 / length - mean**2)**0.5

print(f"Min: {min(fits)}, Max: {max(fits)}, Avg: {mean}, Std: {std}")

best_ind = tools.selBest(population, 1)[0]
print(f"Best individual is {best_ind}, {best_ind.fitness.values}")

if __name__ == "__main__":
main()

ソースコードの詳細を説明します。

初期設定

1
2
3
4
5
random.seed(42)
population = toolbox.population(n=300)
ngen = 40
cxpb = 0.5
mutpb = 0.2
  • random.seed(42):ランダムな値の生成に一貫性を持たせるためにシード値を設定します。
  • population = toolbox.population(n=300):$300$個体の初期集団を生成します。
    toolbox.populationは、遺伝的アルゴリズムのために設定された関数です。
  • ngen = 40:進化の世代数を$40$に設定します。
  • cxpb = 0.5:交叉(crossover)の確率を$0.5$に設定します。
  • mutpb = 0.2:突然変異(mutation)の確率を$0.2$に設定します。

適応度の計算

1
2
3
fitnesses = list(map(toolbox.evaluate, population))
for ind, fit in zip(population, fitnesses):
ind.fitness.values = fit
  • fitnesses = list(map(toolbox.evaluate, population)):初期集団の各個体の適応度を評価し、リストに格納します。
  • for ind, fit in zip(population, fitnesses):各個体の適応度を設定します。

進化のループ

1
2
3
4
5
6
for g in range(ngen):
print(f"-- Generation {g} --")

# 選択
offspring = toolbox.select(population, len(population))
offspring = list(map(toolbox.clone, offspring))
  • for g in range(ngen):40世代の進化を繰り返します。
  • offspring = toolbox.select(population, len(population)):次世代の個体を選択します。
  • offspring = list(map(toolbox.clone, offspring)):選択された個体をクローンします。

交叉と突然変異

1
2
3
4
5
6
7
8
9
10
for child1, child2 in zip(offspring[::2], offspring[1::2]):
if random.random() < cxpb:
toolbox.mate(child1, child2)
del child1.fitness.values
del child2.fitness.values

for mutant in offspring:
if random.random() < mutpb:
toolbox.mutate(mutant)
del mutant.fitness.values
  • for child1, child2 in zip(offspring[::2], offspring[1::2]):2つずつペアにして交叉を行います。
  • if random.random() < cxpb交叉確率に基づいて交叉を実行します。
    交叉した個体の適応度を削除します(再評価のため)。
  • for mutant in offspring突然変異を行います。
  • if random.random() < mutpb突然変異確率に基づいて突然変異を実行します。
    突然変異した個体の適応度を削除します(再評価のため)。

適応度の再計算

1
2
3
4
5
6
invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
fitnesses = map(toolbox.evaluate, invalid_ind)
for ind, fit in zip(invalid_ind, fitnesses):
ind.fitness.values = fit

population[:] = offspring
  • invalid_ind = [ind for ind in offspring if not ind.fitness.valid]:適応度が無効な個体をリストにします。
  • fitnesses = map(toolbox.evaluate, invalid_ind):無効な個体の適応度を再評価します。
  • for ind, fit in zip(invalid_ind, fitnesses):再評価された適応度を個体に設定します。
  • population[:] = offspring:新しい世代の個体で集団を置き換えます。

統計の表示

1
2
3
4
5
6
7
fits = [ind.fitness.values[0] for ind in population]
length = len(population)
mean = sum(fits) / length
sum2 = sum(x*x for x in fits)
std = abs(sum2 / length - mean**2)**0.5

print(f"Min: {min(fits)}, Max: {max(fits)}, Avg: {mean}, Std: {std}")
  • fits = [ind.fitness.values[0] for ind in population]:集団の各個体の適応度を取得します。
  • mean:適応度の平均値を計算します。
  • sum2:適応度の二乗和を計算します。
  • std:標準偏差を計算します。
  • 統計情報(最小値、最大値、平均値、標準偏差)を表示します。

最良個体の表示

1
2
best_ind = tools.selBest(population, 1)[0]
print(f"Best individual is {best_ind}, {best_ind.fitness.values}")
  • best_ind = tools.selBest(population, 1)[0]:最も適応度の高い個体を選択します。
  • 最良個体とその適応度を表示します。

このコード全体は、遺伝的アルゴリズムを用いて集団の適応度を向上させるプロセスを実行しています。

交叉突然変異によって新しい個体を生成し、選択と適応度の再計算を繰り返して、最適解を探索します。

結果解説

[実行結果]

-- Generation 0 --
Min: -23.62122528012338, Max: 63.23216847796773, Avg: 13.94750872866429, Std: 14.117844855355854
-- Generation 1 --
Min: -15.831634197388622, Max: 66.00568167314223, Avg: 25.843928872528465, Std: 12.001841010730187
-- Generation 2 --
Min: 1.571510798629661, Max: 67.7416705994344, Avg: 36.10924261516296, Std: 10.639830457934638
-- Generation 3 --
Min: 17.251545608296084, Max: 70.7377667641021, Avg: 45.18612876125583, Std: 10.224876102168013
-- Generation 4 --
Min: 29.187935928369605, Max: 76.1192756677395, Avg: 53.98213674983369, Std: 8.703585040850301
-- Generation 5 --
Min: 41.67182338605669, Max: 77.93473415749817, Avg: 61.70199699017955, Std: 6.187469589825239
-- Generation 6 --
Min: 52.09535714999056, Max: 78.79192555588493, Avg: 66.62043085395645, Std: 4.748183156856532
-- Generation 7 --
Min: 58.486933075200305, Max: 83.0666536496249, Avg: 70.36244867844377, Std: 4.261643906702605
-- Generation 8 --
Min: 63.62749493842264, Max: 85.16722313166508, Avg: 74.01199704669335, Std: 3.915234083237426
-- Generation 9 --
Min: 67.66567852603893, Max: 87.4679125976825, Avg: 77.51393726261347, Std: 3.6385551917724834
-- Generation 10 --
Min: 70.08903900813345, Max: 90.61910944047975, Avg: 80.6030591219591, Std: 3.427952443409654
-- Generation 11 --
Min: 73.06776847431156, Max: 95.86148332348742, Avg: 83.57800102942038, Std: 2.824878578319968
-- Generation 12 --
Min: 77.12960155059386, Max: 95.86148332348742, Avg: 85.74637442049743, Std: 2.5049934648819057
-- Generation 13 --
Min: 82.22740845194737, Max: 95.86148332348742, Avg: 87.80128821464989, Std: 2.320462412242403
-- Generation 14 --
Min: 83.51896359405892, Max: 97.09479936204475, Avg: 89.75269302806925, Std: 2.438380774085663
-- Generation 15 --
Min: 86.73351504224718, Max: 98.4456927207268, Avg: 91.818120209583, Std: 2.545988726836644
-- Generation 16 --
Min: 88.42499238474882, Max: 100.61389792354642, Avg: 94.09226319826075, Std: 2.285888831052674
-- Generation 17 --
Min: 89.24757995715294, Max: 102.6838520503081, Avg: 96.01564578409, Std: 1.9067732850766526
-- Generation 18 --
Min: 91.86451940577689, Max: 105.13304887155599, Avg: 97.59405316472665, Std: 1.8201684526510389
-- Generation 19 --
Min: 93.09487363840095, Max: 105.37632781629395, Avg: 98.87007574977343, Std: 1.8358356370881441
-- Generation 20 --
Min: 95.24057602651627, Max: 107.56442476215851, Avg: 100.38936719083868, Std: 1.9591513525426185
-- Generation 21 --
Min: 95.29981673148842, Max: 108.63304220299476, Avg: 102.00671712709355, Std: 2.0950703250083142
-- Generation 22 --
Min: 96.9630117534801, Max: 110.31319320264738, Avg: 103.90216297981273, Std: 2.067124226863785
-- Generation 23 --
Min: 100.13912350972645, Max: 113.4146271725918, Avg: 105.59576786559764, Std: 2.047165583227913
-- Generation 24 --
Min: 102.0215365337707, Max: 115.17904279346358, Avg: 107.4314699463842, Std: 2.3470974433739396
-- Generation 25 --
Min: 101.42703170477672, Max: 116.50175575861405, Avg: 109.58206484666111, Std: 2.4236312803244817
-- Generation 26 --
Min: 101.55536223691173, Max: 118.77525231829321, Avg: 111.70036923801922, Std: 2.28186217906567
-- Generation 27 --
Min: 107.3971040840559, Max: 121.77595229614845, Avg: 113.44641244359892, Std: 2.204699713517585
-- Generation 28 --
Min: 108.38391918591364, Max: 121.77595229614845, Avg: 115.29798093245587, Std: 2.071478865376343
-- Generation 29 --
Min: 110.67133487335356, Max: 123.20042811966078, Avg: 116.85867616964943, Std: 2.1699761375316644
-- Generation 30 --
Min: 111.89626831109847, Max: 125.59998028693536, Avg: 118.75443109463095, Std: 2.2065782226142097
-- Generation 31 --
Min: 112.3320132990261, Max: 131.74710678901693, Avg: 120.59690320624509, Std: 2.374984092571527
-- Generation 32 --
Min: 116.66719981073696, Max: 131.74710678901693, Avg: 122.65200536117692, Std: 2.4849500246592715
-- Generation 33 --
Min: 117.25502174406162, Max: 132.02720181664301, Avg: 124.73500348653724, Std: 2.4570849635967265
-- Generation 34 --
Min: 119.70615322562539, Max: 135.4831278203402, Avg: 126.87729177862992, Std: 2.2600095543595398
-- Generation 35 --
Min: 121.96222528366903, Max: 137.06036837117406, Avg: 128.61471631011202, Std: 2.371858244077739
-- Generation 36 --
Min: 123.52343627859534, Max: 140.32944805635685, Avg: 130.8766887732697, Std: 2.5972125655787104
-- Generation 37 --
Min: 124.32891522567849, Max: 140.32944805635685, Avg: 133.02356239681754, Std: 2.7282315847206964
-- Generation 38 --
Min: 128.16134079141833, Max: 142.86296023338744, Avg: 135.31205585347038, Std: 2.590081093067741
-- Generation 39 --
Min: 129.5504839313883, Max: 144.2923424762007, Avg: 137.54003014174438, Std: 2.3489147662185395
Best individual is [16.147497336125568, 18.199133247639097, 18.530466738053317, 14.238184371788577, 11.022451884629906, 11.344285227423816, 12.731661935485569, 13.755150647187078, 16.52025636666822, 11.803254721199565], (144.2923424762007,)

この実行結果は、遺伝的アルゴリズムの進化過程を示しています。

各世代の統計値が出力されており、世代ごとに次の情報が記載されています:

  • Min: その世代における個体の最小評価値
  • Max: その世代における個体の最大評価値
  • Avg: その世代における個体の平均評価値
  • Std: その世代における評価値の標準偏差

以下に、各世代における主な傾向を説明します:

  1. 初期世代 (Generation $0$):

    • 評価値は広範囲に分布しており、最小値は$-23.62$、最大値は$63.23$です。
    • 平均値は$13.95$で、標準偏差は$14.12$と高めで、個体の多様性が高いことを示しています。
  2. 中間世代 (Generation $1~20$):

    • 評価値の最小値と最大値は上昇傾向にあり、個体の適応度が徐々に向上しています。
    • 平均値も一貫して上昇しており、遺伝的アルゴリズムが効果的に良い個体を選択・交配していることを示しています。
    • 標準偏差は徐々に減少しており、個体の評価値が集中的になっていることを示しています。
  3. 後期世代 (Generation $21~39$):

    • 評価値の最小値と最大値がさらに上昇し、最大値はGeneration $39$で$144.29$に達しています。
    • 平均値も上昇し続けており、最終的には$137.54$に達しています。
    • 標準偏差は世代によって若干の変動はあるものの、全体的には減少傾向にあります。

最良の個体:

  • 最良の個体は、評価値144.29を達成した個体であり、具体的な遺伝子(パラメータ)の値は以下の通りです:
    $ [16.15, 18.20, 18.53, 14.24, 11.02, 11.34, 12.73, 13.76, 16.52, 11.80] $

この結果から、遺伝的アルゴリズムが適切に機能し、世代を重ねるごとに適応度の高い個体が生まれていることが分かります。

評価値が安定的に向上しているため、アルゴリズムが良好に収束していると考えられます。

NetworkXのいろいろな使い方

NetworkX

NetworkXは、複雑なネットワークやグラフを扱うための強力なPythonライブラリです。

ここでは、NetworkXの便利な使い方や基本的な操作方法を紹介します。

1. 基本的なグラフ操作

グラフの作成

1
2
3
4
5
6
7
import networkx as nx

# 無向グラフの作成
G = nx.Graph()

# 有向グラフの作成
D = nx.DiGraph()

ノードとエッジの追加

1
2
3
4
5
6
7
# ノードの追加
G.add_node(1)
G.add_nodes_from([2, 3])

# エッジの追加
G.add_edge(1, 2)
G.add_edges_from([(2, 3), (3, 1)])

ノードとエッジの属性設定

1
2
3
4
5
# ノード属性の設定
G.nodes[1]['color'] = 'red'

# エッジ属性の設定
G.edges[1, 2]['weight'] = 4.7

2. グラフの描画

簡単なグラフ描画

1
2
3
4
import matplotlib.pyplot as plt

nx.draw(G, with_labels=True)
plt.show()

[実行結果]


属性に基づいたグラフ描画

1
2
3
4
# ノードの色を属性に基づいて設定
colors = [G.nodes[n]['color'] if 'color' in G.nodes[n] else 'blue' for n in G.nodes]
nx.draw(G, node_color=colors, with_labels=True)
plt.show()

[実行結果]

3. グラフアルゴリズムの使用

最短経路

1
2
3
# 最短経路の計算
shortest_path = nx.shortest_path(G, source=1, target=3)
print(shortest_path)

[実行結果]

[1, 3]

クラスタリング係数

1
2
3
# クラスタリング係数の計算
clustering = nx.clustering(G)
print(clustering)

[実行結果]

{1: 1.0, 2: 1.0, 3: 1.0}

連結成分

1
2
3
# 連結成分の計算
components = list(nx.connected_components(G))
print(components)

[実行結果]

[{1, 2, 3}]

4. ネットワークの生成

様々なタイプのグラフ生成

1
2
3
4
5
6
7
8
# 完全グラフ
K = nx.complete_graph(5)

# エルデシュ・レーニィランダムグラフ
ER = nx.erdos_renyi_graph(100, 0.15)

# バラバシ・アルバートモデルのスケールフリーネットワーク
BA = nx.barabasi_albert_graph(100, 5)

5. データの入出力

グラフの保存

1
2
# グラフをファイルに保存
nx.write_gml(G, "graph.gml")

グラフの読み込み

1
2
# ファイルからグラフを読み込み
G = nx.read_gml("graph.gml")

6. 実際のデータセットを使ったグラフ解析

NetworkXには実際のデータセットを使ってグラフ解析を行うためのサポートも豊富です。

例えば、Facebookの友人ネットワークやTwitterのフォロワーネットワークなどを解析する際に利用できます。

例: Karate Club Graph

1
2
3
4
5
6
7
8
9
10
11
12
# Karate Club Graphの読み込み
G = nx.karate_club_graph()

# 基本情報の表示
print(f"Number of nodes: {G.number_of_nodes()}")
print(f"Number of edges: {G.number_of_edges()}")
print(f"Nodes: {list(G.nodes(data=True))}")
print(f"Edges: {list(G.edges(data=True))}")

# 描画
nx.draw(G, with_labels=True)
plt.show()

[実行結果]


NetworkXは、研究や実際のアプリケーションでネットワーク解析を行うために非常に便利で柔軟性の高いツールです。

上記の基本的な使い方を参考に、さまざまなネットワーク解析を試してみてください。

CVXPYのいろいろな使い方

CVXPY

CVXPYは、Python最適化問題を定義して解くための強力なライブラリです。

以下に、CVXPYの基本的な使い方や便利な機能をいくつか紹介します。

1. インストール

まず、CVXPYをインストールします。

1
pip install cvxpy

2. 基本的な使い方

線形最適化問題を定義して解く基本的な例を示します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import cvxpy as cp

# 変数の定義
x = cp.Variable()
y = cp.Variable()

# 目的関数の定義
objective = cp.Maximize(x + y)

# 制約条件の定義
constraints = [
x + 2*y <= 1,
x - y >= 1,
x >= 0,
y >= 0
]

# 問題の定義
problem = cp.Problem(objective, constraints)

# 問題の解決
problem.solve()

# 結果の表示
print("Optimal value:", problem.value)
print("x:", x.value)
print("y:", y.value)

[実行結果]

Optimal value: 1.0000000001061882
x: 1.0000000002109781
y: -1.0479004738088337e-10

3. 二次最適化

CVXPYを使って二次最適化問題を解く例です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
import numpy as np

# 変数の定義
x = cp.Variable(2)

# 目的関数の定義 (二次関数)
Q = np.array([[2, 0.5], [0.5, 1]])
objective = cp.Minimize(0.5 * cp.quad_form(x, Q))

# 制約条件の定義
constraints = [x[0] + x[1] == 1, x >= 0]

# 問題の定義
problem = cp.Problem(objective, constraints)

# 問題の解決
problem.solve()

# 結果の表示
print("Optimal value:", problem.value)
print("x:", x.value)

[実行結果]

Optimal value: 0.43750000000000006
x: [0.25 0.75]

4. 線形計画問題 (Linear Programming, LP)

線形計画問題を定義して解く例です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# 変数の定義
x = cp.Variable(3)

# 目的関数の定義
c = np.array([1, 2, 3])
objective = cp.Minimize(c.T @ x)

# 制約条件の定義
A = np.array([[1, 1, 0], [0, 1, 1]])
b = np.array([1, 1])
constraints = [A @ x <= b, x >= 0]

# 問題の定義
problem = cp.Problem(objective, constraints)

# 問題の解決
problem.solve()

# 結果の表示
print("Optimal value:", problem.value)
print("x:", x.value)

[実行結果]

Optimal value: 1.3466953587360049e-09
x: [5.70962945e-10 4.91615839e-11 2.25803082e-10]

5. ソフトな制約条件

ペナルティ付きで制約条件を緩和する方法です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# 変数の定義
x = cp.Variable()

# 制約条件に対するペナルティの定義
penalty = cp.Variable(pos=True)
constraints = [x + penalty >= 1, penalty >= 0]

# 目的関数の定義
objective = cp.Minimize(x + 10 * penalty)

# 問題の定義
problem = cp.Problem(objective, constraints)

# 問題の解決
problem.solve()

# 結果の表示
print("Optimal value:", problem.value)
print("x:", x.value)
print("penalty:", penalty.value)

[実行結果]

Optimal value: 1.0000000000163711
x: 1.0000000000163711
penalty: 0.0

6. パラメータの使用

パラメータを使って問題の定義を柔軟に変更する方法です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# 変数とパラメータの定義
x = cp.Variable()
a = cp.Parameter(nonneg=True)
b = cp.Parameter(nonneg=True)

# 目的関数の定義
objective = cp.Minimize(a * x + b)

# 制約条件の定義
constraints = [x >= 0, x <= 10]

# 問題の定義
problem = cp.Problem(objective, constraints)

# パラメータの値を設定して問題を解く
a.value = 2
b.value = 1
problem.solve()

print("Optimal value with a=2, b=1:", problem.value)
print("x:", x.value)

# パラメータの値を変更して再度問題を解く
a.value = 1
b.value = 3
problem.solve()

print("Optimal value with a=1, b=3:", problem.value)
print("x:", x.value)

[実行結果]

Optimal value with a=2, b=1: 1.0000000003054508
x: 1.527254291082821e-10
Optimal value with a=1, b=3: 3.0000000003691296
x: 3.691294636403306e-10

7. ソルバーの指定

特定のソルバーを使用する方法です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 変数の定義
x = cp.Variable()

# 目的関数の定義
objective = cp.Minimize(x**2)

# 制約条件の定義
constraints = [x >= 1]

# 問題の定義
problem = cp.Problem(objective, constraints)

# 特定のソルバーを指定して問題を解く
problem.solve(solver=cp.SCS)

print("Optimal value using SCS solver:", problem.value)
print("x:", x.value)

[実行結果]

Optimal value using SCS solver: 0.9999999999272082
x: 0.9999999999636041

これらの基本的な操作を理解することで、CVXPYを使って幅広い最適化問題を効率的に解くことができます。

3Dグラフ

3Dグラフ

以下の方程式を使用して3Dグラフを描いてみます:

$$
z = \frac{\sin(x^2 + y^2)}{x^2 + y^2}
$$

この方程式は、中心が原点で放射状に渦巻くような形を描くことが知られています。

これをMatplotlibを用いて描画してみましょう。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

# データの生成
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
x, y = np.meshgrid(x, y)
z = np.sin(x**2 + y**2) / (x**2 + y**2)

# プロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z, cmap='viridis')

# ラベルとタイトルの設定
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
ax.set_title('3D Plot of $z = \\frac{\\sin(x^2 + y^2)}{x^2 + y^2}$')

plt.show()

このコードでは、$ ( z = \frac{\sin(x^2 + y^2)}{x^2 + y^2} ) $の関数を使用して、$x$と$y$の範囲で$z$を計算し、それを3D表面プロットとして描画しています。

Matplotlibplot_surface メソッドで表面を描画し、軸ラベルとタイトルも設定しています。

このグラフでは、原点を中心にして放射状に波打つ特徴的な形状が描かれることが期待されます。

[実行結果]

他にもさまざまな数学的な関数や方程式を試して、興味深い3Dグラフを作成してみてください。

ソースコード解説

ソースコードの説明をします。

1. ライブラリのインポート

1
2
3
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
  • matplotlib.pyplot: Matplotlibのプロット機能を利用するためのモジュール。
    pltとしてインポートされています。
  • mpl_toolkits.mplot3d.Axes3D: 3次元のプロットを扱うためのクラスが含まれています。
  • numpy: 数値計算を行うための基本的なモジュール。
    npとしてインポートされています。

2. データの生成

1
2
3
4
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
x, y = np.meshgrid(x, y)
z = np.sin(x**2 + y**2) / (x**2 + y**2)
  • np.linspace(-5, 5, 100): $-5$から$5$までの範囲を等間隔に区切った$100$個の数値を生成し、それらが$x$および$y$の値として使用されます。
  • np.meshgrid(x, y): xとyの配列から格子状の座標を生成します。
    これにより、$x$と$y$の組み合わせがすべて生成され、それに対応する$z$値が計算されます。
  • z = np.sin(x**2 + y**2) / (x**2 + y**2): 指定された数学的な関数を使用して$z$値を計算します。
    ここでは$ ( z = \frac{\sin(x^2 + y^2)}{x^2 + y^2} ) $を計算しています。

3. プロットの設定と描画

1
2
3
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z, cmap='viridis')
  • plt.figure(): 新しい図を生成します。
  • fig.add_subplot(111, projection='3d'): 3D描画用のサブプロットを設定します。
    ここで111は1行1列の1番目の位置を意味し、projection='3d'で3D描画を有効にします。
  • ax.plot_surface(x, y, z, cmap='viridis'): $x$, $y$, $z$の値を使用して3D表面を描画します。
    cmap='viridis'はカラーマップを指定しています。

4. ラベルとタイトルの設定

1
2
3
4
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
ax.set_title('3D Plot of $z = \\frac{\\sin(x^2 + y^2)}{x^2 + y^2}$')
  • ax.set_xlabel('X Label'), ax.set_ylabel('Y Label'), ax.set_zlabel('Z Label'): $x$軸、$y$軸、$z$軸のラベルを設定します。
  • ax.set_title('3D Plot of $z = \\frac{\\sin(x^2 + y^2)}{x^2 + y^2}$'): プロットのタイトルを設定します。
    ここではLaTeXの数式を用いてzの式をタイトルに表示しています。

5. プロットの表示

1
plt.show()
  • plt.show(): 最後にプロットを表示します。

このようにして、与えられた数学的な関数を基にして、3Dの表面プロットをMatplotlibを使って作成することができます。

産業プロセス最適化

産業プロセス最適化

産業プロセス最適化の例題として、混合飲料工場の原材料配合最適化問題を考えてみましょう。

この例では、原材料$A$、$B$、$C$を混ぜて飲料を作る際に、コストを最小化しつつ、必要な栄養成分を満たす配合を見つける問題を扱います。

問題設定

  1. 原材料の単価:

    • 原材料A: $3円/g$
    • 原材料B: $2円/g$
    • 原材料C: $1円/g$
  2. 栄養成分(各$100g$あたり):

    • 原材料A: タンパク質$10g$、脂質$5g$、炭水化物$15g$
    • 原材料B: タンパク質$5g$、脂質$10g$、炭水化物$5g$
    • 原材料C: タンパク質$2g$、脂質$4g$、炭水化物$20g$
  3. 必要な栄養成分量(最小限):

    • タンパク質: $30g$
    • 脂質: $20g$
    • 炭水化物: $50g$

目標

  • 総コストを最小化する原材料$A$、$B$、$C$の配合量を求める。

制約条件

  • 各原材料の配合量は非負($>= 0$)
  • 栄養成分量の制約を満たすこと

最適化モデル

この問題は線形計画法(Linear Programming, LP)を使って解くことができます。

Pythonでは scipy.optimize.linprog を用いることで簡単に解くことができます。

コード例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
import numpy as np
from scipy.optimize import linprog

# 原材料の単価
cost = np.array([3, 2, 1])

# 栄養成分マトリックス (各行が原材料、各列が栄養成分: [タンパク質, 脂質, 炭水化物])
nutrition = np.array([
[10, 5, 15], # 原材料A
[5, 10, 5], # 原材料B
[2, 4, 20] # 原材料C
])

# 必要な栄養成分量
requirements = np.array([30, 20, 50])

# 線形計画法の実行
result = linprog(c=cost, A_ub=-nutrition.T, b_ub=-requirements, bounds=(0, None), method='highs')

# 結果の表示
if result.success:
print("最適な配合量:")
print(f"原材料A: {result.x[0]:.2f} g")
print(f"原材料B: {result.x[1]:.2f} g")
print(f"原材料C: {result.x[2]:.2f} g")
print(f"総コスト: {result.fun:.2f} 円")
else:
print("解が見つかりませんでした")

実行結果の解釈

このコードを実行すると、原材料$A$、$B$、$C$の最適な配合量とそのときの総コストが表示されます。

これにより、工場が最もコスト効率の良い方法で必要な栄養成分を満たすことができます。

最適化の結果は以下の通りです:

最適な配合量

  • 原材料A: $2.67 g$
  • 原材料B: $0.52 g$
  • 原材料C: $0.37 g$

総コスト

  • $9.41 円$

これにより、工場は総コストを$9.41$円に抑えつつ、必要な栄養成分(タンパク質$30g$、脂質$20g$、炭水化物$50g$)を満たすことができます。

この配合が最適な解となります。