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$)を満たすことができます。

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

Privacy Policy

Privacy Policy

This privacy policy applies to the Piano BGM app (hereby referred to as “Application”) for mobile devices that was created by (hereby referred to as “Service Provider”) as a Free service. This service is intended for use “AS IS”.

Information Collection and Use

The Application collects information when you download and use it. This information may include information such as

  • Your device’s Internet Protocol address (e.g. IP address)
  • The pages of the Application that you visit, the time and date of your visit, the time spent on those pages
  • The time spent on the Application
  • The operating system you use on your mobile device

The Application does not gather precise information about the location of your mobile device.

The Application collects your device’s location, which helps the Service Provider determine your approximate geographical location and make use of in below ways:

  • Geolocation Services: The Service Provider utilizes location data to provide features such as personalized content, relevant recommendations, and location-based services.
  • Analytics and Improvements: Aggregated and anonymized location data helps the Service Provider to analyze user behavior, identify trends, and improve the overall performance and functionality of the Application.
  • Third-Party Services: Periodically, the Service Provider may transmit anonymized location data to external services. These services assist them in enhancing the Application and optimizing their offerings.

The Service Provider may use the information you provided to contact you from time to time to provide you with important information, required notices and marketing promotions.

For a better experience, while using the Application, the Service Provider may require you to provide us with certain personally identifiable information. The information that the Service Provider request will be retained by them and used as described in this privacy policy.

Third Party Access

Only aggregated, anonymized data is periodically transmitted to external services to aid the Service Provider in improving the Application and their service. The Service Provider may share your information with third parties in the ways that are described in this privacy statement.

Please note that the Application utilizes third-party services that have their own Privacy Policy about handling data. Below are the links to the Privacy Policy of the third-party service providers used by the Application:

The Service Provider may disclose User Provided and Automatically Collected Information:

  • as required by law, such as to comply with a subpoena, or similar legal process;
  • when they believe in good faith that disclosure is necessary to protect their rights, protect your safety or the safety of others, investigate fraud, or respond to a government request;
  • with their trusted services providers who work on their behalf, do not have an independent use of the information we disclose to them, and have agreed to adhere to the rules set forth in this privacy statement.

Opt-Out Rights

You can stop all collection of information by the Application easily by uninstalling it. You may use the standard uninstall processes as may be available as part of your mobile device or via the mobile application marketplace or network.

Data Retention Policy

The Service Provider will retain User Provided data for as long as you use the Application and for a reasonable time thereafter. If you’d like them to delete User Provided Data that you have provided via the Application, please contact them at angularity@msn.com and they will respond in a reasonable time.

Children

The Service Provider does not use the Application to knowingly solicit data from or market to children under the age of 13.

The Application does not address anyone under the age of 13. The Service Provider does not knowingly collect personally identifiable information from children under 13 years of age. In the case the Service Provider discover that a child under 13 has provided personal information, the Service Provider will immediately delete this from their servers. If you are a parent or guardian and you are aware that your child has provided us with personal information, please contact the Service Provider (angularity@msn.com) so that they will be able to take the necessary actions.

Security

The Service Provider is concerned about safeguarding the confidentiality of your information. The Service Provider provides physical, electronic, and procedural safeguards to protect information the Service Provider processes and maintains.

Changes

This Privacy Policy may be updated from time to time for any reason. The Service Provider will notify you of any changes to the Privacy Policy by updating this page with the new Privacy Policy. You are advised to consult this Privacy Policy regularly for any changes, as continued use is deemed approval of all changes.

This privacy policy is effective as of 2024-06-18

Your Consent

By using the Application, you are consenting to the processing of your information as set forth in this Privacy Policy now and as amended by us.

Contact Us

If you have any questions regarding privacy while using the Application, or have questions about the practices, please contact the Service Provider via email at angularity@msn.com.


Rastrigin関数

Rastrigin関数

Rastrigin関数の最小化を考えます。

Rastrigin関数は、多次元の最適化問題でよく使われる非線形かつ非凸な関数で、複数の局所最小値を持つため、グローバル最適化の良い例題です。

Rastrigin関数の定義は以下の通りです:

