トラフィックエンジニアリング

問題:

ISP（インターネットサービスプロバイダー）は、特定の時間帯におけるネットワークトラフィックを効率的に管理し、ボトルネックを防ぎたいと考えています。
各リンクのトラフィック容量と予測されるトラフィック量が与えられたとき、トラフィックをどのように分散すれば全体の遅延を最小化できるかを求めてください。

入力:

• ノードとエッジのリスト。
• 各エッジのトラフィック容量。
• 各ノード間の予測トラフィック量。

出力:

• トラフィックの分配方法（各経路にどれだけのトラフィックを流すか）。

解法:

この問題は「多品種フロー問題」として知られ、線形計画法やネットワークシンプレックス法を用いて解決できます。

これらの例題は、通信網最適化の典型的な問題であり、現実のネットワーク設計や管理において非常に重要です。
各問題に対する解法は、適切なアルゴリズムとデータ構造を利用して実装できます。

例題

具体例として、簡単なトラフィックエンジニアリング問題を設定し、Pythonで解決し、グラフ化します。

具体例では、リンクのトラフィック容量と予測されるトラフィック量を考慮して、トラフィックを最適に分配します。

具体例

以下のようなネットワークがあるとします：

• ノード: A, B, C, D
• エッジとその容量:
  • A-B: $10$
  • A-C: $15$
  • B-D: $10$
  • C-D: $10$
  • B-C: $5$

予測されるトラフィック量は以下の通りです：

• AからDへのトラフィック: $10$
• BからCへのトラフィック: $5$

解法

この問題を線形計画法を使って解決します。

Pythonのscipy.optimize.linprogを使用します。

Pythonコード

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

1

pip install scipy networkx matplotlib

次に、トラフィックエンジニアリングの具体例を解くPythonコードを示します。

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
62
63
64
65
66
67

import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
from scipy.optimize import linprog

# 定義されたエッジと容量
edges = [('A', 'B', 10), ('A', 'C', 15), ('B', 'D', 10), ('C', 'D', 10), ('B', 'C', 5)]
nodes = set(sum(([u, v] for u, v, _ in edges), []))
node_list = sorted(nodes)

# ノード間の予測トラフィック量
traffic = {('A', 'D'): 10, ('B', 'C'): 5}

# エッジの数とノードの数
num_edges = len(edges)
num_nodes = len(node_list)

# 線形計画法のセットアップ
c = np.ones(num_edges)  # 最小化するコスト関数の係数（単純に全て1）
A_eq = np.zeros((num_nodes, num_edges))
b_eq = np.zeros(num_nodes)

# 各エッジのインデックスを保持
edge_indices = {edge[:2]: idx for idx, edge in enumerate(edges)}

# ノードのインデックスを保持
node_indices = {node: idx for idx, node in enumerate(node_list)}

# 各エッジの制約をセットアップ
for (u, v, capacity) in edges:
    idx = edge_indices[(u, v)]
    A_eq[node_indices[u], idx] = 1
    A_eq[node_indices[v], idx] = -1

# トラフィックの制約をセットアップ
for (src, dst), demand in traffic.items():
    b_eq[node_indices[src]] += demand
    b_eq[node_indices[dst]] -= demand

# 各エッジの容量制約をセットアップ
bounds = [(0, capacity) for u, v, capacity in edges]

# 線形計画法を解く
res = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs')

# 結果の表示
if res.success:
    print("Optimal traffic distribution:")
    for (u, v, capacity), flow in zip(edges, res.x):
        print(f"Flow on edge {u}-{v}: {flow:.2f}")
else:
    print("Optimization failed.")

# グラフの描画
G = nx.DiGraph()
for (u, v, capacity), flow in zip(edges, res.x):
    G.add_edge(u, v, capacity=capacity, flow=flow)

pos = nx.spring_layout(G)
edge_labels = {(u, v): f"{flow:.2f}/{capacity}" for (u, v, capacity), flow in zip(edges, res.x)}

plt.figure(figsize=(10, 8))
nx.draw(G, pos, with_labels=True, node_size=3000, node_color="lightblue", font_size=16, font_weight="bold", arrowsize=20)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, font_size=14)
nx.draw_networkx_edges(G, pos, width=2, edge_color='black', arrows=True)
plt.title("Optimized Traffic Distribution in the Network")
plt.show()

ソースコード解説

1. ライブラリのインポート: numpy, networkx, matplotlib, scipy.optimizeをインポートします。
2. ネットワークの定義: ノード、エッジ、容量、予測トラフィック量を定義します。
3. 線形計画法のセットアップ: コスト関数、制約行列、制約ベクトル、境界を設定します。
4. 線形計画法の解決: linprogを使用して問題を解きます。
5. 結果の表示: 各エッジに対する最適なトラフィック分配を表示します。
6. グラフの描画: networkxを使用してネットワークグラフを描画し、各エッジのフローを表示します。

このコードを実行すると、最適なトラフィック分配が計算され、その結果がグラフとして視覚化されます。

結果解説

表示される結果は、最適なトラフィック分配を示しています。

[実行結果]

以下に各エッジのフローについて説明します。

フロー結果

1. Flow on edge A-B: 5.00

   • AからBへのフローが$5.00$です。
     これは、BからCへのトラフィック量を処理するためにAからBに送られるトラフィックです。

2. Flow on edge A-C: 5.00

   • AからCへのフローが$5.00$です。
     これは、AからDへのトラフィック量の一部がAからCを経由してDに送られることを示しています。

3. Flow on edge B-D: 10.00

   • BからDへのフローが$10.00$です。
     これは、AからDへのトラフィック量がA-B-D経由で送られていることを示しています。
     つまり、AからBに送られたトラフィックのうち、$5.00$はBからDに直接送られています。
     また、B-Cへのトラフィック量も含まれています。

4. Flow on edge C-D: 0.00

   • CからDへのフローは$0.00$です。
     これは、AからDへのトラフィックがすべてA-B-D経由で送られているため、C-D経由でのトラフィックは発生していないことを示しています。

5. Flow on edge B-C: 0.00

   • BからCへのフローは$0.00$です。
     BからCへのトラフィック量$5.00$は、AからCを経由して送られているため、B-C間にはトラフィックが発生していません。

全体の流れ

• AからDへのトラフィック（10）:

  • AからBを経由して$5.00$が送られ、残りの$5.00$はAからCを経由して送られています。
  • BからDへのフローは$10.00$で、これはAからDへのトラフィック全体を運ぶために使われています。

• BからCへのトラフィック（5）:

  • AからBに送られた$5.00$がBからDに送られていますが、最終的にはAからCを経由してCに到達する形で処理されています。

この結果は、全体のネットワークトラフィックが最適に分配され、エッジの容量制限内でトラフィックが処理されていることを示しています。

二次計画法（Quadratic Programming, QP）

CVXPYを使用して二次計画法（Quadratic Programming, QP）の問題を解く方法を説明します。

二次計画法の一般的な形式は以下のようになります。

目的関数：
$$
\text{minimize} \quad \frac{1}{2} x^T P x + q^T x
$$

制約：
$$
Gx \leq h
$$
$$
Ax = b
$$

ここで、$(x)$は決定変数のベクトル、$(P)$は正定値対称行列、$(q)$はベクトル、$(G)$、$(h)$、$(A)$、および$(b)$はそれぞれ行列またはベクトルです。

