[PYTHON] Méthode de planification secondaire par méthode de point interne
Problème de planification secondaire
Dernière fois a résolu le problème de programmation linéaire en utilisant la méthode du point interne (méthode du marqueur de voiture). Cette fois, je vais essayer de résoudre le problème de planification secondaire en utilisant la méthode du point interne de la même manière. Le problème de planification quadratique est exprimé comme un problème qui minimise le vecteur $ x $ comme indiqué ci-dessous, dans lequel la fonction objectif est exprimée sous forme quadratique et les contraintes sont exprimées sous forme linéaire.
\min \frac{1}{2} {\bf x}^T {\bf D} {\bf x} + {\bf c}^T {\bf x}\\
subject \ \ \ {\bf A}{\bf x} \ge {\bf b}
Ici, si la dimension de $ {\ bf x} $ est $ n $ et le nombre de contraintes est $ m $, alors $ {\ bf x}, {\ bf c}, {\ bf b} $ est un vecteur de dimension $ n $. Ainsi, $ {\ bf A} $ est une matrice de $ n \ fois m $, et $ {\ bf D} $ est une matrice cible positive de $ n \ fois n $.
Méthode du point interne dans les problèmes de planification secondaires
Cette fois, il a été introduit dans Optimisation numérique comme méthode pour résoudre le problème de planification secondaire par la méthode du point interne. J'ai essayé d'implémenter la méthode du modificateur de prédicteur. Il est plus complexe que le problème de planification linéaire.
#!/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):
"""
Méthode du point interne(Predictor-Collector)Résoudre le problème de planification secondaire en
object min z = 0.5 * x^T * G * x + cT * x
subject Ax >= b
"""
ndim = gmat.shape[0] #Nombre de dimensions
ncnst = amat.shape[0] #Nombre de contraintes
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
#Calcul de l'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])
#Paramètres centraux
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_Double calcul
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,mise à jour lmd
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]])
#dessin
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()
Comme la dernière fois, je vais essayer de dessiner avec un problème à deux dimensions. La zone de la solution qui satisfait la contrainte se trouve dans les trois lignes et la direction dans laquelle le bleu devient plus foncé est la valeur optimale. Les points jaunes dans le graphique indiquent la solution optimale finale, et on peut voir que la solution appropriée peut être obtenue même si la valeur initiale est donnée de l'extérieur de la contrainte.
Puisque le code créé cette fois-ci est Python, la vitesse de calcul est lente, mais j'ai pensé qu'il serait un peu plus facile de comprendre l'algorithme. Dans ce code, nous avons utilisé une pseudo matrice inverse lors de la recherche des étapes, mais il semble qu'elle puisse être accélérée en utilisant la méthode du gradient conjugué ou en utilisant la parcimonie de ʻexmat`.
Recommended Posts