$$
f(\mathbf{x}) = An + \sum_{i=1}^{n} \left[ x_i^2 - A \cos(2 \pi x_i) \right]
$$

ここで、通常$ ( A = 10 ) $で、$ ( n ) $は次元数です。

この関数のグローバル最小値は$ (\mathbf{x} = \mathbf{0}) $であり、そのときの関数値は$0$です。

プログラム例

以下は、PythonRastrigin関数を最適化するためのコードです。

scipy.optimizeモジュールのdual_annealing関数を使用してグローバル最適化を行います。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import numpy as np
from scipy.optimize import dual_annealing

# Rastrigin関数の定義
def rastrigin(x):
A = 10
return A * len(x) + sum([(xi**2 - A * np.cos(2 * np.pi * xi)) for xi in x])

# 最適化の設定
bounds = [(-5.12, 5.12) for _ in range(2)] # 2次元のRastrigin関数

# dual_annealingによる最適化
result = dual_annealing(rastrigin, bounds)

# 結果の表示
print("最適化結果:")
print("最小値:", result.fun)
print("最適解:", result.x)
print("反復回数:", result.nit)

ソースコード解説

  1. 関数の定義

    • Rastrigin関数を定義します。
      この関数は入力ベクトル x に対して評価されます。
  2. 最適化の設定

    • 最適化する範囲(ここでは、各次元が ([-5.12, 5.12]) の範囲にあります)を設定します。
  3. dual_annealingによる最適化

    • SciPyのdual_annealing関数を使用して、Rastrigin関数の最小値を求めます。
  4. 結果の表示

    • 最適化結果として、最小値最適解、および反復回数を表示します。

このコードを実行することで、Rastrigin関数のグローバル最小値を求めることができます。

Rastrigin関数の性質上、最適化アルゴリズムが局所最小値に陥らないようにするための良いテストケースとなります。

結果解説

[実行結果]

最適化結果:
最小値: 1.8474111129762605e-13
最適解: [-3.00458545e-08 -5.08182682e-09]
反復回数: 1000
  1. 最小値: 1.8474111129762605e-13

    • これはRastrigin関数の最小値で、最適化アルゴリズムが見つけた関数の評価値です。
      この値は非常に小さいことから、ほぼゼロに近いことを示しています。
      Rastrigin関数の理論的なグローバル最小値は$0$なので、この結果は非常に良い最適化が行われたことを示しています。
  2. 最適解: [-3.00458545e-08 -5.08182682e-09]

    • これは関数が最小値を取るときの入力ベクトル(解)です。
      ここで示されている値は、非常に小さい数値(ほぼゼロ)です。
      理論的な最適解はすべての変数がゼロのベクトル$ (\mathbf{0}) $なので、この結果は理論的な最適解に非常に近いことを示しています。
    • 具体的には、$ (-3.00458545 \times 10^{-8}) $と$ (-5.08182682 \times 10^{-9}) $はどちらもゼロに非常に近い値です。
  3. 反復回数: 1000

    • これは最適化アルゴリズムが終了するまでに行われた反復(イテレーション)の回数です。
      dual_annealingアルゴリズムが設定された最大の反復回数に達したことを示しています。
      最適化が成功し、非常に小さい最小値を見つけた一方で、アルゴリズムは最大反復回数まで探索を続けました。

最適化アルゴリズムと結果の評価

  • dual_annealing アルゴリズム:

    • これは焼きなまし法(Simulated Annealing)に基づくグローバル最適化アルゴリズムです。
      多くの局所最小値を持つ関数(非凸関数)の最適化に適しています。
    • このアルゴリズムは、温度を徐々に下げながら探索範囲を狭めていくことで、局所最小値に陥ることを防ぎつつ、グローバル最小値を探索します。
  • 結果の評価:

    • 最小値がほぼゼロであり、最適解が理論的な最適解(ゼロベクトル)に非常に近いことから、アルゴリズムは非常に効果的に動作したと評価できます。
    • 最大反復回数に達したものの、最適化の目的は十分に達成されたと考えられます。
      探索プロセスをさらに続ける必要はないでしょう。

まとめ

この結果は、Rastrigin関数の最小化においてdual_annealingアルゴリズムが有効に機能したことを示しています。