プログラム例

以下はPythonのコードで、CVXPYを使用してこのような問題を解く方法の一例です。

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

import cvxpy as cp
import numpy as np

# 定義する二次計画法の問題
P = np.array([[2, 0], [0, 2]])    # 目的関数の2次項の係数
q = np.array([-2, -5])            # 目的関数の1次項の係数
G = np.array([[1, 2], [1, -4], [-1, 0], [0, -1]]) # 不等式制約の係数
h = np.array([3, 3, 0, 0])        # 不等式制約の右辺
A = np.array([[1, 1]])            # 等式制約の係数
b = np.array([1])                 # 等式制約の右辺

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

# 目的関数を定義
objective = cp.Minimize((1/2) * cp.quad_form(x, P) + q @ x)

# 制約を定義
constraints = [G @ x <= h, A @ x == b]

# 問題を定義
prob = cp.Problem(objective, constraints)

# 問題を解く
result = prob.solve()

# 結果を出力
print(f"最適値: {result}")
print(f"最適解: {x.value}")

説明

1. 行列およびベクトルの定義:

   • P、q、G、h、A、b それぞれ問題の行列およびベクトルを定義します。

2. 変数の定義:

   • x = cp.Variable(2) で、最適化変数 $( x ) $を定義します。
     ここでは2次元の変数を考えています。

3. 目的関数の設定:

   • objective = cp.Minimize((1/2) * cp.quad_form(x, P) + q @ x) で、目的関数を定義しています。

4. 制約の設定:

   • constraints = [G @ x <= h, A @ x == b] で不等式制約と等式制約を定義しています。

5. 問題の解決:

   • prob = cp.Problem(objective, constraints) で問題を作成し、result = prob.solve() で解を求めます。

6. 結果の出力:

   • 最適値と最適解を出力します。

このようにして、CVXPYを使用して二次計画法の問題をシンプルに解くことができます。

CVXPYは非常に強力であり、もっと複雑な制約や目的関数にも対応することができます。

[実行結果]

最適値: -4.0
最適解: [1.78280298e-24 1.00000000e+00]

結果解説

CVXPYを使用した二次計画法の解の結果について説明します。

解の結果

• 最適値 (Objective Value): $-4.0$
• 最適解 (Optimal Solution): $[1.78280298e-24, 1.00000000e+00]$

1. 最適値 (Objective Value)

最適値: -4.0 という結果は、設定した二次計画法の目的関数を最小化したときの最小値が$ -4.0 $になることを示しています。
つまり、与えられた制約の範囲内で目的関数の評価が最も低くなる点（最適解）での値が$ -4.0 $です。

2. 最適解 (Optimal Solution)

最適解: [1.78280298e-24, 1.00000000e+00] という結果は、$( x ) $の最適値が$ [1.78280298e-24, 1.00000000e+00] $であることを示しています。
これが、与えられた目的関数を最小にする$ ( x ) $の値です。

結果の詳細

• 最適解$ ( x )$:
  • $( x_1 ) = 1.78280298e-24$
  • $( x_2 ) = 1.00000000e+00$

ここで、$( x_1 ) $が非常に小さい値（ほぼゼロ）になっていますが、これは最適解の計算過程での丸め誤差や数値計算の性質からくるものです。

実質的には$ ( x_1 \approx 0 ) $と解釈しても問題ありません。

制約の確認

この結果が制約を満たしていることを確認してみましょう。

【不等式制約】
$[ Gx \leq h ]$

ここで$ ( G ) $と$ ( h ) $は次のように定義されます。
$[ G = \begin{pmatrix} 1 & 2 \ 1 & -4 \ -1 & 0 \ 0 & -1 \end{pmatrix}, \quad h = \begin{pmatrix} 3 \ 3 \ 0 \ 0 \end{pmatrix} ]$

最適解$ ( x ) $が $[0, 1] $付近であることを考慮します。

• 第1の制約: $( 1 \cdot 0 + 2 \cdot 1 \leq 3 ) $ は$ ( 2 \leq 3 ) $で成立。
• 第2の制約: $( 1 \cdot 0 - 4 \cdot 1 \leq 3 ) $は$ ( -4 \leq 3 ) $で成立。
• 第3の制約: $( -1 \cdot 0 \leq 0 ) $は$ ( 0 \leq 0 ) $で成立。
• 第4の制約: $( 0 - 1 \cdot 1 \leq 0 ) $は$ ( -1 \leq 0 ) $で成立。

これらの不等式制約はすべて満たされています。

【等式制約】
$[ Ax = b ]$

ここで $( A ) $と$ ( b ) $は次のように定義されます。
$[ A = \begin{pmatrix} 1 & 1 \end{pmatrix}, \quad b = 1 ]$

• 制約: $( 1 \cdot 0 + 1 \cdot 1 = 1 )$

こちらの等式制約も満たされています。

まとめ

今回求めた最適解と最適値は、与えられた二次計画法の目的関数および制約条件に対して正しく計算されています。

最適解$ [1.78280298e-24, 1.00000000e+00] $が正しい最適解であり、これにより目的関数の最小値 $ -4.0 $が得られました。

• 最適値 (目的関数の最小値) は$ -4.0$。
• 最適解は$ ( x_1 \approx 0 )$、$( x_2 = 1 )$。

バーンズリーのシダ（Barnsley's fern）

バーンズリーのシダ（Barnsley’s fern）は、確率的な方法で生成されるフラクタル構造の一種であり、数学的な美しさと興味深さで知られています。

このシダは、数学者であるMichael Barnsleyによって導入されました。

バーンズリーのシダは、次のような確率的な方法で生成されます。

まず、初期条件として $(0, 0) $の点を選びます。

その後、一連の確率的な変換を適用して新しい点を生成します。

これらの変換は、特定の確率で選択され、シダの形を生成します。

バーンズリーのシダの生成には、以下のような変換が使用されます：

1. $1%$の確率で、点$ (x, y) $は$ (0, 0.16y) $に変換されます。
   これにより、葉の部分が生成されます。

2. $85%$の確率で、点$ (x, y) $は$ (0.85x + 0.04y, -0.04x + 0.85y + 1.6) $に変換されます。
   これにより、茎と枝が生成されます。

3. $7%$の確率で、点$ (x, y) $は$ (0.20x - 0.26y, 0.23x + 0.22y + 1.6) $に変換されます。
   これにより、さらに複雑な枝が生成されます。

4. $7%$の確率で、点$ (x, y) $は$ (-0.15x + 0.28y, 0.26x + 0.24y + 0.44) $に変換されます。
   これにより、より細かい枝が生成されます。

これらの変換を多数回繰り返すことで、バーンズリーのシダの構造が形成されます。

このフラクタル構造は、自然界に見られるシダの葉や枝の形状に類似しており、その美しさと複雑さから数学的な芸術の一形態として広く愛されています。

プログラム例

以下は、NetworkXを使用してバーンズリーのシダを生成し、グラフで表現するサンプルコードです。

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

import random
import networkx as nx
import matplotlib.pyplot as plt

