2021年8月26日木曜日

組み合わせ最適化1 (ナップザック問題)

 組み合わせ最適化は、基本的に説明変数が整数値をとり最適化を行います。説明変数は0または1をとるときも多いです。典型的な例でナップザック問題をやってみます。

ナップザック問題の定式化

例題として、お菓子を値段制限があるなかで、どれだけ満足度を高められるか考えます。(Pythonで始めるプログラミング入門 コロナ社 P127 プログラム10-2 を改変)

name = ['チョコ', 'ポテトチップス', 'クッキー', 'ラムネ', 'ガム']
weight = [130, 120, 80, 30, 20] #お菓子の価格
value = [18, 15, 12, 4, 2] #お菓子の満足度
capacity = 300 #お菓子の値段合計上限

プログラムは再帰呼び出しを利用して、最大の満足度になるように計算をしています。途中、リストに入っている場合はcontinueをしているので、お菓子の重複が無いようにしています。定式化でも重複が無い式になっています。

ソースコード

ソースコードの一番下に、ソルバーを使った場合のナップザック問題を入れているので、それも実行してみると同様の答えが得られます。

今回のプログラムは少ないお菓子でしたので、総当たりで計算ができますが、商品数が増えるととても計算が終わりません。その場合は、分枝限定法やメタヒューリスティックの方法で解いていきます。

復習として、下記も参考にしてください。

線形計画法1

非線形計画法1(無制約)

非線形計画法2(制約付き)

2021/09/06 ナップサック問題をdeapを使ったものを上記のソースコードに追加。遺伝アルゴリズムのフレームワークで、形式的に利用できるが、カスタマイズをするにはそれなりに内容を理解しないと使いこなせないものです。

2021年8月25日水曜日

非線形計画法2 (制約つき)

 非線形計画法での無制約のやり方は、

https://mizunolab.sist.ac.jp/2021/08/1.html

でした。ここでは、非線形計画法の中で制約つきをscipy.optimize.minimize を使ったやり方でやっていきます。この場合、数学的な知識がそれなりに必要ですが、まずはツールとして使えるようにしていきます。制約つきの中で、(i) 線形制約、(ii)非線形制約 でツールの使い方が少し変わってきます。線形計画法と同様に制約の行列表現が重要です。また非線形制約の場合は、制約条件を関数で返しています。

ソースコード

2021年8月23日月曜日

線形計画法 : LP(scipy.optimize.linprog)を利用

 線形計画法は線形関数を最適化する方法です。LPは様々なソルバーが揃っていますが、非線形計画をやることを考えてscipy.opimize.linprogを使ってみました。scipy.optimizeで非線形計画としても解けますが、微分などwarningが出てきますので、ここではlinprogを使っていきます。

scipy.opimize.linprog の仕様 : https://docs.scipy.org/doc/scipy/reference/optimize.linprog-interior-point.html

ソースコード

ポイントとしては、線形計画問題を Ax <= b の形で行列表現をして、それをソルバーに渡してあげます。等式制約の場合も同じです。x >= 0の場合もマイナスをかけて、不等式制約に当てはまるように変形します。

一つ例題3であげたもので、変数に制約が無い場合、boundsで制約をつけないとうまく計算ができませんでした。

このようなソルバーを使ったとき大規模化が問題になるので、大規模問題では検証が必要です。

2021年8月22日日曜日

非線形計画法1 (無制約)

 非線形関数を最適化することを非線形最適化と言って、利用する場面は多いです。ここでは理論的な説明は省き、ツールを使って最適化を行う方法を紹介します。

利用するソルバー:scipy.optimize

最適化手法:Nelder-Mead Simplex algorithm

pythonコード

pythonコードを確認してもらうと、いくつか例題があります。目的関数値を返す関数を作成して、その関数を呼び出してツールの書式に従って最適化を行います。パターン化されますので理解もしやすいと思います。最適化手法は今回Nelder-Meadという基本的なものを利用しています。他の最適化手法でも試して欲しいです。

20210907 追記

ローゼンブロック関数に対して、遺伝アルゴリズム(GA)のdeapライブラリを使って実施したものを追加しました。2変数での精度は出ましたが、3変数以上は難しいです。

