このページは、nAGライブラリのJupyterノートブックExampleの日本語翻訳版です。オリジナルのノートブックはインタラクティブに操作することができます。
クリスマスデモ: クリスマスツリーの飾り付け方法
クリスマスツリーの飾り付けは難しい場合があります！この問題に最適化がどのように役立つか見てみましょう。
# 関連するライブラリをインポートする
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from naginterfaces.base import utils
from naginterfaces.library import opt
from naginterfaces.library import mip1. ツリーの定義
まず、関数を最適化したい領域としてツリーを定義する必要があります。
クリスマスらしい色のスキームを定義することから始めましょう：
cdict_green = {'red': ((0, 0, 0),
(1, 0, 0)),
'green': ((0, 0, 0),
(1, 1, 1)),
'blue': ((0, 0, 0),
(1, 0, 0))}
cdict_red = {'red': ((0, 0.5, 0.5),
(1, 1, 1)),
'green': ((0, 0, 0),
(1, 0, 0)),
'blue': ((0, 0, 0),
(1, 0, 0))}
green_cm = LinearSegmentedColormap('green_cm', cdict_green)
red_cm = LinearSegmentedColormap('red_cm', cdict_red)次に、2次元の木を描画するのに役立つ関数をいくつか定義します
def objf_plot(X,Y):
return X**2 + X + np.cos(2*np.pi*Y)**2 + 0.5*Ydef plot_domain(orig,lb,ub,_cm, axis, ax):
xmax = orig - lb
X, Y = np.meshgrid(np.arange(-xmax,xmax,0.01), np.arange(lb,orig,0.01))
Z = objf_plot(X,Y)
mask = np.zeros_like(Z, dtype=bool)
row, col = X.shape
for i in range(row):
for j in range(col):
if X[i,j]>=0 and X[i,j]+Y[i,j]>orig:
mask[i,j] = 1
elif X[i,j]<=0 and Y[i,j]-X[i,j]>orig:
mask[i,j] = 1
elif Y[i,j] > ub:
mask[i,j] = 1
else:
mask[i,j] = 0
Z = np.ma.array(Z, mask=mask)
ax.axis(axis)
ax.contourf(X, Y, Z, cmap=green_cm, corner_mask=True)結果を見てみましょう：
# Jupyter の表示バックエンドを選択:
%matplotlib inline
fig = plt.figure()
ax = fig.add_subplot(111)
axis = [-3.5, 3.5, 1, 6.5]
# メインツリー
plot_domain(5.5, 4, 5, green_cm, axis, ax)
plot_domain(6, 5, 999, green_cm, axis, ax)
plot_domain(5, 2.5, 4, green_cm, axis, ax)
# 幹
X, Y = np.meshgrid(np.arange(-0.75,0.75,0.01), np.arange(2,2.5,0.01))
Z = objf_plot(X,Y)
ax.contourf(X, Y, Z, cmap=green_cm)
# 現在
X, Y = np.meshgrid(np.arange(-1.,0.25,0.01), np.arange(1.3,1.8,0.01))
Z = objf_plot(X,Y)
ax.contourf(X, Y, Z, cmap=red_cm)
plt.show()
一見して木の最小値を見つけるのは簡単ではありません: 1. ドメインが分断されている 2. 木の各部分は線形制約に従っていますが、異なる部分を単一の制約セットにマージする簡単な方法はありません
2. ドメインの分割
まずは単純な分割
解決策の1つは、木の異なる部分を個別に見ることです。例えば、木の頂上の最小値を求めたい場合、3つの線形制約で定義されるドメインがあります:
orig = 6
lb = 5
ub = 999
axis = [-1.5, 1.5, 4.5, 6.5]
fig = plt.figure()
ax = fig.add_subplot(111)
plot_domain(orig, lb, ub, green_cm, axis, ax)
plt.show()
この部分の最適化は簡単であるはずです：nAGの非線形プログラミングソルバーの1つを使用しましょう。
# ソルバーから異なるコールバック関数に渡されるユーザー構造体を定義する
# ここでは、ソルバーによって評価されたポイントを記録するために使用します
class usr_data:
def __init__(self):
self.npts = 0
self.pts = []
self.present = False# SQPソルバーのメインコールバック関数
def objfun_e04wd(mode, x, grad, _nstate, data):
if mode in (0, 2):
objf = x[0]**2 + x[0] + np.cos(2*np.pi*x[1])**2 + 0.5*x[1]
data.npts += 1
data.pts.append((x[0], x[1]))
if mode in (1, 2):
grad[0] = 2*x[0] + 1
grad[1] = -4*np.pi*np.sin(2*np.pi*x[1])*np.cos(2*np.pi*x[1]) + 0.5
return mode, objf, grad# モデルを定義する: ドメインは2つの線形制約と1つの境界制約によって制限されています
bigbnd = 1.0e20
n = 2
nclin = 2
ncnln = 0
a = [[1.,1.],[-1.,1.]]
bl = [-bigbnd,5,-bigbnd,-bigbnd]
bu = [bigbnd,bigbnd,6,6]
# ソルバー変数の残りを初期化します: ウォームスタート情報や2次導関数情報は必要ありません
istate = [0,0,0,0]
ccon = np.zeros(1)
cjac = np.zeros((0,1))
clamda = np.zeros(n+nclin+ncnln)
h = np.zeros((2,2))
iom = utils.FileObjManager(locus_in_output=False)
# 開始点を定義し、ユーザーデータを初期化する
xstart1 = [0.,6.]
x = xstart1
comm = opt.nlp2_init(io_manager=iom)
data = usr_data()
# モデルを解きます。解と収束の指標はslnに格納されます
sln = opt.nlp2_solve(a, bl, bu, objfun_e04wd, istate, ccon, cjac, clamda, h, x, comm, data=data, io_manager=iom)これで最適値がどこにあるかを発見し、ソルバーが評価したポイントでツリーの頂上を飾り始めることができます。
# 木の頂点と評価されたポイントの履歴をプロットする
fig = plt.figure()
ax = fig.add_subplot(111)
plot_domain(orig, lb, ub, green_cm, axis, ax)
ax.plot(sln.x[0], sln.x[1], marker='*', markersize=12)
ax.scatter(*zip(*data.pts), s=25, marker='x', c='red')
plt.show()
2つのドメインの結合
ここで、木の中央部分を追加したい場合、ドメインの定義はより複雑になります:
fig = plt.figure()
ax = fig.add_subplot(111)
axis = [-2., 2., 3.5, 6.5]
plot_domain(5.5, 4, 5, green_cm, axis, ax)
plot_domain(6, 5, 999, green_cm, axis, ax)
plt.show()
中央部分と頂部を同じモデルに結合することは、問題に追加の決定変数を導入することで可能です。目標は、各部分領域 \(K\) に対して二値変数 \(\delta_K\) を導入し、\(\delta_K=1 \Rightarrow (x,y) \in K\) となるようにすることです。一般に、\(ax+by \leq c\) の形の制約に対して、元の不等式を次のように置き換えることで、そのような決定変数 \(\delta\) を作成できます：\[ax+by+M\delta \leq M + c\] ここで、\(M\) は \(ax+by-c\) の上限です。これにより、\(\delta=1\) の場合は元の不等式が適用され、\(\delta=0\) の場合は不等式が \(ax+by-c \leq M\) となり、常に満たされるため影響を与えません。
各部分領域は3つの線形不等式で定義されているため、この手法をすべての制約に適用し、2つの二値変数 \(\delta_{top}\) と \(\delta_{mid} \in \{0,1\}\) を導入できます。モデルは以下のように記述できます： \[ \begin{equation*} \min_{(x,y)\in \mathbb{R}^n} F(x,y)\\ \text{s.t. } \left\{ \begin{aligned} x + y + \delta_{top}M &\leq M + 6\\ -x + y + \delta_{top}M &\leq M + 6\\ y - \delta_{top}M &\geq -M + 5\\ x + y + \delta_{mid}M &\leq M + 5.5\\ -x + y + \delta_{mid}M &\leq M + 5.5\\ y - \delta_{mid}M &\geq -M + 4\\ \end{aligned}\right. \end{equation*} \] 簡単のため、これらすべての不等式に対する上限として単一の \(M=9\) を選択しました。
現状では、モデルはまだ不完全です。ソルバーは \(\delta_{mid}=\delta_{top}=0\) を選択し、すべての制約を無視する可能性があります。モデルに少なくとも1つの領域内の点を選択させるために、以下の拘束制約を追加する必要があります： \[ \delta_{mid}+\delta_{top} \geq 1 \] これは、少なくとも1つの制約セットが満たされなければならないことを指定しています。
このモデルは現在二値変数を含むため、古典的なNLPソルバーはもはや適していません。我々はこのモデルをnAGのMINLPソルバーmip.sqp（h02da）で解きました。
# MINLPソルバーのメインコールバック関数
def objfun_h02da(mode, _varcon, x, objgrd, _nstate, data):
objmip = x[0]**2 + x[0] + np.cos(2*np.pi*x[1])**2 + 0.5*x[1]
data.npts = data.npts + 1
data.pts.append((x[0], x[1]))
if x[1] <= 1.8 and data.present:
objmip = (x[0]+0.625)**2 + (x[1]-1.55)**2
if mode == 1:
objgrd[0] = 2*x[0] + 1
objgrd[1] = -4*np.pi*np.sin(2*np.pi*x[1])*np.cos(2*np.pi*x[1]) + 0.5
if x[1] <= 1.8 and data.present:
objgrd[0] = 2*(x[0]+0.625)
objgrd[1] = 2*(x[1]-1.55)
return objmip, objgrd# モデルを定義する
ncnln = 0
M = 9.
a = np.array([[0.,0.,1.,1.],
[1.,-1.,-M,0.],
[-1.,-1.,-M,0.],
[0.,1.,-M,0.],
[1.,-1.,0.,-M],
[-1.,-1.,0.,-M],
[0.,1.,0.,-M]])
d = [1.,-6.-M,-6.-M,-M+5.,-M-5.5,-M-5.5,-M+4.]
bl = [-5.,-2.,0.,0.]
bu = [5.,10.,1.,1.]
varcon = [0,0,1,1,4,4,4,4,4,4,4]
# 開始点を選択し、ソルバーとユーザーデータを初期化する
xstart = [0.5,4.5,0.,0.]
x = xstart
comm = {}
data = usr_data()
mip.optset('Initialize = sqp', comm)
# モデルを解く
sln = mip.sqp(ncnln, bl, bu, varcon, x, objfun_h02da, comm, a = a, d=d, maxit=500, acc=0.0001, data=data, io_manager=iom)# ドメインをプロットし、評価されたポイントを追加する
fig = plt.figure()
ax = fig.add_subplot(111)
plot_domain(5.5, 4, 5, green_cm, axis, ax)
plot_domain(6, 5, 999, green_cm, axis, ax)
ax.plot(sln.x[0], sln.x[1], marker='*', markersize=12)
ax.scatter(*zip(*data.pts), s=25, marker='x', c='red')
plt.show()
完全なモデル
これで、より分離した部分領域に同じテクニックを使用して、完全なモデルを考えることができます： \[ \begin{equation*} \min_{(x,y)\in \mathbb{R}^n} F(x,y)\\ \text{s.t. } \left\{ \begin{aligned} \left. \delta_{top} + \delta_{mid} + \delta_{bot} + \delta_{tr} + \delta_{pr} \geq 1 \right\}\text{binding constraint}\\ \left.\begin{aligned} x + y + \delta_{top}M &\leq M + 6\\ -x + y + \delta_{top}M &\leq M + 6\\ y - \delta_{top}M &\geq -M + 5\\ \end{aligned}\right\}\text{top of the tree}\\ \left.\begin{aligned} x + y + \delta_{mid}M &\leq M + 5.5\\ -x + y + \delta_{mid}M &\leq M + 5.5\\ y - \delta_{mid}M &\geq -M + 4\\ \end{aligned}\right\}\text{middle of the tree}\\ \left.\begin{aligned} x + y + \delta_{bot}M &\leq M + 5\\ -x + y + \delta_{bot}M &\leq M + 5\\ y - \delta_{bot}M &\geq -M + 2.5\\ \end{aligned}\right\}\text{bottom of the tree}\\ \left.\begin{aligned} x - \delta_{tr}M &\geq -M - 0.75\\ x + \delta_{tr}M &\leq M + 0.75\\ y - \delta_{tr}M &\geq -M + 2\\ y + \delta_{tr}M &\leq M + 2.5\\ \end{aligned}\right\}\text{trunk}\\ \left.\begin{aligned} x - \delta_{pr}M &\geq -M - 1\\ x + \delta_{pr}M &\leq M - 0.25\\ y - \delta_{pr}M &\geq -M + 1.8\\ y + \delta_{pr}M &\leq M + 1.3\\ \end{aligned}\right\}\text{present}\\ \end{aligned}\right. \end{equation*} \]
これはローカルソルバーであり、開始点に応じて領域内の任意のローカル最小値で収束を宣言する可能性があることに注意することが依然として重要です。
# モデルを拡張する
row, col = a.shape
a2 = np.zeros((row, col+3))
a2[:,:-3] = a
a2[0,4:6] = 1.
add_rows = [[1.,-1.,0.,0.,-M,0.,0.],
[-1.,-1.,0.,0.,-M,0.,0.],
[0.,1.,0.,0.,-M,0.,0.],
[0.,1.,0.,0.,0.,-M,0.],
[0.,-1.,0.,0.,0.,-M,0.],
[1.,0.,0.,0.,0.,-M,0.],
[-1.,0.,0.,0.,0.,-M,0.],
[0.,1.,0.,0.,0.,0.,-M],
[0.,-1.,0.,0.,0.,0.,-M],
[1.,0.,0.,0.,0.,0.,-M],
[-1.,0.,0.,0.,0.,0.,-M]]
a2 = np.vstack((a2, add_rows))
d.extend([-M-5,-M-5,-M+2.5])
d.extend([-M+2,-M-2.5,-M-0.75,-M-0.75])
d.extend([-M+1.3,-M-1.8,-M-1,-M+0.25])
bl.extend([0.,0.,0.])
bu.extend([0.,0.,0.])# モデルを最適化する
xstart = [1.,2.5,0.,0.,1.,0.,0.]
x = xstart
varcon = [0,0,1,1,1,1,1,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4]
comm = {}
data = usr_data()
mip.optset('Initialize = sqp', comm)
sln = mip.sqp(ncnln, bl, bu, varcon, x, objfun_h02da, comm, a=a2, d=d,maxit=500, acc=0.0001, data=data,
io_manager=iom)# ツリーと評価されたポイントをプロットする
fig = plt.figure()
ax = fig.add_subplot(111)
axis = [-3.5, 3.5, 0.5, 6.5]
# メインツリー
plot_domain(5.5, 4, 5, green_cm, axis, ax)
plot_domain(6, 5, 999, green_cm, axis, ax)
plot_domain(5, 2.5, 4, green_cm, axis, ax)
# 幹
X, Y = np.meshgrid(np.arange(-0.75,0.75,0.01), np.arange(2,2.5,0.01))
Z = objf_plot(X,Y)
ax.contourf(X, Y, Z, cmap=green_cm)
# 現在
X, Y = np.meshgrid(np.arange(-1.,0.25,0.01), np.arange(1.3,1.8,0.01))
Z = objf_plot(X,Y)
ax.contourf(X, Y, Z, cmap=red_cm)
ax.plot(sln.x[0], sln.x[1], marker='*', markersize=12)
ax.scatter(*zip(*data.pts), s=25, marker='x', c='red')
plt.show()