# バーンズリーのシダ生成関数
def barnsley_fern(n):
    # 初期条件
    x, y = 0, 0
    points = [(x, y)]
    
    # バーンズリーのシダ生成
    for _ in range(n):
        r = random.random()
        if r <= 0.01:
            x, y = 0, 0.16 * y
        elif r <= 0.86:
            x, y = 0.85 * x + 0.04 * y, -0.04 * x + 0.85 * y + 1.6
        elif r <= 0.93:
            x, y = 0.20 * x - 0.26 * y, 0.23 * x + 0.22 * y + 1.6
        else:
            x, y = -0.15 * x + 0.28 * y, 0.26 * x + 0.24 * y + 0.44
        points.append((x, y))
    
    return points

# バーンズリーのシダを生成
fern_points = barnsley_fern(10000)

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

# ノードの追加
for i, (x, y) in enumerate(fern_points):
    G.add_node(i, pos=(x, y))

# エッジの追加
for i in range(len(fern_points) - 1):
    G.add_edge(i, i+1)

# グラフの描画
pos = nx.get_node_attributes(G, 'pos')
nx.draw(G, pos, node_size=0, edge_color='green', alpha=0.5, linewidths=0.5)
plt.title("Barnsley Fern using NetworkX")
plt.show()

このコードでは、バーンズリーのシダを生成するための関数 barnsley_fern を定義し、その結果をNetworkXのグラフに変換して描画しています。

生成されたグラフは、バーンズリーのシダの形を模倣した美しい構造を持っています。

[実行結果]

ソースコード解説

このソースコードは、Pythonを使用してバーンズリーのシダ（Barnsley’s fern）を生成し、NetworkXを使用してそのグラフを描画するものです。

以下では、コードの構造と各部分の役割について詳しく説明します。

1. バーンズリーのシダ生成関数の定義

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

def barnsley_fern(n):
    x, y = 0, 0  # 初期条件
    points = [(x, y)]  # 座標を格納するリスト

    # バーンズリーのシダ生成
    for _ in range(n):
        r = random.random()  # 0から1の乱数を生成
        if r <= 0.01:
            x, y = 0, 0.16 * y
        elif r <= 0.86:
            x, y = 0.85 * x + 0.04 * y, -0.04 * x + 0.85 * y + 1.6
        elif r <= 0.93:
            x, y = 0.20 * x - 0.26 * y, 0.23 * x + 0.22 * y + 1.6
        else:
            x, y = -0.15 * x + 0.28 * y, 0.26 * x + 0.24 * y + 0.44
        points.append((x, y))  # 新しい座標をリストに追加

    return points

この部分では、バーンズリーのシダを生成する関数 barnsley_fern が定義されています。

この関数は、生成する点の数 n を引数として受け取ります。

バーンズリーのシダは、一連の確率的な変換を繰り返すことで生成されます。

各ステップで、ランダムに選択された変換に基づいて新しい座標が計算され、リスト points に追加されます。

2. バーンズリーのシダを生成

1

fern_points = barnsley_fern(10000)

この行では、barnsley_fern 関数を呼び出して、$10000$個の点からなるバーンズリーのシダを生成し、fern_points に格納しています。

3. グラフの作成

1

G = nx.Graph()

この行では、NetworkXのグラフオブジェクト G を作成しています。

4. ノードの追加

1
2

for i, (x, y) in enumerate(fern_points):
    G.add_node(i, pos=(x, y))

この部分では、バーンズリーのシダで生成された点をグラフのノードとして追加しています。

各ノードには、座標を属性として持つことができます。

ここでは、各ノードに座標情報を pos という名前の属性として追加しています。

5. エッジの追加

1
2

for i in range(len(fern_points) - 1):
    G.add_edge(i, i+1)

この部分では、連続する点の間にエッジを追加しています。

これにより、シダの枝の構造が表現されます。

6. グラフの描画

1
2
3
4

pos = nx.get_node_attributes(G, 'pos')
nx.draw(G, pos, node_size=0, edge_color='green', alpha=0.5, linewidths=0.5)
plt.title("Barnsley Fern using NetworkX")
plt.show()

最後に、NetworkXの draw 関数を使用してグラフを描画しています。

ノードの座標は pos から取得し、それに基づいてノードとエッジが描画されます。

結果として得られるグラフは、バーンズリーのシダの形状を美しく表現しています。

輸送問題 (Transportation Problem)

最適化問題の例として、輸送問題 (Transportation Problem) を考えます。

この問題は、複数の供給地から複数の需要地に商品を輸送する際の輸送コストを最小化する問題です。

以下は簡単な輸送問題の例です。

問題設定:

• 供給地A、B、Cがあり、それぞれの供給量は次の通りです:
  • A: 20
  • B: 30
  • C: 25
• 需要地1、2、3があり、それぞれの需要量は次の通りです:
  • 1: 30
  • 2: 35
  • 3: 10
• 各供給地から各需要地への輸送コストは次のように設定されています:

この輸送問題を PuLP で解いてみましょう。

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

import pulp

# 供給量
supply = {'A': 20, 'B': 30, 'C': 25}

# 需要量
demand = {1: 30, 2: 35, 3: 10}

# 輸送コスト
cost = {'A': {1: 8, 2: 6, 3: 10},
        'B': {1: 9, 2: 12, 3: 13},
        'C': {1: 14, 2: 9, 3: 16}}

# 問題の定義
prob = pulp.LpProblem("Transportation_Problem", pulp.LpMinimize)

# 変数の定義
routes = [(i,j) for i in supply for j in demand]
vars = pulp.LpVariable.dicts("Route", (supply, demand), lowBound=0, cat='Continuous')

# 目的関数
prob += pulp.lpSum([vars[i][j] * cost[i][j] for (i, j) in routes])

# 制約条件
# 各供給地から出る供給量がその供給地の供給量に等しい
for i in supply:
    prob += pulp.lpSum([vars[i][j] for j in demand]) == supply[i]

# 各需要地に入る供給量がその需要地の需要量に等しい
for j in demand:
    prob += pulp.lpSum([vars[i][j] for i in supply]) == demand[j]

# 最適化の実行
prob.solve()

# 結果の表示
print("Status:", pulp.LpStatus[prob.status])
print("Total Cost:", pulp.value(prob.objective))
for v in prob.variables():
    print(v.name, "=", v.varValue)

このコードでは、供給地と需要地の間の輸送量を決定し、全体の輸送コストを最小化するように最適化しています。

最適化の結果に基づいて、どの供給地からどの需要地にいくらの商品を輸送すればコストが最小になるかを表示します。

[実行結果]

Status: Optimal
Total Cost: 655.0
Route_A_1 = 0.0
Route_A_2 = 10.0
Route_A_3 = 10.0
Route_B_1 = 30.0
Route_B_2 = 0.0
Route_B_3 = 0.0
Route_C_1 = 0.0
Route_C_2 = 25.0
Route_C_3 = 0.0

ソースコード解説

ソースコードを詳しく解説します。

1. 概要

このプログラムは、Pulpライブラリを使用して輸送問題を最適化する例です。

輸送問題は、供給地から需要地へ物資を輸送する際のコストを最小化するために、各ルートに割り当てる物資の量を決定する問題です。

2. 必要なライブラリのインポート

1

import pulp

• Pulpは、線形計画問題（LP）や整数計画問題（ILP）を解くためのPythonライブラリです。

3. データの定義

1
2
3
4
5
6
7
8
9
10

# 供給量
supply = {'A': 20, 'B': 30, 'C': 25}

