二次計画問題 (Quadratic Programming, QP)

二次計画問題 (Quadratic Programming, QP)

凸最適化問題の具体例として、二次計画問題 (Quadratic Programming, QP) を挙げます。
この問題は、凸二次関数を最小化するもので、以下の形式を持ちます:

$$
\text{minimize} \quad \frac{1}{2}x^T Q x + c^T x
$$
$$
\text{subject to} \quad Ax \leq b
$$

ここで、$( Q ) $は対称かつ正定値行列、$( c ) $はベクトル、$( A ) $は行列、$( b ) $はベクトルです。

具体例:
$$ Q = \begin{bmatrix} 2 & 0 \ 0 & 2 \end{bmatrix}, \quad c = \begin{bmatrix} -4 \ -6 \end{bmatrix}, \quad A = \begin{bmatrix} 1 & 1 \ -1 & 2 \ 2 & 1 \end{bmatrix}, \quad b = \begin{bmatrix} 3 \ 2 \ 4 \end{bmatrix} $$

この例では、二次関数$ (\frac{1}{2}x^T Q x + c^T x) $を最小化し、線形制約$ (Ax \leq b) $を満たす$ (x) $を見つけます。

以下は、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
import cvxpy as cp
import numpy as np

# 定義されたパラメータ
Q = np.array([[2, 0], [0, 2]])
c = np.array([-4, -6])
A = np.array([[1, 1], [-1, 2], [2, 1]])
b = np.array([3, 2, 4])

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

# 目的関数の定義
objective = cp.Minimize(0.5 * cp.quad_form(x, Q) + c.T @ x)

# 制約の定義
constraints = [A @ x <= b]

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

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

# 結果の表示
print(f"Optimal value: {result}")
print(f"Optimal solution: x = {x.value}")

このコードでは、次の手順を踏みます:

  1. cvxpyライブラリをインポートします。
  2. パラメータ$ (Q)$、$(c)$、$(A)$、および$ (b) $を定義します。
  3. 最適化変数$ (x) $を定義します。
  4. 目的関数を定義します。
  5. 制約を定義します。
  6. 問題を定義し、解きます。
  7. 最適解と最適値を表示します。

実行すると、最適な解$ (x) $とそのときの最小値を得ることができます。

[実行結果]

Optimal value: -10.399999999999999
Optimal solution: x = [1.2 1.6]

ソースコード解説

コードの各部分を詳しく説明します。

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

1
2
import cvxpy as cp
import numpy as np

最初に、最適化ライブラリであるCVXPYと数値計算ライブラリであるNumPyをインポートしています。
CVXPY凸最適化問題を簡潔に定義し解くためのライブラリです。

2. パラメータの定義

1
2
3
4
5
# 定義されたパラメータ
Q = np.array([[2, 0], [0, 2]])
c = np.array([-4, -6])
A = np.array([[1, 1], [-1, 2], [2, 1]])
b = np.array([3, 2, 4])

ここでは、最適化問題に必要なパラメータを定義しています。

  • Q: 二次形式の係数行列(正定値行列)
  • c: 線形項の係数ベクトル
  • A: 制約条件の係数行列
  • b: 制約条件の右辺ベクトル

3. 変数の定義

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

次に、最適化変数 x を定義しています。
この変数は2次元のベクトルです。

4. 目的関数の定義

1
2
# 目的関数の定義
objective = cp.Minimize(0.5 * cp.quad_form(x, Q) + c.T @ x)

目的関数を定義しています。
この問題の目的は、以下の形式の関数を最小化することです:

$$
\frac{1}{2} x^T Q x + c^T x
$$

cp.quad_form(x, Q) は$ (x^T Q x) $を計算し、c.T @ x は$ (c^T x) $を計算します。
これらを足し合わせて目的関数を作成し、cp.Minimize を使って最小化します。

5. 制約の定義

1
2
# 制約の定義
constraints = [A @ x <= b]

制約条件を定義しています。
ここでは、以下の形式の線形不等式制約を設定しています:

$$
A x \leq b
$$

6. 問題の定義

1
2
# 問題の定義
prob = cp.Problem(objective, constraints)

最適化問題全体を定義しています。
cp.Problem 関数は、目的関数制約条件を受け取り、最適化問題を作成します。

7. 問題を解く

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

定義した最適化問題を解きます。
prob.solve() メソッドは、問題を解いて最適値を返します。

8. 結果の表示

