PYTHONによる多目的最適化プログラム例: 新道ブログ

PYTHONはいろいろなところで使われているが、なかなか概念が理解できない。関数、オブジェクト、クラスなど抽象的な対象物がプログラム言語化されるので、用語の相対関係、引用、参照関係がプログラムリストから読み取れないのである。しかし、データ処理機能のパッケージ化が世界中で行われているので、それを利用すれば簡単に高度なデータ処理ができる。
　下記は、多目的最適化の例（目的関数２個，制約関数１個，設計変数３個）を３次元グラフ表示できるようにしたものである。
　ただし、下記コードでは、アサブロの改行ルールの影響を受け、PYTHON文の各行の頭のブランクが削除される。このため、本コードをPYTHONで流す前に、Def行およびFor行に続く行（#行又はブランク行が現れる範囲）の頭に半角ブランクを4桁入力する必要がある。（2022.5．6改良版アップ）
（なお、PYTHON及びJUPYTER　NOTEBOOKのインストールはネット情報を用いている。）

#0506-一般化PYTHONN14-27ケースサーベイエクセル回帰3D棒グラフColorZ.txt
#https://sidestory.pandanote.info/6890bis.html
#日本語グラフ表示フォント利用版
from platypus import NSGAII, Problem, Real, nondominated
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
#行数無制限指定forＮｏｔｅｂｏｏｋ
from IPython.core.display import display, HTML

display(HTML("<style>div.output_scroll { height: unset; }</style>"))
np.random.seed(0)

#%matplotlib inline
#3次元グラフパレート解＋回転可能
#%matplotlib tk
%matplotlib notebook

# 目的関数 by 0409ケースサーベイエクセル回帰
#x1:設計変数１
#x2:設計変数２
#x3:設計変数３
def objective(vars):
x1 = vars[0]
x2 = vars[1]
x3 = vars[2]
f1 = 0.001284257*x1 -0.004223327*x2 +9.943736814*x3 -9.502343321 #設計変数１
f2 = 0.073949761*x1 -0.006730252*x2 +5.427614607*x3 -7.110353214 #設計変数２
v1 =0.000652886*x1 +0.000100475*x2 +0.340144444*x3 +0.618259667 #制約関数
return [f1, f2], [v1]

# 最適化計算を実行する関数
def optimization(n_var, n_obj, n_con, vars, n_run):
# 設計変数と目的関数の数を設定
problem = Problem(n_var, n_obj, n_con)
# 最適化計算の条件
problem.directions[0] = Problem.MINIMIZE
problem.directions[1] = Problem.MINIMIZE
problem.constraints[:] = ">=1.002"
# 設計変数と目的関数を設定
problem.types[:] = vars
problem.function = objective
# 最適化アルゴリズムを設定して計算を実行する
algorithm = NSGAII(problem)
algorithm.run(n_run)

return algorithm

# 変数を設定する
var1 = Real(0., 30.)
var2 = Real(0., 60.)
var3 = Real(1.00, 1.20)
vars = [var1, var2, var3]

# 最適化計算を実行する
algorithm = optimization(n_var=3, n_obj=2, n_con=1, vars=vars, n_run=3000)

# ここからグラフ描画
# フォントの種類とサイズを設定する。日本語用にMS Gothicを設定した。
plt.rcParams['font.size'] = 14
plt.rcParams['font.family'] = 'MS Gothic'