# 需要量
demand = {1: 30, 2: 35, 3: 10}

# 輸送コスト
cost = {'A': {1: 8, 2: 6, 3: 10},
        'B': {1: 9, 2: 12, 3: 13},
        'C': {1: 14, 2: 9, 3: 16}}

• 供給量: 各供給地の供給量を表します。
• 需要量: 各需要地の需要量を表します。
• 輸送コスト: 各供給地から需要地への輸送コストを表します。

4. 問題の定義

1

prob = pulp.LpProblem("Transportation_Problem", pulp.LpMinimize)

• 輸送問題を最小化問題として定義します。ここでは、輸送コストを最小化することが目的です。

5. 変数の定義

1
2

routes = [(i,j) for i in supply for j in demand]
vars = pulp.LpVariable.dicts("Route", (supply, demand), lowBound=0, cat='Continuous')

• routes: 全ての供給地と需要地の組み合わせをリストにします。
• vars: 各ルートにおける輸送量を表す変数を定義します。ここでは、連続変数として定義しています。

6. 目的関数の定義

1

prob += pulp.lpSum([vars[i][j] * cost[i][j] for (i, j) in routes])

• 目的関数は、全てのルートにおける輸送コストの総和を最小化することです。

7. 制約条件の定義

1
2
3
4
5
6
7

# 各供給地から出る供給量がその供給地の供給量に等しい
for i in supply:
    prob += pulp.lpSum([vars[i][j] for j in demand]) == supply[i]

# 各需要地に入る供給量がその需要地の需要量に等しい
for j in demand:
    prob += pulp.lpSum([vars[i][j] for i in supply]) == demand[j]

• 供給制約: 各供給地から出る供給量が、その供給地の供給量に等しくなるようにします。
• 需要制約: 各需要地に入る供給量が、その需要地の需要量に等しくなるようにします。

8. 最適化の実行

1

prob.solve()

• Pulpのsolveメソッドを使用して、定義した問題を解きます。

9. 結果の表示

1
2
3
4

print("Status:", pulp.LpStatus[prob.status])
print("Total Cost:", pulp.value(prob.objective))
for v in prob.variables():
    print(v.name, "=", v.varValue)

• 最適化のステータス（成功・失敗）を表示します。
• 最小化された総輸送コストを表示します。
• 各ルートにおける最適な輸送量を表示します。

結果解説

この実行結果を詳しく解説します。

[実行結果]

Status: Optimal
Total Cost: 655.0
Route_A_1 = 0.0
Route_A_2 = 10.0
Route_A_3 = 10.0
Route_B_1 = 30.0
Route_B_2 = 0.0
Route_B_3 = 0.0
Route_C_1 = 0.0
Route_C_2 = 25.0
Route_C_3 = 0.0

結果:

• 最適化ステータス: Optimal
  • この結果は、PuLPが問題に対して最適な解を見つけたことを示しています。
• 総輸送コスト: 655.0
  • この値は、最適な輸送計画に基づいて計算された最小の総コストです。

各供給地から各需要地への輸送量:

• Route_A_1 = 0.0
  • 供給地Aから需要地1への輸送量は 0 です。
• Route_A_2 = 10.0
  • 供給地Aから需要地2への輸送量は 10 です。
• Route_A_3 = 10.0
  • 供給地Aから需要地3への輸送量は 10 です。
• Route_B_1 = 30.0
  • 供給地Bから需要地1への輸送量は 30 です。
• Route_B_2 = 0.0
  • 供給地Bから需要地2への輸送量は 0 です。
• Route_B_3 = 0.0
  • 供給地Bから需要地3への輸送量は 0 です。
• Route_C_1 = 0.0
  • 供給地Cから需要地1への輸送量は 0 です。
• Route_C_2 = 25.0
  • 供給地Cから需要地2への輸送量は 25 です。
• Route_C_3 = 0.0
  • 供給地Cから需要地3への輸送量は 0 です。

解釈

この結果から、以下の輸送プランがコストを最小化するために用いられます。

• 供給地A:
  • 需要地2に 10 単位を輸送します。 (コスト: $(6 \times 10 = 60)$)
  • 需要地3に 10 単位を輸送します。 (コスト: $(10 \times 10 = 100)$)
• 供給地B:
  • 需要地1に 30 単位を輸送します。 (コスト: $(9 \times 30 = 270)$)
• 供給地C:
  • 需要地2に 25 単位を輸送します。 (コスト: $(9 \times 25 = 225)$)

総コストの計算

上記の輸送プランに基づいた輸送コストは次のように計算できます。

• 供給地Aのコスト: $(60 + 100 = 160)$
• 供給地Bのコスト: $(270)$
• 供給地Cのコスト: $(225)$

したがって、合計の輸送コストは下記のようになります。
$$
160 + 270 + 225 = 655
$$

エンネーパーの面（Enneper Surface）

エンネーパーの面（Enneper Surface）は、微分幾何学における特異な極小曲面の一例であり、以下の特徴を持ちます：

特徴

1. 極小曲面：

   • エンネーパーの面は、曲面上の任意の点で曲率が最小となる極小曲面です。
     極小曲面は、物理的には石鹸膜のように、表面積を最小化する形状を取る面です。

2. 自己交差：

   • エンネーパーの面は、自己交差する点を持ちます。
     つまり、曲面の一部が別の部分と交差する点が存在します。

3. パラメトリック方程式：

   • エンネーパーの面は、以下のパラメトリック方程式で定義されます：
     $$
     x(u, v) = u - \frac{u^3}{3} + u v^2
     $$
     $$
     y(u, v) = v - \frac{v^3}{3} + v u^2
     $$
     $$
     z(u, v) = u^2 - v^2
     $$
     ここで、$(u) $と$ (v) $はパラメータです。

4. 対称性：

   • エンネーパーの面は、対称性を持つ曲面であり、中心対称的な形状をしています。
     曲面は中央部分から外側に向かって曲がりくねる形状を持ち、対称的に広がります。

5. 幾何学的美しさ：

   • エンネーパーの面は、その複雑で滑らかな形状から、数学的にも美しい曲面として知られています。

例

エンネーパーの面は、数学的研究や教育においてよく使用される曲面で、極小曲面の典型例として知られています。

コンピュータグラフィックスや芸術的なデザインにおいても、その独特の形状が利用されることがあります。

視覚的な説明

エンネーパーの面を視覚化すると、以下のような特徴が見られます：

• 曲面は中央部で交差し、内側にくびれた形状を形成します。
• 高さの違いによって色が変わるカラーマップを適用することで、曲面の高低差が視覚的に強調されます。

エンネーパーの面は、その独特の形状と数学的特性から、微分幾何学における重要な研究対象の一つです。

プログラム例

以下にPythonコードを示します。

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

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns

# seabornのスタイルを使用
sns.set(style='darkgrid')

# パラメータの範囲を設定
u = np.linspace(-2, 2, 100)
v = np.linspace(-2, 2, 100)
u, v = np.meshgrid(u, v)

# エンネーパーの面のパラメトリック方程式
x = u - (u**3) / 3 + u * v**2
y = v - (v**3) / 3 + v * u**2
z = u**2 - v**2

# プロットを作成
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

# 3Dサーフェスプロット
surf = ax.plot_surface(x, y, z, cmap='plasma', edgecolor='none')