20210908 追記

deapのベンチマークであるackley関数を追加しました。下記は参考です。

https://deap.readthedocs.io/en/master/api/benchmarks.html#deap.benchmarks.ackley

https://qiita.com/tyoshitake/items/e76f6f8e4110731606bc


2021年6月19日土曜日

Google Colaboratoryでのpythonファイル(.py)の利用について

 Google Colaboratoryはpythonを動かすのに便利な環境ですが、Jupyterの環境ですのでJupyterファイル(.ipynb)を動かすことになります。クラスやライブラリを作成してpythonファイル(.py)にしたときに、取り込みをできるようにしていきます。

(1) pythonファイルを用意する

今回は、下記のような2点を与えてその2点間距離と中点を求めるクラスのプログラムを作成してあるとします。これを自分のPCで「PointDrive.py」として保存して下さい(Colab上でやるのではなく、自分のローカルPC上でやる)。

[PointDrive.py]

(2) pythonファイル(PointDrive.py)をGoogle Driveの指定場所にアップロード

作成したPointDrive.pyをGoogle Drive上の決められた場所にアップロードします。

(3) Google ColaboratoryでJupyterファイルを開き、PointDrive.pyを利用する
Colab上でJupyterファイルを開きます。Jupyterファイルの保存場所はどこでも良いです。
(i) Google Driveをマウント

from google.colab import drive
drive.mount('/content/drive')

最初マウントを許可するための認証が入るかもしれません。

マウントができると、Jupyter上からDriveのファイルが見えます。

(ii) PointDrive.pyのパスを設定
PointDrive.pyが入っているディレクトリ(今回の例ではlibrary)を右クリックして、Copy pathをクリックします。そのパスを使って、下記のようにパスの設定をします。パスを設定後、PointDriveをインポートしてmdlとして使えるようにしています。

import sys
sys.path.append('/content/drive/MyDrive/public/python/library')
import PointDrive as mdl

(4) PointDriveクラスのインスタンスを作成して、実行する。

次のように座標を2つ与えてインスタンスを生成します。

x1 = [5, 7]
x2 = [2, 3]
p = mdl.PointDrive(x1,x2)

print('Distance = {0}, MidPoint = {1}'.format(p.get_distance(), p.get_midpoint()))

このように、クラスを利用して実行ができました。

main全体のコード


2021/08/25 追記
pythonファイルでクラスが書いてあるファイルにメインの呼び出しが書いてあるときは、もっと簡単に実行ができます。ファイルにメインの呼び出しがあるとは、ここを参考にして下さい。

このように (1) Google Driveをマウントする。(2) pythonコマンドでファイルを実行する。この2段階で実行できます。ファイルのフルパスを得る必要がありますが、それができれば簡単です。

その他の注意
Google Colab上で入力や結果でcsvを使う場合、これもGoogle Drive上でフルパスで指定しなければならず、意外に面倒です。

2021年6月8日火曜日

動的計画法(2) M/M/1の計算

 数式で漸化式で表現されれば、再帰関数を利用して動的計画法で計算できるということでした。今回は待ち行列のM/M/1の定常分布を動的計画法で求めていきます。

(1) 待ち行列 M/M/1の表現

式の表現は、[1]やhttps://www.win.tue.nl/~iadan/que/h4.pdfを見て欲しいのですが、リンクの式(2), (3), (4)がM/M/1の定常分布を出すための平衡方程式です。これをDPで解いていくのですが、その前にこの平衡方程式を解けば、定常分布の理論解が得られます。理論解は式(6)です。要するにここでは、理論解で得られた値とDPで得られた結果を比較して精度を確認していきます。

(2) 理論解の導出

この条件で、系内人数が10人となる確率は0.014078378677368164となります。

(2) DPでの算出
次に漸化式を利用して、再帰関数を作り、DPで同じ確率を算出します。漸化式を整理すると、3項間漸化式になり、pi(0)とpi(1)を初期値で設定します。
これをプログラムにします。
それなりの精度で出てきます。0.014078378677368164です。

このように漸化式であらわされているものは、再帰関数を用いて動的計画法として解いていくと便利です。