最適化アルゴリズムは、非常に小さな関数値と理論的な最適解に近い解を見つけることができました。

医療リソース最適化

医療リソース最適化

医療リソース最適化の問題を考え、Pythonで解いてみましょう。

ここでは、病院内のベッドの配分を最適化する問題を例にします。

具体的には、限られた数のベッドを複数の患者グループにどのように配分すれば、全体の患者満足度を最大化できるかを考えます。

問題設定

  • 病院には合計$ (N) $台のベッドがある。
  • 患者グループは$ (G) $グループあり、それぞれのグループには異なる数の患者がいる。
  • 各患者グループ$ (i) $の患者数は$ (P_i) $であり、満足度は$ (U_i) $で与えられる。
  • ベッドの割り当てを最適化して、患者満足度の総和を最大化する。

数理モデル

この問題は整数線形計画問題としてモデル化できます。

変数

  • $(x_i)$: 患者グループ$ (i) $に割り当てるベッドの数(整数)

目的関数

最大化する目的関数は全体の満足度の総和です。
$$
\text{maximize} \quad \sum_{i=1}^{G} U_i \cdot x_i
$$

制約条件

  • 各グループに割り当てるベッドの総数は病院のベッド数を超えてはいけません。
    $$
    \sum_{i=1}^{G} x_i \leq N
    $$
  • 各グループに割り当てるベッドの数は、そのグループの患者数を超えてはいけません。
    $$
    0 \leq x_i \leq P_i \quad \forall i
    $$

Python実装

Pythonでこの問題を解くために、数理最適化ライブラリである PuLP を使用します。

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

1
pip install pulp

次に、以下の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
import pulp

# パラメータの設定
N = 50 # 病院のベッド数
P = [10, 20, 30, 40] # 各患者グループの患者数
U = [5, 4, 3, 2] # 各患者グループの満足度
G = len(P) # 患者グループの数

# 問題の定義
problem = pulp.LpProblem("Bed_Allocation_Optimization", pulp.LpMaximize)

# 変数の定義
x = [pulp.LpVariable(f"x{i}", lowBound=0, upBound=P[i], cat='Integer') for i in range(G)]

# 目的関数の設定
problem += pulp.lpSum([U[i] * x[i] for i in range(G)]), "Total_Satisfaction"

# 制約条件の設定
problem += pulp.lpSum(x) <= N, "Total_Beds_Constraint"
for i in range(G):
problem += x[i] <= P[i], f"Patient_Group_{i}_Constraint"

# 問題の解決
problem.solve()

# 結果の表示
print(f"Status: {pulp.LpStatus[problem.status]}")
for i in range(G):
print(f"Group {i + 1}: Allocate {x[i].varValue} beds")

print(f"Total Satisfaction: {pulp.value(problem.objective)}")

# 結果の可視化
import matplotlib.pyplot as plt

allocated_beds = [x[i].varValue for i in range(G)]
group_labels = [f"Group {i + 1}" for i in range(G)]

plt.bar(group_labels, allocated_beds)
plt.xlabel("Patient Groups")
plt.ylabel("Allocated Beds")
plt.title("Optimal Bed Allocation")
plt.show()

結果の説明

このコードは、病院のベッド配分問題を整数線形計画問題としてモデル化し、最適なベッド配分を求めます。

最後に、各患者グループに割り当てられたベッドの数と、合計の患者満足度を表示し、グラフで可視化します。

[実行結果]

ソースコード解説

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

1
2
import pulp
import matplotlib.pyplot as plt
  • pulp は線形計画法(LP)を解くためのライブラリです。
  • matplotlib.pyplot はグラフの描画に使用されます。

2. パラメータの設定

1
2
3
4
N = 50  # 病院のベッド数
P = [10, 20, 30, 40] # 各患者グループの患者数
U = [5, 4, 3, 2] # 各患者グループの満足度
G = len(P) # 患者グループの数
  • N は病院の総ベッド数です。
  • P は各患者グループの患者数を表します。
  • U は各患者グループの満足度を表します。
  • G は患者グループの数です。

3. 問題の定義