# カラーバーを追加
fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)

# 軸ラベル
ax.set_xlabel('X軸')
ax.set_ylabel('Y軸')
ax.set_zlabel('Z軸')

# タイトルを追加
ax.set_title('Enneper Surface')

# グラフを表示
plt.show()

このコードを実行すると、エンネーパーの面の3Dプロットが描画されます。

学術的に興味深い形状を持つエンネーパーの面を使って、複雑な3Dグラフを表示します。

[実行結果]

ソースコード解説

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

1. 必要なライブラリのインポート

1
2
3
4

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns

このセクションでは、以下のライブラリをインポートしています。

• numpy: 数値計算を効率的に行うためのライブラリ。
• matplotlib.pyplot: グラフ作成のためのライブラリ。
• mpl_toolkits.mplot3d: 3Dプロットをサポートするためのライブラリ。
• seaborn: 見栄えの良いグラフを作成するためのライブラリ。

2. Seabornスタイルの設定

1
2

# seabornのスタイルを使用
sns.set(style='darkgrid')

Seabornのスタイルをdarkgridに設定しています。
これにより、プロットの背景がダークグリッドスタイルになります。

3. パラメータの範囲を設定

1
2
3
4

# パラメータの範囲を設定
u = np.linspace(-2, 2, 100)
v = np.linspace(-2, 2, 100)
u, v = np.meshgrid(u, v)

ここでは、uとvの範囲を-2から2まで$100$等分した配列を作成しています。
その後、np.meshgridを使用して2次元のグリッドデータを作成しています。

4. エンネーパーの面のパラメトリック方程式

1
2
3
4

# エンネーパーの面のパラメトリック方程式
x = u - (u**3) / 3 + u * v**2
y = v - (v**3) / 3 + v * u**2
z = u**2 - v**2

エンネーパーの面のパラメトリック方程式を使って、x, y, zの座標を計算しています。

5. プロットの作成

1
2
3

# プロットを作成
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

figオブジェクトを作成し、サイズを設定しています。
次に、3Dプロット用のaxオブジェクトを追加しています。

6. 3Dサーフェスプロット

1
2

# 3Dサーフェスプロット
surf = ax.plot_surface(x, y, z, cmap='plasma', edgecolor='none')

ax.plot_surfaceを使用して、x, y, zのデータを基に3Dサーフェスプロットを作成しています。
色のマッピングにはplasmaカラーマップを使用し、エッジの色を無くしています。

7. カラーバーの追加

1
2

# カラーバーを追加
fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)

サーフェスプロットに対応するカラーバーを追加しています。
shrinkでカラーバーのサイズを調整し、aspectでアスペクト比を設定しています。

8. 軸ラベルの設定

1
2
3
4

# 軸ラベル
ax.set_xlabel('X軸')
ax.set_ylabel('Y軸')
ax.set_zlabel('Z軸')

3Dプロットの各軸にラベルを設定しています。

9. タイトルの追加

1
2

# タイトルを追加
ax.set_title('Enneper Surface')

3Dプロットにタイトルを追加しています。

10. グラフの表示

1
2

# グラフを表示
plt.show()

最終的に、plt.show()を使用してグラフを表示しています。

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

各ステップがどのように3Dプロットを作成するために機能しているかを理解することで、同様の3Dプロットを作成する際に応用することができます。

結果解説

[実行結果]

このグラフは、エンネーパーの面（Enneper Surface）を3Dで表示したものです。

エンネーパーの面は、数学や幾何学で知られる特異な曲面で、以下のような特徴を持っています：

1. 形状：

エンネーパーの面は、曲線的で非線形な形状を持ち、複雑な対称性を示します。
表面は中央部から外側に向かって曲がりくねり、波打つような形状になっています。

2. カラーマップ：

plasmaというカラーマップが使用されており、表面の高さや深さに応じて色が変化します。
これにより、3D構造の高低差が視覚的にわかりやすくなっています。

3. 軸：

• X軸：横方向の軸。
• Y軸：縦方向の軸。
• Z軸：高さ方向の軸。

4. カラーバー：

グラフの右側にカラーバーが表示されており、色と高さの対応関係を示しています。
これにより、特定の色がどの高さを表しているかがわかります。

このグラフは、エンネーパーの面の複雑な幾何学的特徴を直感的に理解できるように可視化しています。

非線形計画法（Nonlinear Programming、NLP）

非線形計画法（Nonlinear Programming、NLP）は、目的関数と制約条件が非線形となる最適化問題を扱います。

今回は、非線形最適化の典型的な例題として、以下の問題を考えます。

例題

最小化する関数と条件は以下の通りです。

$$
\text{minimize } f(x, y) = (x-1)^2 + (y-1)^2
$$

制約条件:

$$
g_1(x, y) = x^2 + y^2 - 1 \leq 0
$$
$$
g_2(x, y) = x \geq 0
$$
$$
g_3(x, y) = y \geq 0
$$

この問題をPythonのscipy.optimizeを使って解き、その結果をグラフ化します。

必要なライブラリのインストール

Google Colab環境で以下のコマンドを実行して必要なライブラリをインストールします。

1

!pip install scipy matplotlib

コード

以下がこの問題を解き、その結果をグラフ化するPythonコードです。

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

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# 最適化する関数の定義
def objective_function(x):
    return (x[0] - 1)**2 + (x[1] - 1)**2

# 制約条件の定義
def constraint1(x):
    return 1 - (x[0]**2 + x[1]**2)  # x^2 + y^2 <= 1 に対応

# 境界条件の定義
bounds = [(0, None), (0, None)]  # x >= 0, y >= 0

# 制約条件のリスト
constraints = [{'type': 'ineq', 'fun': constraint1}]

# 初期点
x0 = [0.5, 0.5]

# 最適化の実行
result = minimize(objective_function, x0, method='SLSQP', bounds=bounds, constraints=constraints)

# 結果の出力
print('最適解:', result.x)
print('最小値:', result.fun)

# 結果のグラフ化
x = np.linspace(-0.5, 1.5, 400)
y = np.linspace(-0.5, 1.5, 400)
X, Y = np.meshgrid(x, y)
Z = (X - 1)**2 + (Y - 1)**2

plt.contour(X, Y, Z, levels=30)
plt.plot(result.x[0], result.x[1], 'ro', label='Optimal Point')  # 最適解を赤い点でプロット

# 制約条件の境界線をプロット
theta = np.linspace(0, 2*np.pi, 100)
x_circle = np.cos(theta)
y_circle = np.sin(theta)
plt.plot(x_circle, y_circle, 'b-', label='$x^2 + y^2 \leq 1$')

plt.xlim(-0.5, 1.5)
plt.ylim(-0.5, 1.5)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Nonlinear Programming Optimization')
plt.grid(True)
plt.legend()
plt.fill_between(x_circle, y_circle, 1.5, where=(y_circle >= 0), alpha=0.3, color='gray')  # 可行領域
plt.show()

解説

1. 最適化する関数の定義:
   目的関数 $ ( f(x, y) = (x-1)^2 + (y-1)^2 ) $を定義します。

2. 制約条件の定義:
   非線形制約 $ ( g_1(x, y) = x^2 + y^2 - 1 \leq 0 ) $を定義します。

3. 境界条件の定義:
   $( x \geq 0 ) $および$ ( y \geq 0 ) $に対応する境界条件を定義します。