1
2
3
# 結果の表示
print(f"Optimal value: {result}")
print(f"Optimal solution: x = {x.value}")

最適化の結果を表示します。

  • result最適値(目的関数の最小値)です。
  • x.value最適解(変数 x の値)です。

全体の流れ

  1. ライブラリをインポートし、パラメータを定義。
  2. 最適化変数を定義。
  3. 目的関数を定義。
  4. 制約条件を定義。
  5. 目的関数と制約条件を組み合わせて最適化問題を定義。
  6. 問題を解く。
  7. 結果を表示する。

このコードは、二次計画問題を解くための一般的な方法を示しており、CVXPYライブラリを使うことで問題の定義と解決を簡潔に行えることを示しています。

結果解説

Optimal value: -10.399999999999999

最適値(目的関数の最小値)は$ (-10.4) $です。
これは、与えられた二次計画問題において、制約条件 $(Ax \leq b) $を満たす$ (x) $に対して目的関数が最小となる値です。

Optimal solution: x = [1.2, 1.6]

最適解$ (x = [1.2, 1.6]) $です。
これは、目的関数$ (\frac{1}{2}x^T Q x + c^T x) $が最小値をとる$ (x) $の値です。
この$ (x) $に対して、制約条件$ (Ax \leq b) $がすべて満たされます。

計算の詳細

目的関数の計算

与えられた$ (x = [1.2, 1.6]) $に対して目的関数を計算します。

$$
\frac{1}{2}x^T Q x + c^T x
$$

まず、二次項を計算します:
$$
x^T Q x = \begin{bmatrix} 1.2 & 1.6 \end{bmatrix} \begin{bmatrix} 2 & 0 \ 0 & 2 \end{bmatrix} \begin{bmatrix} 1.2 \ 1.6 \end{bmatrix} = 1.2 \times 2 \times 1.2 + 1.6 \times 2 \times 1.6 = 2.88 + 5.12 = 8.0
$$

次に、線形項を計算します:
$$
c^T x = \begin{bmatrix} -4 & -6 \end{bmatrix} \begin{bmatrix} 1.2 \ 1.6 \end{bmatrix} = -4 \times 1.2 + -6 \times 1.6 = -4.8 - 9.6 = -14.4
$$

したがって、目的関数の値は:
$$
\frac{1}{2} \times 8.0 + (-14.4) = 4.0 - 14.4 = -10.4
$$

これは、表示された最適値と一致します。

制約条件の確認

次に、与えられた最適解 $(x = [1.2, 1.6]) $が制約を満たすことを確認します:

$$
Ax = \begin{bmatrix} 1 & 1 \ -1 & 2 \ 2 & 1 \end{bmatrix} \begin{bmatrix} 1.2 \ 1.6 \end{bmatrix} = \begin{bmatrix} 1 \times 1.2 + 1 \times 1.6 \ -1 \times 1.2 + 2 \times 1.6 \ 2 \times 1.2 + 1 \times 1.6 \end{bmatrix} = \begin{bmatrix} 2.8 \ 2.0 \ 4.0 \end{bmatrix}
$$

これに対して、制約条件 $(b = [3, 2, 4]) $と比較します:
$$
\begin{bmatrix} 2.8 \ 2.0 \ 4.0 \end{bmatrix} \leq \begin{bmatrix} 3 \ 2 \ 4 \end{bmatrix}
$$

すべての不等式が満たされていることが確認できます。

結論

この問題に対する最適解は$ (x = [1.2, 1.6]) $であり、最適値は$ (-10.4) $です。
与えられた制約条件を満たしつつ、目的関数の値を最小化する解が求められました。

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

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

問題:

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)

二次計画法(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. 行列およびベクトルの定義:

    • PqGhAb それぞれ問題の行列およびベクトルを定義します。
  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)

バーンズリーのシダ(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()

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

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

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

輸送問題 (Transportation Problem)

輸送問題 (Transportation Problem)

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

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

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

問題設定:

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

この輸送問題を 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)

エンネーパーの面(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)

ここでは、uvの範囲を-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)

非線形計画法(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)

目的最適化の一つの例として、線形計画問題(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.successTrue の場合、最適解目的関数値を表示します。

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

粒子群最適化(Particle Swarm Optimization, PSO)

粒子群最適化(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

PyGMONumPyMatplotlibライブラリをインポートします。

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) ) $の位置にあります。
      これにより、最適化によって見つけた最適解がグラフ上のどこにあるかが視覚的にわかります。