1
problem = pulp.LpProblem("Bed_Allocation_Optimization", pulp.LpMaximize)
  • ここでは、問題を「Bed_Allocation_Optimization」という名前で定義し、目的が最大化(満足度の合計)であることを指定しています。

4. 変数の定義

1
x = [pulp.LpVariable(f"x{i}", lowBound=0, upBound=P[i], cat='Integer') for i in range(G)]
  • 各患者グループに割り当てるベッド数を表す変数 x[i] を定義します。
  • lowBound=0 は変数の下限を0に設定し、upBound=P[i] は各グループの患者数を上限とします。
  • cat='Integer' は変数が整数であることを指定しています。

5. 目的関数の設定

1
problem += pulp.lpSum([U[i] * x[i] for i in range(G)]), "Total_Satisfaction"
  • 各グループの満足度 U[i] に割り当てるベッド数 x[i] を掛け合わせた総和を目的関数として定義しています。
  • pulp.lpSum はリストの総和を計算します。

6. 制約条件の設定

1
2
3
problem += pulp.lpSum(x) <= N, "Total_Beds_Constraint"
for i in range(G):
problem += x[i] <= P[i], f"Patient_Group_{i}_Constraint"
  • pulp.lpSum(x) <= N は、全体のベッド数が N 以下であるという制約を設定します。
  • ループ内では、各グループのベッド割り当てがそのグループの患者数を超えないように制約を追加しています。

7. 問題の解決

1
problem.solve()
  • problem.solve() は定義された線形計画問題を解きます。

8. 結果の表示

1
2
3
4
print(f"Status: {pulp.LpStatus[problem.status]}")
for i in range(G):
print(f"Group {i + 1}: Allocate {x[i].varValue} beds")
print(f"Total Satisfaction: {pulp.value(problem.objective)}")
  • pulp.LpStatus[problem.status] は解のステータス(最適、不可解など)を表示します。
  • ループ内では各グループに割り当てられたベッド数を表示します。
  • pulp.value(problem.objective) は最適な満足度の合計を表示します。

9. 結果の可視化

1
2
3
4
5
6
7
8
allocated_beds = [x[i].varValue for i in range(G)]
group_labels = [f"Group {i + 1}" for i in range(G)]

plt.bar(group_labels, allocated_beds)
plt.xlabel("Patient Groups")
plt.ylabel("Allocated Beds")
plt.title("Optimal Bed Allocation")
plt.show()
  • allocated_beds は各グループに割り当てられたベッド数をリストにします。
  • group_labels は各グループのラベルを作成します。
  • plt.bar は棒グラフを作成します。
  • plt.xlabel, plt.ylabel, plt.title はグラフのラベルとタイトルを設定します。
  • plt.show() はグラフを表示します。

このスクリプトは、病院の限られたベッド数を最も効率的に配分することで患者満足度を最大化するための線形計画問題を解き、その結果を視覚化します。

結果解説

[実行結果]

最適なベッド配分

問題を解いた結果、以下のようなベッド配分が最適であることが示されました:

  • Group 1: $10$台のベッドを割り当てる
  • Group 2: $20$台のベッドを割り当てる
  • Group 3: $20$台のベッドを割り当てる
  • Group 4: ベッドを割り当てない

合計の満足度

各グループに割り当てたベッドの数に応じて、全体の満足度の合計は$190$となります。

詳細な解析

  • Group 1:
    このグループには$10$人の患者がいて、満足度は$5$です。
    ベッドを10台割り当てると、満足度は$ (10 \times 5 = 50) $になります。
  • Group 2:
    このグループには$20$人の患者がいて、満足度は$4$です。
    ベッドを$20$台割り当てると、満足度は$ (20 \times 4 = 80) $になります。
  • Group 3:
    このグループには$30$人の患者がいて、満足度は$3$です。
    ベッドを$20$台割り当てると、満足度は$ (20 \times 3 = 60) $になります。
  • Group 4:
    このグループには$40$人の患者がいて、満足度は$2$です。
    しかし、このグループにはベッドが割り当てられていないため、満足度は$0$です。

これらを合計すると、全体の満足度は$ (50 + 80 + 60 + 0 = 190) $になります。