4. 制約条件のリスト constraints:
   制約条件をリスト形式で定義します。

5. 初期点の設定 x0:
   最適化アルゴリズムの初期点を設定します。

6. 最適化の実行 result:
   scipy.optimize.minimize を使用して最適化を実行します。
   methodにはSLSQP（Sequential Least Squares Programming）を使用します。

7. 結果の出力:
   得られた最適解と最小値を表示します。

8. グラフ化:

   • 目的関数の等高線を描画。
   • 制約条件 $ ( x^2 + y^2 \leq 1 ) $の境界線を描画。
   • 最適解を赤い点でプロット。
   • 制約条件を満たす領域を灰色で塗りつぶします。

以上で、非線形最適化問題の解と制約条件を含む結果をグラフで視覚的に確認できます。

実行結果

弱肉強食シミュレーション

弱肉強食の様子をアニメーションで視覚化するために、プレデター（捕食者）とプレイ（被食者）の関係をシミュレートするプログラムを作成します。

matplotlib.animationモジュールを使用して、シミュレーションの進行状況をリアルタイムで描画します。

以下は、プレデターとプレイのシミュレーションを行い、その様子をアニメーションとして可視化するPythonコードです。

（以下のコードは、Google Colab上で動作確認できます。）

必要なパッケージのインストール

まず、必要なパッケージをインストールします。

1

!pip install matplotlib numpy

Pythonコード

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
62
63
64
65
66
67
68
69
70

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

# パラメータ
num_prey = 50  # 被食者の数
num_predator = 5  # 捕食者の数
world_size = (100, 100)  # 世界のサイズ
prey_speed = 1  # 被食者の移動速度
predator_speed = 1.5  # 捕食者の移動速度

# 被食者クラス
class Prey:
    def __init__(self, x, y):
        self.x = x
        self.y = y

    def move(self):
        self.x = (self.x + np.random.uniform(-prey_speed, prey_speed)) % world_size[0]
        self.y = (self.y + np.random.uniform(-prey_speed, prey_speed)) % world_size[1]

# 捕食者クラス
class Predator:
    def __init__(self, x, y):
        self.x = x
        self.y = y

    def move(self, prey_list):
        if prey_list:
            # 最も近い被食者を探す
            prey = min(prey_list, key=lambda p: (p.x - self.x)**2 + (p.y - self.y)**2)
            dx = prey.x - self.x
            dy = prey.y - self.y
            dist = np.sqrt(dx**2 + dy**2)
            if dist > 0:
                self.x = (self.x + predator_speed * dx / dist) % world_size[0]
                self.y = (self.y + predator_speed * dy / dist) % world_size[1]

        # 捕食チェック
        for prey in prey_list.copy():
            if np.sqrt((prey.x - self.x)**2 + (prey.y - self.y)**2) < 1:
                prey_list.remove(prey)

# 初期化
preys = [Prey(np.random.uniform(0, world_size[0]), np.random.uniform(0, world_size[1])) for _ in range(num_prey)]
predators = [Predator(np.random.uniform(0, world_size[0]), np.random.uniform(0, world_size[1])) for _ in range(num_predator)]

fig, ax = plt.subplots()
scat_prey = ax.scatter([prey.x for prey in preys], [prey.y for prey in preys], c='blue')
scat_pred = ax.scatter([pred.x for pred in predators], [pred.y for pred in predators], c='red')
plt.xlim(0, world_size[0])
plt.ylim(0, world_size[1])
plt.gca().set_aspect('equal', adjustable='box')
plt.title('Weak Meat and Strong Eat Simulation')

def update(frame):
    for prey in preys:
        prey.move()
    for pred in predators:
        pred.move(preys)
        
    scat_prey.set_offsets(np.c_[[prey.x for prey in preys], [prey.y for prey in preys]])
    scat_pred.set_offsets(np.c_[[pred.x for pred in predators], [pred.y for pred in predators]])
    return scat_prey, scat_pred

ani = FuncAnimation(fig, update, frames=200, interval=100, blit=True)

# アニメーションをHTMLとして表示する
HTML(ani.to_jshtml())

コードの説明

1. インストールとインポート:
   必要なパッケージをインストールし、インポートします。

2. パラメータの設定:

   1 2 3 4 5

   num_prey = 50 # 被食者の数 num_predator = 5 # 捕食者の数 world_size = (100, 100) # 世界のサイズ prey_speed = 1 # 被食者の移動速度 predator_speed = 1.5 # 捕食者の移動速度

3. 被食者と捕食者のクラス定義:
   前述の通り、被食者と捕食者の移動ロジックを定義します。

4. 初期化:
   被食者と捕食者の位置をランダムに初期化します。

   1 2	preys = [Prey(np.random.uniform(0, world_size[0]), np.random.uniform(0, world_size[1])) for _ in range(num_prey)] predators = [Predator(np.random.uniform(0, world_size[0]), np.random.uniform(0, world_size[1])) for _ in range(num_predator)]

5. アニメーションの設定:
   アニメーション更新関数を設定し、FuncAnimationを使用してアニメーションを作成します。

   1	ani = FuncAnimation(fig, update, frames=200, interval=100, blit=True)

6. アニメーションの表示:
   HTML形式でアニメーションを表示します。

   1	HTML(ani.to_jshtml())

このコードをGoogle Colabで実行すると、弱肉強食の様子がアニメーションとして表示されます。

結果解説

[実行結果]

アニメーションの様子を具体的に理解するためには、シミュレーションがどのように進行し、どのような視覚効果があるのかを詳しく説明します。

アニメーションの概要

このシミュレーションとアニメーションは、シンプルなエコシステムをモデル化しています。
シミュレーション内で、青い点は被食者（プレイ）、赤い点は捕食者（プレデター）を表しています。
これらの点が指定されたルールに従って移動し、捕食の様子がリアルタイムで表示されます。

初期状態

• 被食者（青い点）：
  • 数量：$50$
  • ランダムな位置に配置され、移動速度は$1$です。
• 捕食者（赤い点）：
  • 数量：$5$
  • ランダムな位置に配置され、移動速度は$1.5$です。

アニメーションの進行

1. 被食者の移動：

   • 被食者は毎フレームごとにランダムな方向に少しずつ移動します。
   • 移動範囲は、世界サイズ（100x100）の中でラップアラウンド（画面端から反対側に移動）します。

2. 捕食者の移動：

   • 捕食者は最も近い被食者を追いかけます。
     具体的には、次のステップで被食者との距離が最も小さくなる方向に移動します。
   • 捕食者も同様に世界サイズ（100x100）の中でラップアラウンドします。

3. 捕食チェック：

   • 捕食者が被食者に非常に近づいた場合（距離が1未満）、その被食者は捕食されてリストから削除されます。

アニメーションのビジュアル効果

• 被食者は青い点として、捕食者は赤い点として表示されます。
• 被食者が捕食されると、その青い点は消え、アニメーション進行中に点の数が減少します。
• アニメーションは200フレーム間で進行し、各フレームごとに被食者と捕食者の位置が更新され、描画されます。

実行の様子

1. 開始時：

   • 画面にはランダムに配置された多くの青い点（被食者）と少数の赤い点（捕食者）が表示されます。