参考

[1] Bolch, G., Greiner, S., De Meer, H., & Trivedi, K. S. (2006). Queueing networks and Markov chains: modeling and performance evaluation with computer science applications. John Wiley & Sons.   (P105 Birth Death(λ_i、μ_i) P214でMM1(λ、μ)を与えている)


2021年6月6日日曜日

動的計画法(1)

 プログラムで再帰関数を利用することで、プログラムを簡略して書くことができます。特に漸化式で表されている数式をプログラミングする場合は、非常に簡単にプログラムができます。ここでは再帰関数が動的計画法で利用する場面を紹介します。プログラムの環境はGoogle Colabです。

(1) フィボナッチ数列をプログラムする ([1]参照)

フィボナッチ数列は前の2項(n-1とn-2)の和を現在の項(n)にするという数列です。これをプログラミングすると、下記になります。

これはフィボナッチ数列の第10項までを実行したものですので簡単に終わります。しかし、40項まで実行してみると、1分程度かかってしまいます。下記は3回の実行結果です。
40項 = 102334155 calclation_time:62.73430848121643[sec]
40項 = 102334155 calclation_time:62.58229374885559[sec]
40項 = 102334155 calclation_time:63.6257905960083[sec]

再帰計算を行うとき、同じ計算過程を何度も行う場合があり、計算上無駄が多くなっています。

(2) メモ化を利用する ([2参照])
メモ化は既に計算した値を確保しておき、必要になったときに利用していきます。既に計算されたものは値があるので、計算が驚くほど早くなります。下記は100000件を動かした結果です。
calclation_time:1.049574851989746[sec]
100000件でも1秒程度で計算できています。1,000,000件にするとメモリエラーになりました。

(3) 課題
メモ化でかなりの高速化ができるということがわかったが、メモリも消費する。MPIで並列計算を行い、メモリ分散ができるようにしてみたい。

[参考]
[1] 大和田 勇人, 金盛 克俊, Pythonで始めるプログラミング入門, コロナ社, (2015)
[2] Python言語によるプログラミングイントロダクション第2版 (著者省略), (2017)

2021年2月5日金曜日

matplotlib日本語化対応

matplotlibで日本語データをグラフ表示すると文字化けしてしまいます。

そこで日本語化への対応方法を示します。




フォントファイルの取得
まず、フォントファイルを取得していきます。
今回は、記事作成時点で最新のIPAexフォント(Ver.004.01)のIPAexゴシックを利用してきます。
以下のリンクにアクセスしてフォントファイルIPAexゴシック(Ver.004.01)をダウンロードします。
https://moji.or.jp/ipafont/ipaex00401/


ダウンロードしたら展開します。



フォントファイルを配置
次にmatplotlibでこのフォントファイルを利用できるように、フォントファイルを配置していきます。

matplotlibの設定ファイルの場所は、以下のコードにて確認することができます。
import matplotlib as mpl
print(mpl.matplotlib_fname())
私の場合は、「/opt/anaconda3/lib/python3.8/site-packages/matplotlib/mpl-data/matplotlibrc」でした。

フォントファイルは、設定ファイルの1階層上のfontsフォルダ内のtff内に入れていきます。
私の場合は、「/opt/anaconda3/lib/python3.8/site-packages/matplotlib/mpl-data/fonts/ttf」にipaexg.tffを入れます。


キャッシュファイルの削除
使用しているフォントのキャッシュが残っている場合があります。
そこでキャッシュファイルを削除していきます。
windowsの場合は、「C:\Users\<ユーザ名>\.matplotlib」
Macの場合は、「/Users/<ユーザ名>/.matplotlib」にあるjsonファイルを削除します。

そして、一度現在使用しているVSCodeやJupterNotebookなどは落として再起動します。


matplotlibの設定
フォントファイルがフォルダに入り、キャッシュも削除されました。
このままではまだ日本語化されません。matplotlibを実行するときに、使用するフォントファイルを指定してあげる必要があります。
ソースコード内に
plt.rcParams['font.family'] = 'IPAexGothic'
を記載することで、IPAexゴシックを利用するようになります。

これにて終了です。


参考