[PYTHON] Sekundärplanungsmethode nach interner Punktmethode
Sekundäres Planungsproblem
Letztes Mal löste das lineare Programmierproblem mithilfe der internen Punktmethode (Automarkierungsmethode). Dieses Mal werde ich versuchen, das sekundäre Planungsproblem mit der internen Punktmethode auf die gleiche Weise zu lösen. Das quadratische Planungsproblem wird als ein Problem ausgedrückt, das den Vektor $ x $ wie unten gezeigt minimiert, wobei die Zielfunktion in der quadratischen Form und die Einschränkungen in der linearen Form ausgedrückt werden.
\min \frac{1}{2} {\bf x}^T {\bf D} {\bf x} + {\bf c}^T {\bf x}\\
subject \ \ \ {\bf A}{\bf x} \ge {\bf b}
Wenn hier die Dimension von $ {\ bf x} $ $ n $ ist und die Anzahl der Einschränkungen $ m $ ist, dann sind $ {\ bf x}, {\ bf c}, {\ bf b} $ $ n $ Dimensionsvektoren. $ {\ Bf A} $ ist also eine Matrix von $ n \ mal m $, und $ {\ bf D} $ ist eine positive Zielmatrix von $ n \ mal n $.
Interne Punktmethode bei sekundären Planungsproblemen
Dieses Mal wurde es als Methode zur Lösung des sekundären Planungsproblems mit der internen Punktmethode in Numerical Optimization eingeführt. Ich habe versucht, die Prädiktor-Modifikator-Methode zu implementieren. Es ist komplexer als das lineare Planungsproblem.
#!/usr/bin/python
#coding:utf-8
import numpy as np
EPS_ZERO = 1.0e-9
def select_alpha(alphas):
return 1.0 if alphas.size == 0 or np.min(alphas) > 1.0 else np.min(alphas)
def inv_vec(vec):
inv = np.matrix(zeros((vec.shape[0], 1)))
cond = np.where(np.abs(vec) > EPS_ZERO)
inv[cond] = 1.0 / vec[cond]
return inv
def less_zero(vec):
return vec <= -EPS_ZERO
def quadratic_programming(x, gmat, c, amat, b,
tau=0.5, eps=1.0e-3, nloop=30):
"""
Interne Punktmethode(Predictor-Collector)Lösen des sekundären Planungsproblems durch
object min z = 0.5 * x^T * G * x + cT * x
subject Ax >= b
"""
ndim = gmat.shape[0] #Anzahl der Dimensionen
ncnst = amat.shape[0] #Anzahl der Einschränkungen
x_idx, y_idx, l_idx = np.split(np.arange(0, ndim + 2 * ncnst),
[ndim, ndim + ncnst])
y_l_idx = np.r_[y_idx, l_idx]
zeros_d_c = np.zeros((ndim, ncnst))
zeros_c2 = np.zeros((ncnst, ncnst))
eye_c = np.identity(ncnst)
f = np.bmat([[-gmat * x - c],
[zeros((ncnst, 1))],
[amat * x - b]])
exmat = np.bmat([[gmat, zeros_d_c, -amat.T],
[zeros_d_c.T, zeros_c2, eye_c],
[-amat, eye_c, zeros_c2]])
delta_xyl0 = np.linalg.pinv(exmat) * f
abs_yl0 = np.abs(delta_xyl0[y_l_idx])
abs_yl0[abs_yl0 < 1.0] = 1.0
y, lmd = np.vsplit(abs_yl0, 2)
rd = gmat * x - amat.T * lmd + c
rp = amat * x - y - b
for n in range(nloop):
err = np.linalg.norm(rd)**2 + np.linalg.norm(rp)**2 + np.asscalar(y.T * lmd)
if err < eps:
break
# affine scaling step
iy = inv_vec(y)
f = np.bmat([[-rd],
[-lmd],
[rp]])
exmat = np.bmat([[gmat, zeros_d_c, -amat.T],
[zeros_d_c.T, np.diag(np.multiply(lmd, iy).A1), eye_c],
[-amat, eye_c, zeros_c2]])
d_xyl_aff = np.linalg.pinv(exmat) * f
#Berechnung von Alpha
y_l_cnb = np.bmat([[y],
[lmd]])
cnd = less_zero(d_xyl_aff[y_l_idx])
alpha_aff = select_alpha(-y_l_cnb[cnd] / d_xyl_aff[y_l_idx][cnd])
#Zentrale Parameter
mu = np.asscalar(y.T * lmd / ncnst)
mu_aff = np.asscalar((y + alpha_aff * d_xyl_aff[y_idx]).T * (lmd + alpha_aff * d_xyl_aff[l_idx]) / ncnst)
sig = (mu_aff / mu)**3
# corrector step
dl_y_aff = np.multiply(np.multiply(d_xyl_aff[y_idx],
d_xyl_aff[l_idx]),
iy)
f = np.bmat([[-rd],
[-lmd - dl_y_aff + sig * mu * iy],
[rp]])
d_xyl = np.linalg.pinv(exmat) * f
# alpha_pri, alpha_Doppelte Berechnung
cnd = less_zero(d_xyl)
alpha_pri = select_alpha(-tau * y[cnd[y_idx]] / d_xyl[y_idx][cnd[y_idx]])
alpha_dual = select_alpha(-tau * lmd[cnd[l_idx]] / d_xyl[l_idx][cnd[l_idx]])
alpha_hat = min([alpha_pri, alpha_dual])
# x, y,lmd Update
x, y, lmd = np.vsplit(np.r_[x, y, lmd] + alpha_hat * d_xyl,
[ndim, ndim + ncnst])
rp = (1.0 - alpha_pri) * rp
rd = (1.0 - alpha_dual) * rd + (alpha_pri - alpha_dual) * gmat * d_xyl[x_idx]
return x, y, lmd
if __name__ == "__main__":
c = np.matrix([[-2.0],
[-4.0]])
dmat = np.matrix([[2.0, 0.0],
[0.0, 2.0]])
amat = np.matrix([[1.0, 1.0],
[1.0, -1.0],
[-3.0, 1.0]])
b = np.matrix([[-0.5],
[1.0],
[-3.0]])
#Zeichnung
from pylab import *
ax = subplot(111, aspect='equal')
x = np.arange(-3.0, 3.01, 0.15)
y = np.arange(-3.0, 3.01, 0.15)
X,Y = meshgrid(x, y)
t = arange(-3.0, 3.01, 0.15)
func = lambda x, y : 0.5 * (dmat[0, 0] * x**2 + dmat[1, 1] * y**2) + c[0, 0] * x + c[1, 0] * y
const = [lambda x : -amat[i, 0] / amat[i, 1] * x + b[i, 0] / amat[i, 1] for i in range(amat.shape[0])]
Z = func(X, Y)
s = [const[i](t)foriinrange(amat.shape[0])]
pcolor(X, Y, Z)
for i in range(amat.shape[0]):
ax.plot(t, s[i], 'k')
x, y, lmd = quadratic_programming(np.matrix([[0.0], [0.0]]),
dmat, c, amat, b,
tau=0.5, eps=1.0e-3, nloop=25)
print x
ax.plot([x[0, 0]], [x[1, 0]],'yo')
axis([-3, 3, -3, 3])
show()
Wie beim letzten Mal werde ich versuchen, mit einem zweidimensionalen Problem zu zeichnen. Der Bereich der Lösung, der die Bedingung erfüllt, liegt innerhalb der drei Linien, und die Richtung, in der das Blau dunkler wird, ist der optimale Wert. Die gelben Punkte in der Grafik geben die endgültige optimale Lösung an, und es ist ersichtlich, dass die geeignete Lösung erhalten werden kann, selbst wenn der Anfangswert von außerhalb der Beschränkung angegeben wird.
Da der diesmal erstellte Code Python ist, ist die Berechnungsgeschwindigkeit langsam, aber ich dachte, es wäre etwas einfacher, den Algorithmus zu verstehen. In diesem Code haben wir beim Auffinden der Schritte eine pseudo-inverse Matrix verwendet, aber es scheint, dass sie durch Verwendung der konjugierten Gradientenmethode oder durch Verwendung der Spärlichkeit von "exmat" beschleunigt werden kann.
Recommended Posts