2. 中盤：

   • 青い点がランダムに動き回る一方で、赤い点が最も近い青い点を追いかける様子が確認できます。
   • 捕食が成功すると、青い点が消える様子が見えます。

3. 終了時：

   • 多くの青い点が捕食され、減少しています。
   • 赤い点が依然として動き回り、残った青い点を追いかけ続けます。

線形計画問題（Linear Programming Problem）

目的最適化の一つの例として、線形計画問題（Linear Programming Problem）を挙げることができます。

以下に、シンプルな線形計画問題をPythonで解く方法を示します。

ここではscipy.optimize.linprogを使います。

例題

問題の目的は、以下の制約条件を満たしつつ、次の目的関数を最小化することです。

目的関数:

$
f(x, y) = 3x + 4y
$

制約条件:

1. $( x + 2y \geqq 4 )$
2. $( 3x + y \geqq 3 )$
3. $( x \geqq 0 )$
4. $( y \geqq 0 )$

この問題をPythonで解く方法を以下に示します。

Pythonコード:

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

import numpy as np
from scipy.optimize import linprog

# 目的関数の係数
c = np.array([3, 4])

# 制約条件の係数行列
A = np.array([[-1, -2], [-3, -1]])

# 制約条件の右辺ベクトル
b = np.array([-4, -3])

# 変数の下限定義 (x ≥ 0, y ≥ 0)
x_bounds = (0, None)
y_bounds = (0, None)

# linprog関数を使って問題を解く
result = linprog(c, A_ub=A, b_ub=b, bounds=[x_bounds, y_bounds], method='highs')

# 結果を表示
if result.success:
    print("最適解:", result.x)
    print("目的関数値:", result.fun)
else:
    print("最適化は成功しませんでした。")

[実行結果]

解説

1. 目的関数の係数: c = [3, 4]

   • 目的関数 $ ( f(x, y) = 3x + 4y ) $の係数です。

2. 制約条件の係数行列: A = [[-1, -2], [-3, -1]]

   • 制約条件を標準形$（ (Ax \leqq b) ）$にするために、不等式の符号を反転させています。

3. 制約条件の右辺ベクトル: b = [-4, -3]

   • 同様に符号を反転させた右辺ベクトルです。

4. 変数の下限定義: x_bounds = (0, None), y_bounds = (0, None)

   • 変数$ (x, y) $が非負であることを示しています。

5. linprog関数の実行:

   • 目的関数、制約条件、変数の範囲を指定して、線形計画問題を解きます。
   • method='highs' は最先端のアルゴリズムを使って最適化を行うオプションです。

6. 結果の表示:

   • result.success が True の場合、最適解と目的関数値を表示します。

このようにして、この線形計画問題の最適化をPythonコードで解くことができます。

粒子群最適化（Particle Swarm Optimization, PSO）

粒子群最適化（Particle Swarm Optimization, PSO）は、進化的アルゴリズムの一種で、群知能を利用して最適解を見つける手法です。

具体的な問題として、ここでは以下の二次関数の最小化問題を解いてみます：

$$
f(x, y) = (x - 3)^2 + (y - 2)^2
$$

この関数の最小値を探すことが目的です。

この関数はパラボロイドであり、最小値は点$ ( (3, 2) ) $で$ 0 $となります。

以下にPythonを使ったPSOの具体的なコード例を示します。

このコードでは、PyGMOライブラリを使ってPSOを実装します。

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

import pygmo as pg
import numpy as np
import matplotlib.pyplot as plt

# 最適化したい関数の定義
def objective_function(x):
    return (x[0] - 3)**2 + (x[1] - 2)**2

# 問題の定義
class MyProblem:
    def fitness(self, x):
        return [objective_function(x)]
    def get_bounds(self):
        return ([-10, -10], [10, 10])

# 問題のインスタンスを作成
prob = pg.problem(MyProblem())

# アルゴリズムの選択 (PSO)
algo = pg.algorithm(pg.pso(gen=100))  # 100世代

# アルゴリズムのインスタンスを作成
algo.set_verbosity(1)

# 集団 (swarm) の定義
pop = pg.population(prob, size=20)

# 最適化の実行
pop = algo.evolve(pop)

# 最適化結果の取得
best_solution = pop.champion_x
best_fitness = pop.champion_f

print(f"Best solution: {best_solution}")
print(f"Best fitness: {best_fitness}")

# グラフの描画
# 評価関数の値をカラーマップでプロット
x = np.linspace(-10, 10, 400)
y = np.linspace(-10, 10, 400)
X, Y = np.meshgrid(x, y)
Z = objective_function([X, Y])

plt.contourf(X, Y, Z, levels=50, cmap='viridis')
plt.colorbar()
plt.plot(best_solution[0], best_solution[1], 'ro')  # 最適解を赤い点でプロット
plt.title("Particle Swarm Optimization Result")
plt.xlabel("x")
plt.ylabel("y")
plt.show()

ソースコード解説

1. 必要なライブラリのインポート

1
2
3

import pygmo as pg
import numpy as np
import matplotlib.pyplot as plt

PyGMO、NumPy、Matplotlibライブラリをインポートします。

2. 最適化したい関数の定義

1
2

def objective_function(x):
    return (x[0] - 3)**2 + (x[1] - 2)**2

ここで、最小化したい二次関数を定義します。

3. 最適化問題の定義

1
2
3
4
5

class MyProblem:
    def fitness(self, x):
        return [objective_function(x)]
    def get_bounds(self):
        return ([-10, -10], [10, 10])

PyGMOが理解できる形式で問題を定義します。
この例では、解の範囲を$ -10 $から$ 10 $とします。

4. 問題インスタンスの作成

1

prob = pg.problem(MyProblem())

最適化問題のインスタンスを作成します。

5. アルゴリズムの選択

1

algo = pg.algorithm(pg.pso(gen=100))  # 100世代

PSOアルゴリズムを選択し、世代数を設定します。

6. アルゴリズムインスタンスの作成

1

algo.set_verbosity(1)

アルゴリズムの詳細な実行状況を表示するように設定します。

7. 集団の定義

1

pop = pg.population(prob, size=20)

解の集団（swarm）を定義します。
ここでは $20$ 個体で初期化します。

8. 最適化の実行

1

pop = algo.evolve(pop)

最適化を実行します。

9. 最適化結果の取得

1
2
3
4
5

best_solution = pop.champion_x
best_fitness = pop.champion_f

print(f"Best solution: {best_solution}")
print(f"Best fitness: {best_fitness}")

最適な解とその評価値を取得し、表示します。

10. グラフの描画

1
2
3
4
5
6
7
8
9
10
11
12

x = np.linspace(-10, 10, 400)
y = np.linspace(-10, 10, 400)
X, Y = np.meshgrid(x, y)
Z = objective_function([X, Y])

plt.contourf(X, Y, Z, levels=50, cmap='viridis')
plt.colorbar()
plt.plot(best_solution[0], best_solution[1], 'ro')  # 最適解を赤い点でプロット
plt.title("Particle Swarm Optimization Result")
plt.xlabel("x")
plt.ylabel("y")
plt.show()

評価関数のカラーマップをプロットし、最適解を赤い点で示します。

結果解説

[実行結果]

最適化結果の解釈