# グラフの入れ物を用意する。
fig1 = plt.figure(figsize=(2,2))
ax1 = fig1.add_subplot(222,projection='3d')
ax2 = fig1.add_subplot(223,projection='3d')
ax3 = fig1.add_subplot(224,projection='3d')
ax4 = fig1.add_subplot(221)
# 軸のラベルを設定する。
ax1.set_xlabel('目的関数１')
ax1.set_ylabel('目的関数２')
ax1.set_zlabel('設計変数１)
ax2.set_xlabel('目的関数１')
ax2.set_ylabel('目的関数２')
ax2.set_zlabel('設計変数２')
ax3.set_xlabel('目的関数１')
ax3.set_ylabel('目的関数２')
ax3.set_zlabel('制約関数１')
ax4.set_ylabel('目的関数１')
ax4.set_xlabel('目的関数２')

# 軸の最大値・最小値の設定
ax1.set_xlim(1.5, 0.5)
ax1.set_ylim(-1.5, 0.5)
ax1.set_zlim(0.0, 30.0)
ax2.set_xlim(1.5, 0.5)
ax2.set_ylim(-1.5, 0.5)
ax2.set_zlim(0.0, 60.0)
ax3.set_xlim(1.5, 0.5)
ax3.set_ylim(-1.5, 0.5)
ax3.set_zlim(1.0, 1.1)
plt.xlim(-1.5, 0.5)
plt.ylim(0.5, 1.5)


#
# 非劣解を抽出する
n_dominated = nondominated(algorithm.result)
# パレート解をプリントする
# パレート解をプロットする。

var1= []
obj1 = []
obj2 = []
for solution in n_dominated:
obj1.append(solution.objectives[0])
obj2.append(solution.objectives[1])
var1.append(solution.variables[0])
# print('Variables=', solution.variables[:])
# print('Objectives=', solution.objectives[:])
# print('Constraints=', solution.constraints[:])
#
ax1.bar3d(obj1, obj2, 0.0, 0.01, 0.01, var1,color='red', label='Pareto front1' )
ax1.legend([obj1,obj2],["A","B"])
# パレート解をプロットする。

var2= []
obj1 = []
obj2 = []
for solution in n_dominated:
obj1.append(solution.objectives[0])
obj2.append(solution.objectives[1])
var2.append(solution.variables[1])
#
ax2.bar3d(obj1, obj2, 0.0, 0.01, 0.01, var2, color='red', label='Pareto front2' )
ax2.legend([obj1,obj2],["C","D"])

# パレート解をプロットする。

var3= []
obj1 = []
obj2 = []
for solution in n_dominated:
obj1.append(solution.objectives[0])
obj2.append(solution.objectives[1])
var3.append(solution.variables[2]-1.0)
#
ax3.bar3d(obj1, obj2, 1.0, 0.01, 0.01, var3, color='red', label='Pareto front3' )
ax3.legend([obj1,obj2],["E","F"])


# パレート解をプロットする。(Dummy Fig. for plotting the 3rd Fig. normally)

var3= []
obj1 = []
obj2 = []
for solution in n_dominated:
obj1.append(solution.objectives[0])
obj2.append(solution.objectives[1])
var3.append(solution.variables[2])
#
# 実現可能解プロット X,Y 逆　to ax1.bar3d--ax3.bar3d
obj1 = []
obj2 = []
for solution in algorithm.result:
obj1.append(solution.objectives[0])
obj2.append(solution.objectives[1])
ax4.scatter(obj2, obj1,5.0, color='black', edgecolor='black', label='Solution')
# パレート解プロット
n_dominated = nondominated(algorithm.result) # 非劣解
obj1 = []
obj2 = []
for solution in n_dominated:
obj1.append(solution.objectives[0])
obj2.append(solution.objectives[1])
line1="Variables={0:.5f} {1:.5f} {2:.5f} Objectives={3:.5f} {4:.5f} {5:.5f}".\
format(solution.variables[0], solution.variables[1],solution.variables[2]\
,solution.objectives[0], solution.objectives[1],solution.constraints[0])
print(line1)

ax4.scatter(obj2, obj1,5.0, color='red', edgecolor='red', label='Pareto front')
ax4.legend()


#dummy ax4 for 3dplot
ax4.bar3d(obj1, obj2, 0.0, 0.01, 0.01, var3, color='black', label='Pareto front3' )
ax4.legend()#(Dummy legend. for plotting the 3rd Fig. normally)


#
# グラフを表示する。
plt.show()
plt.close()




表示した左上図で、黒丸が実行可能解、赤丸がパレートフロントとなっている。その他の３つの図は、パレート解における３つの説明変数の値を表示している。
（ただし、3次元図での各解の数値は、プリントされた数値で確認する必要がある。）