You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
419 lines
14 KiB
419 lines
14 KiB
|
|
import numpy as np
|
|
import scipy.sparse
|
|
import scipy.sparse.linalg
|
|
|
|
def mmasub(m,n,iter,xval,xmin,xmax,xold1,xold2,
|
|
f0val,df0dx,fval,dfdx,low,upp,a0,a,c,d):
|
|
|
|
# Version September 2007 (and a small change August 2008)
|
|
#
|
|
# Krister Svanberg <krille@math.kth.se>
|
|
# Department of Mathematics, SE-10044 Stockholm, Sweden.
|
|
#
|
|
# This function mmasub performs one MMA-iteration, aimed at
|
|
# solving the nonlinear programming problem:
|
|
#
|
|
# Minimize f_0(x) + a_0*z + sum( c_i*y_i + 0.5*d_i*(y_i)^2 )
|
|
# subject to f_i(x) - a_i*z - y_i <= 0, i = 1,...,m
|
|
# xmin_j <= x_j <= xmax_j, j = 1,...,n
|
|
# z >= 0, y_i >= 0, i = 1,...,m
|
|
# *** INPUT:
|
|
#
|
|
# m = The number of general constraints.
|
|
# n = The number of variables x_j.
|
|
# iter = Current iteration number ( =1 the first time mmasub is called).
|
|
# xval = Column vector with the current values of the variables x_j.
|
|
# xmin = Column vector with the lower bounds for the variables x_j.
|
|
# xmax = Column vector with the upper bounds for the variables x_j.
|
|
# xold1 = xval, one iteration ago (provided that iter>1).
|
|
# xold2 = xval, two iterations ago (provided that iter>2).
|
|
# f0val = The value of the objective function f_0 at xval.
|
|
# df0dx = Column vector with the derivatives of the objective function
|
|
# f_0 with respect to the variables x_j, calculated at xval.
|
|
# fval = Column vector with the values of the constraint functions f_i,
|
|
# calculated at xval.
|
|
# dfdx = (m x n)-matrix with the derivatives of the constraint functions
|
|
# f_i with respect to the variables x_j, calculated at xval.
|
|
# dfdx(i,j) = the derivative of f_i with respect to x_j.
|
|
# low = Column vector with the lower asymptotes from the previous
|
|
# iteration (provided that iter>1).
|
|
# upp = Column vector with the upper asymptotes from the previous
|
|
# iteration (provided that iter>1).
|
|
# a0 = The constants a_0 in the term a_0*z.
|
|
# a = Column vector with the constants a_i in the terms a_i*z.
|
|
# c = Column vector with the constants c_i in the terms c_i*y_i.
|
|
# d = Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
|
|
#
|
|
# *** OUTPUT:
|
|
#
|
|
# xmma = Column vector with the optimal values of the variables x_j
|
|
# in the current MMA subproblem.
|
|
# ymma = Column vector with the optimal values of the variables y_i
|
|
# in the current MMA subproblem.
|
|
# zmma = Scalar with the optimal value of the variable z
|
|
# in the current MMA subproblem.
|
|
# lam = Lagrange multipliers for the m general MMA constraints.
|
|
# xsi = Lagrange multipliers for the n constraints alfa_j - x_j <= 0.
|
|
# eta = Lagrange multipliers for the n constraints x_j - beta_j <= 0.
|
|
# mu = Lagrange multipliers for the m constraints -y_i <= 0.
|
|
# zet = Lagrange multiplier for the single constraint -z <= 0.
|
|
# s = Slack variables for the m general MMA constraints.
|
|
# low = Column vector with the lower asymptotes, calculated and used
|
|
# in the current MMA subproblem.
|
|
# upp = Column vector with the upper asymptotes, calculated and used
|
|
# in the current MMA subproblem.
|
|
|
|
epsimin = 1e-10
|
|
raa0=0.01
|
|
albefa=0.4
|
|
asyinit=0.1
|
|
asyincr=0.8
|
|
asydecr=0.6
|
|
|
|
# epsimin = sqrt(m+n)*10^(-9);
|
|
#epsimin = 1e-7
|
|
#raa0 = 1e-5
|
|
move = 1.0
|
|
# move = 0.05 # MODIFICATION BY DK
|
|
#albefa = 0.1
|
|
#asyinit = 0.5
|
|
#asyincr = 1.2
|
|
#asydecr = 0.7
|
|
# asyincr = 1.05 # MODIFICATION BY DK
|
|
# asydecr = 0.65 # MODIFICATION BY DK
|
|
eeen = np.ones((n,1))
|
|
eeem = np.ones((m,1))
|
|
zeron = np.zeros((n,1))
|
|
|
|
# Calculation of the asymptotes low and upp :
|
|
if iter < 2.5:
|
|
low = xval - asyinit*(xmax-xmin)
|
|
upp = xval + asyinit*(xmax-xmin)
|
|
else:
|
|
zzz = (xval-xold1)*(xold1-xold2)
|
|
factor = eeen
|
|
factor[zzz > 0] = asyincr
|
|
factor[zzz < 0] = asydecr
|
|
low = xval - factor*(xold1 - low)
|
|
upp = xval + factor*(upp - xold1)
|
|
lowmin = xval - 10*(xmax-xmin)
|
|
lowmax = xval - 0.01*(xmax-xmin)
|
|
uppmin = xval + 0.01*(xmax-xmin)
|
|
uppmax = xval + 10*(xmax-xmin)
|
|
low = np.maximum(low,lowmin)
|
|
low = np.minimum(low,lowmax)
|
|
upp = np.minimum(upp,uppmax)
|
|
upp = np.maximum(upp,uppmin)
|
|
|
|
# Calculation of the bounds alfa and beta :
|
|
zzz1 = low + albefa*(xval-low)
|
|
zzz2 = xval - move*(xmax-xmin)
|
|
zzz = np.maximum(zzz1,zzz2)
|
|
alfa = np.maximum(zzz,xmin)
|
|
zzz1 = upp - albefa*(upp-xval)
|
|
zzz2 = xval + move*(xmax-xmin)
|
|
zzz = np.minimum(zzz1,zzz2)
|
|
beta = np.minimum(zzz,xmax)
|
|
|
|
# Calculations of p0, q0, P, Q and b.
|
|
xmami = xmax-xmin
|
|
xmamieps = 0.00001*eeen
|
|
xmami = np.maximum(xmami,xmamieps)
|
|
xmamiinv = eeen/xmami
|
|
#print(xmamiinv.shape) #112 x 1
|
|
ux1 = upp-xval
|
|
ux2 = ux1*ux1
|
|
xl1 = xval-low
|
|
xl2 = xl1*xl1
|
|
uxinv = eeen/ux1
|
|
xlinv = eeen/xl1
|
|
|
|
# p0 = zeron
|
|
# q0 = zeron
|
|
|
|
p0 = np.maximum(df0dx,0)
|
|
q0 = np.maximum(-df0dx,0)
|
|
|
|
# p0(find(df0dx > 0)) = df0dx(find(df0dx > 0))
|
|
# q0(find(df0dx < 0)) = -df0dx(find(df0dx < 0))
|
|
|
|
pq0 = 0.001*(p0 + q0) + raa0*xmamiinv
|
|
p0 = p0 + pq0
|
|
q0 = q0 + pq0
|
|
p0 = p0*ux2
|
|
q0 = q0*xl2
|
|
|
|
#P = coo_matrix((m,n))
|
|
#Q = coo_matrix((m,n))
|
|
|
|
P = np.maximum(dfdx,0)
|
|
Q = np.maximum(-dfdx,0)
|
|
|
|
# P(find(dfdx > 0)) = dfdx(find(dfdx > 0))
|
|
# Q(find(dfdx < 0)) = -dfdx(find(dfdx < 0))
|
|
#print("--",(eeem @ xmamiinv.T).shape)
|
|
PQ = 0.001*(P + Q) + raa0 * eeem @ xmamiinv.T
|
|
#print(PQ.shape)
|
|
P = P + PQ
|
|
Q = Q + PQ
|
|
P = P @ scipy.sparse.spdiags(ux2.T,0,n,n)
|
|
Q = Q @ scipy.sparse.spdiags(xl2.T,0,n,n)
|
|
b = P @ uxinv + Q @ xlinv - fval
|
|
#print("--",P.shape)
|
|
|
|
# %%% Solving the subproblem by a primal-dual Newton method
|
|
xmma,ymma,zmma,lam,xsi,eta,mu,zet,s = subsolve(m,n,epsimin,low,upp,alfa,
|
|
beta,p0,q0,P,Q,a0,a,b,c,d)
|
|
|
|
return xmma, ymma, zmma, lam, xsi, eta, mu, zet, s, low, upp
|
|
|
|
|
|
def subsolve(m,n,epsimin,low,upp,alfa,beta,p0,q0,P,Q,a0,a,b,c,d):
|
|
# -------------------------------------------------------------
|
|
# This is the file subsolv.m
|
|
#
|
|
# Version Dec 2006.
|
|
# Krister Svanberg <krille@math.kth.se>
|
|
# Department of Mathematics, KTH,
|
|
# SE-10044 Stockholm, Sweden.
|
|
#
|
|
#
|
|
# This function subsolv solves the MMA subproblem:
|
|
#
|
|
# minimize SUM[ p0j/(uppj-xj) + q0j/(xj-lowj) ] + a0*z +
|
|
# + SUM[ ci*yi + 0.5*di*(yi)^2 ],
|
|
#
|
|
# subject to SUM[ pij/(uppj-xj) + qij/(xj-lowj) ] - ai*z - yi <= bi,
|
|
# alfaj <= xj <= betaj, yi >= 0, z >= 0.
|
|
#
|
|
# Input: m, n, low, upp, alfa, beta, p0, q0, P, Q, a0, a, b, c, d.
|
|
# Output: xmma,ymma,zmma, slack variables and Lagrange multiplers.
|
|
#
|
|
een = np.ones((n,1))
|
|
eem = np.ones((m,1))
|
|
epsi = 1
|
|
# epsvecn = epsi*een
|
|
# epsvecm = epsi*eem
|
|
x = 0.5*(alfa+beta)
|
|
y = eem
|
|
z = 1
|
|
lam = eem
|
|
xsi = een/(x-alfa)
|
|
xsi = np.maximum(xsi,een)
|
|
eta = een/(beta-x)
|
|
eta = np.maximum(eta,een)
|
|
mu = np.maximum(eem,0.5*c)
|
|
zet = 1
|
|
s = eem
|
|
itera = 0
|
|
|
|
while epsi > epsimin:
|
|
|
|
epsvecn = epsi*een
|
|
epsvecm = epsi*eem
|
|
|
|
ux1 = upp-x
|
|
xl1 = x-low
|
|
ux2 = ux1*ux1
|
|
xl2 = xl1*xl1
|
|
|
|
uxinv1 = een/ux1
|
|
xlinv1 = een/xl1
|
|
|
|
#print(p0.shape)
|
|
#print(P.T.shape)
|
|
#print(lam.shape)
|
|
plam = p0 + P.T @ lam
|
|
qlam = q0 + Q.T @ lam
|
|
|
|
gvec = P @ uxinv1 + Q @ xlinv1
|
|
dpsidx = plam/ux2 - qlam/xl2
|
|
|
|
rex = dpsidx - xsi + eta
|
|
rey = c + d*y - mu - lam
|
|
rez = a0 - zet - a.T @ lam
|
|
|
|
relam = gvec - a*z - y + s - b
|
|
rexsi = xsi*(x-alfa) - epsvecn
|
|
reeta = eta*(beta-x) - epsvecn
|
|
remu = mu*y - epsvecm
|
|
rezet = zet*z - epsi
|
|
res = lam*s - epsvecm
|
|
|
|
residu1 = np.concatenate((rex,rey,rez), axis = 0)
|
|
residu2 = np.concatenate((relam, rexsi, reeta, remu, np.array([[rezet]]), res), axis = 0)
|
|
residu = np.concatenate((residu1, residu2), axis = 0)
|
|
residunorm = np.linalg.norm(residu)
|
|
residumax = np.amax(np.absolute(residu))
|
|
|
|
|
|
ittt = 0
|
|
|
|
while (residumax > 0.9*epsi) and (ittt < 200):
|
|
|
|
ittt = ittt + 1
|
|
itera = itera + 1
|
|
|
|
ux1 = upp-x
|
|
xl1 = x-low
|
|
ux2 = ux1*ux1
|
|
xl2 = xl1*xl1
|
|
ux3 = ux1*ux2
|
|
xl3 = xl1*xl2
|
|
|
|
uxinv1 = een/ux1
|
|
xlinv1 = een/xl1
|
|
uxinv2 = een/ux2
|
|
xlinv2 = een/xl2
|
|
|
|
plam = p0 + P.T @ lam
|
|
qlam = q0 + Q.T @ lam
|
|
|
|
gvec = P @ uxinv1 + Q @ xlinv1
|
|
GG = P @ scipy.sparse.spdiags(uxinv2.T,0,n,n) - Q @ scipy.sparse.spdiags(xlinv2.T,0,n,n)
|
|
dpsidx = plam/ux2 - qlam/xl2
|
|
# print('line 261, ittt=',ittt)
|
|
# print(np.amin(np.absolute(x-alfa)))
|
|
# print(np.amin(np.absolute(beta-x)))
|
|
# print()
|
|
|
|
delx = dpsidx - epsvecn/(x-alfa) + epsvecn/(beta-x)
|
|
dely = c + d*y - lam - epsvecm/y
|
|
delz = a0 - a.T@lam - epsi/z
|
|
dellam = gvec - a*z - y - b + epsvecm/lam
|
|
diagx = plam/ux3 + qlam/xl3
|
|
diagx = 2*diagx + xsi/(x-alfa) + eta/(beta-x)
|
|
diagxinv = een/diagx
|
|
diagy = d + mu/y
|
|
diagyinv = eem/diagy
|
|
diaglam = s/lam
|
|
diaglamyi = diaglam + diagyinv
|
|
|
|
if m < n:
|
|
blam = dellam + dely/diagy - GG @ (delx/diagx)
|
|
bb = np.concatenate((blam, delz), axis = 0)
|
|
Alam = scipy.sparse.spdiags(diaglamyi.T,0,m,m) + GG @ scipy.sparse.spdiags(diagxinv.T,0,n,n) @ GG.T
|
|
AA = np.concatenate((np.concatenate((Alam, a), axis = 1),
|
|
np.concatenate((a.T, np.array([[-zet/z]])), axis = 1)), axis = 0)
|
|
solut = scipy.sparse.linalg.spsolve(AA,bb)[:,np.newaxis]
|
|
dlam = solut[0:m]
|
|
dz = solut[m,0]
|
|
dx = -delx/diagx - (GG.T @ dlam)/diagx
|
|
|
|
|
|
else:
|
|
diaglamyiinv = eem/diaglamyi
|
|
dellamyi = dellam + dely/diagy
|
|
Axx = scipy.sparse.spdiags(diagx.T,0,n,n) + GG.T @ scipy.sparse.spdiags(diaglamyiinv.T,0,m,m) @ GG
|
|
azz = zet/z + a.T @ (a/diaglamyi)
|
|
axz = -GG.T @ (a/diaglamyi)
|
|
bx = delx + GG.T @ (dellamyi/diaglamyi)
|
|
bz = delz - a.T @ (dellamyi/diaglamyi)
|
|
AA = np.concatenate((np.concatenate((Axx, axz), axis = 1),
|
|
np.concatenate((axz.T, azz), axis = 1)), axis = 0)
|
|
bb = np.concatenate((-bx, -bz),axis = 0)
|
|
solut = scipy.sparse.linalg.spsolve(AA,bb)
|
|
dx = solut[0:n]
|
|
dz = solut[n,0]
|
|
dlam = (GG @ dx)/diaglamyi - dz*(a/diaglamyi) + dellamyi/diaglamyi
|
|
|
|
|
|
dy = -dely/diagy + dlam/diagy
|
|
dxsi = -xsi + epsvecn/(x-alfa) - (xsi*dx)/(x-alfa)
|
|
deta = -eta + epsvecn/(beta-x) + (eta*dx)/(beta-x)
|
|
dmu = -mu + epsvecm/y - (mu*dy)/y
|
|
dzet = -zet + epsi/z - zet*dz/z
|
|
ds = -s + epsvecm/lam - (s*dlam)/lam
|
|
xx = np.concatenate((y, np.array([[z]]), lam, xsi, eta, mu, np.array([[zet]]), s), axis = 0)
|
|
dxx = np.concatenate((dy, np.array([[dz]]), dlam, dxsi, deta, dmu, np.array([[dzet]]), ds), axis = 0)
|
|
|
|
|
|
stepxx = -1.01*dxx/xx
|
|
stmxx = np.amax(stepxx)
|
|
stepalfa = -1.01*dx/(x-alfa)
|
|
stmalfa = np.amax(stepalfa)
|
|
stepbeta = 1.01*dx/(beta-x)
|
|
stmbeta = np.amax(stepbeta)
|
|
stmalbe = np.maximum(stmalfa,stmbeta)
|
|
stmalbexx = np.maximum(stmalbe,stmxx)
|
|
stminv = np.maximum(stmalbexx,1)
|
|
steg = 1/stminv
|
|
|
|
xold = x
|
|
yold = y
|
|
zold = z
|
|
lamold = lam
|
|
xsiold = xsi
|
|
etaold = eta
|
|
muold = mu
|
|
zetold = zet
|
|
sold = s
|
|
|
|
itto = 0
|
|
resinew = 2*residunorm
|
|
|
|
while (resinew > residunorm) and (itto < 50):
|
|
itto = itto+1
|
|
x = xold + steg*dx
|
|
y = yold + steg*dy
|
|
z = zold + steg*dz
|
|
lam = lamold + steg*dlam
|
|
xsi = xsiold + steg*dxsi
|
|
eta = etaold + steg*deta
|
|
mu = muold + steg*dmu
|
|
zet = zetold + steg*dzet
|
|
s = sold + steg*ds
|
|
|
|
ux1 = upp-x
|
|
xl1 = x-low
|
|
ux2 = ux1*ux1
|
|
xl2 = xl1*xl1
|
|
|
|
uxinv1 = een/ux1
|
|
xlinv1 = een/xl1
|
|
plam = p0 + P.T @ lam
|
|
qlam = q0 + Q.T @ lam
|
|
gvec = P @ uxinv1 + Q @ xlinv1
|
|
|
|
dpsidx = plam/ux2 - qlam/xl2
|
|
rex = dpsidx - xsi + eta
|
|
rey = c + d*y - mu - lam
|
|
rez = a0 - zet - a.T @ lam
|
|
relam = gvec - a*z - y + s - b
|
|
rexsi = xsi*(x-alfa) - epsvecn
|
|
reeta = eta*(beta-x) - epsvecn
|
|
remu = mu*y - epsvecm
|
|
rezet = zet*z - epsi
|
|
res = lam*s - epsvecm
|
|
|
|
residu1 = np.concatenate((rex, rey, rez), axis = 0)
|
|
residu2 = np.concatenate((relam, rexsi, reeta, remu, np.array([[rezet]]), res), axis = 0)
|
|
residu = np.concatenate((residu1, residu2), axis = 0)
|
|
resinew = np.linalg.norm(residu)
|
|
|
|
steg = steg/2
|
|
|
|
residunorm = resinew
|
|
residumax = np.amax(np.absolute(residu))
|
|
steg = 2*steg
|
|
|
|
if ittt > 198:
|
|
print(epsi)
|
|
print(ittt)
|
|
|
|
epsi = 0.1*epsi
|
|
|
|
|
|
# store final output
|
|
xmma = x
|
|
ymma = y
|
|
zmma = z
|
|
lamma = lam
|
|
xsimma = xsi
|
|
etamma = eta
|
|
mumma = mu
|
|
zetmma = zet
|
|
smma = s
|
|
|
|
return xmma, ymma, zmma, lamma, xsimma, etamma, mumma, zetmma, smma
|
|
|