• Best solution: [2.99999768, 2.00000504]

  • これはPSOアルゴリズムによって見つかった最適解です。
    関数$ ( f(x, y) = (x - 3)^2 + (y - 2)^2 ) $の最小値を見つけるために、PSOが探索した結果、最適解が$ ( x \approx 3 ) と ( y \approx 2 ) $に非常に近い値であることが示されています。

• Best fitness: [3.07979981e-11]

  • これは最適解における関数の値（評価値）です。
    最小化したい関数の最小値は$ 0 $ですが、見つけた最適解での関数値は$ ( 3.07979981 \times 10^{-11} ) $となっています。
    これは非常に小さい値であり、関数の理論的な最小値に非常に近いことを示しています。

グラフの説明

グラフは評価関数の値をカラーマップとしてプロットしたものです。
以下のような特徴があります：

• カラーマップ:

  • 色の濃淡で関数の値を表しています。

  濃い色の部分が関数の値が低い（小さい）領域を示し、薄い色の部分が関数の値が高い（大きい）領域を示します。

• 等高線:

  • 関数の等高線をプロットしています。

  等高線は、同じ関数値を持つ点を結んだ線であり、最小値に向かって収束するように見えます。

• 最適解のプロット:

  • 赤い点で最適解を示しています。
    この点が$ ( (2.99999768, 2.00000504) ) $の位置にあります。
    これにより、最適化によって見つけた最適解がグラフ上のどこにあるかが視覚的にわかります。

食事計画（線形計画問題）

CVXPYを使って線形計画問題を解きます。

ここでは、食事計画の問題を考えます。

この問題では、複数の食品から適切な組み合わせを選んで、特定の栄養要件を満たしつつ、コストを最小限に抑えることを目的とします。

問題設定

食品の種類：

1. チキン (Chicken)
2. 牛肉 (Beef)
3. 豆 (Beans)

各食品のコストと栄養素含有量は以下の通り：

• チキン: コスト 2ドル、タンパク質 20g、脂肪 5g、炭水化物 0g
• 牛肉: コスト 3ドル、タンパク質 15g、脂肪 10g、炭水化物 0g
• 豆: コスト 1ドル、タンパク質 10g、脂肪 2g、炭水化物 20g

栄養要件：

• タンパク質: 最低 30g
• 脂肪: 最低 10g、最大 20g
• 炭水化物: 最低 20g

解法

この問題をCVXPYを使って解きます。

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

import cvxpy as cp
import numpy as np

# 各食品のコストと栄養素含有量
cost = np.array([2, 3, 1])
protein = np.array([20, 15, 10])
fat = np.array([5, 10, 2])
carbs = np.array([0, 0, 20])

# 栄養要件
min_protein = 30
min_fat = 10
max_fat = 20
min_carbs = 20

# 食品の量 (変数)
x = cp.Variable(3, nonneg=True)

# コスト最小化
objective = cp.Minimize(cost @ x)

# 制約条件
constraints = [
    protein @ x >= min_protein,
    fat @ x >= min_fat,
    fat @ x <= max_fat,
    carbs @ x >= min_carbs
]

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

# 結果の表示
print(f"Status: {problem.status}")
print(f"Optimal value: {problem.value:.2f} dollars")
print(f"Optimal solution: {x.value}")
print(f"Chicken: {x.value[0]:.2f} units")
print(f"Beef: {x.value[1]:.2f} units")
print(f"Beans: {x.value[2]:.2f} units")

このコードでは、以下の手順で食事計画問題を解きます：

1. 各食品のコストと栄養素含有量を定義します。
2. 栄養要件を設定します。
3. 変数 x を定義し、各食品の量を表します（非負制約）。
4. 目的関数をコストの最小化として定義します。
5. 栄養要件に基づく制約を設定します。
6. 問題を定義し、解決します。
7. 最適なコストと食品の量を表示します。

このコードを実行すると、栄養要件を満たしつつコストを最小限に抑える最適な食品の組み合わせが得られます。

[実行結果]

Status: optimal
Optimal value: 3.72 dollars
Optimal solution: [0.64 0.48 1.  ]
Chicken: 0.64 units
Beef: 0.48 units
Beans: 1.00 units

ソースコード解説

ソースコードの詳細は以下の通りです。

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

1
2

import cvxpy as cp
import numpy as np

ここでは、最適化問題を解くためのCVXPYライブラリと、数値計算を行うためのNumPyライブラリをインポートしています。

2. 各食品のコストと栄養素含有量の定義

1
2
3
4

cost = np.array([2, 3, 1])
protein = np.array([20, 15, 10])
fat = np.array([5, 10, 2])
carbs = np.array([0, 0, 20])

ここでは、各食品（チキン、ビーフ、豆）のコストと栄養素（タンパク質、脂肪、炭水化物）の含有量を配列として定義しています。

例えば、チキンのコストは2ドルで、タンパク質が20グラム、脂肪が5グラム含まれています。

3. 栄養要件の定義

1
2
3
4

min_protein = 30
min_fat = 10
max_fat = 20
min_carbs = 20

ここでは、最低限必要なタンパク質、脂肪の最小および最大量、最低限必要な炭水化物の量を定義しています。

4. 変数の定義

1

x = cp.Variable(3, nonneg=True)

ここでは、最適化問題の変数である各食品の量を表すベクトル x を定義しています。
この変数は非負（負の値を取らない）であることが指定されています。

5. コスト最小化の目的関数

1

objective = cp.Minimize(cost @ x)

ここでは、総コストを最小化する目的関数を定義しています。
cost @ x はコストと食品の量の内積を計算し、これを最小化します。

6. 制約条件の定義

1
2
3
4
5
6

constraints = [
    protein @ x >= min_protein,
    fat @ x >= min_fat,
    fat @ x <= max_fat,
    carbs @ x >= min_carbs
]

ここでは、栄養要件に基づく制約条件を定義しています。
具体的には、次の制約を課しています：

• 総タンパク質が30グラム以上
• 総脂肪が10グラム以上20グラム以下
• 総炭水化物が20グラム以上

7. 最適化問題の定義と解決

1
2

problem = cp.Problem(objective, constraints)
problem.solve()

ここでは、目的関数と制約条件をもとに最適化問題を定義し、これを解決しています。
problem.solve() によって、最適解が計算されます。

8. 結果の表示

1
2
3
4
5
6

print(f"Status: {problem.status}")
print(f"Optimal value: {problem.value:.2f} dollars")
print(f"Optimal solution: {x.value}")
print(f"Chicken: {x.value[0]:.2f} units")
print(f"Beef: {x.value[1]:.2f} units")
print(f"Beans: {x.value[2]:.2f} units")

ここでは、最適化のステータス、最適なコスト、各食品の最適な量を表示しています。

結果解説

[実行結果]

Status: optimal
Optimal value: 3.72 dollars
Optimal solution: [0.64 0.48 1.  ]
Chicken: 0.64 units
Beef: 0.48 units
Beans: 1.00 units

• Status: optimal:
  問題が最適に解決されたことを示しています。
• Optimal value: 3.72 dollars:
  最適な食品の組み合わせのコストが3.72ドルであることを示しています。
• Optimal solution:
  各食品の最適な量を示しています。
  具体的には、チキン0.64単位、ビーフ0.48単位、豆1.00単位です。

これにより、指定された栄養要件を満たしつつ、コストを最小限に抑えるための最適な食事計画が得られます。