pypower/pips.py

   1 # Copyright (c) 1996-2015 PSERC. All rights reserved.
   2 # Use of this source code is governed by a BSD-style
   3 # license that can be found in the LICENSE file.
   4
   5 """Python Interior Point Solver (PIPS).
   6 """
   7
   8 from numpy import array, Inf, any, isnan, ones, r_, finfo, \
   9     zeros, dot, absolute, log, flatnonzero as find
  10
  11 from numpy.linalg import norm
  12
  13 from scipy.sparse import vstack, hstack, eye, csr_matrix as sparse
  14
  15 from pypower.pipsver import pipsver
  16 from pypower.pplinsolve import pplinsolve
  17
  18
  19 EPS = finfo(float).eps
  20
  21
  22 def pips(f_fcn, x0=None, A=None, l=None, u=None, xmin=None, xmax=None,
  23          gh_fcn=None, hess_fcn=None, opt=None):
  24     """Primal-dual interior point method for NLP (nonlinear programming).
  25     Minimize a function F(X) beginning from a starting point M{x0}, subject to
  26     optional linear and nonlinear constraints and variable bounds::
  27
  28             min f(x)
  29              x
  30
  31     subject to::
  32
  33             g(x) = 0            (nonlinear equalities)
  34             h(x) <= 0           (nonlinear inequalities)
  35             l <= A*x <= u       (linear constraints)
  36             xmin <= x <= xmax   (variable bounds)
  37
  38     Note: The calling syntax is almost identical to that of FMINCON from
  39     MathWorks' Optimization Toolbox. The main difference is that the linear
  40     constraints are specified with C{A}, C{L}, C{U} instead of C{A}, C{B},
  41     C{Aeq}, C{Beq}. The functions for evaluating the objective function,
  42     constraints and Hessian are identical.
  43
  44     Example from U{http://en.wikipedia.org/wiki/Nonlinear_programming}:
  45         >>> from numpy import array, r_, float64, dot
  46         >>> from scipy.sparse import csr_matrix
  47         >>> def f2(x):
  48         ...     f = -x[0] * x[1] - x[1] * x[2]
  49         ...     df = -r_[x[1], x[0] + x[2], x[1]]
  50         ...     # actually not used since 'hess_fcn' is provided
  51         ...     d2f = -array([[0, 1, 0], [1, 0, 1], [0, 1, 0]], float64)
  52         ...     return f, df, d2f
  53         >>> def gh2(x):
  54         ...     h = dot(array([[1, -1, 1],
  55         ...                    [1,  1, 1]]), x**2) + array([-2.0, -10.0])
  56         ...     dh = 2 * csr_matrix(array([[ x[0], x[0]],
  57         ...                                [-x[1], x[1]],
  58         ...                                [ x[2], x[2]]]))
  59         ...     g = array([])
  60         ...     dg = None
  61         ...     return h, g, dh, dg
  62         >>> def hess2(x, lam, cost_mult=1):
  63         ...     mu = lam["ineqnonlin"]
  64         ...     a = r_[dot(2 * array([1, 1]), mu), -1, 0]
  65         ...     b = r_[-1, dot(2 * array([-1, 1]), mu),-1]
  66         ...     c = r_[0, -1, dot(2 * array([1, 1]), mu)]
  67         ...     Lxx = csr_matrix(array([a, b, c]))
  68         ...     return Lxx
  69         >>> x0 = array([1, 1, 0], float64)
  70         >>> solution = pips(f2, x0, gh_fcn=gh2, hess_fcn=hess2)
  71         >>> round(solution["f"], 11) == -7.07106725919
  72         True
  73         >>> solution["output"]["iterations"]
  74         8
  75
  76     Ported by Richard Lincoln from the MATLAB Interior Point Solver (MIPS)
  77     (v1.9) by Ray Zimmerman.  MIPS is distributed as part of the MATPOWER
  78     project, developed at the Power System Engineering Research Center (PSERC) (PSERC),
  79     Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more info.
  80     MIPS was ported by Ray Zimmerman from C code written by H. Wang for his
  81     PhD dissertation:
  82       - "On the Computation and Application of Multi-period
  83         Security-Constrained Optimal Power Flow for Real-time
  84         Electricity Market Operations", Cornell University, May 2007.
  85
  86     See also:
  87       - H. Wang, C. E. Murillo-Sanchez, R. D. Zimmerman, R. J. Thomas,
  88         "On Computational Issues of Market-Based Optimal Power Flow",
  89         IEEE Transactions on Power Systems, Vol. 22, No. 3, Aug. 2007,
  90         pp. 1185-1193.
  91
  92     All parameters are optional except C{f_fcn} and C{x0}.
  93     @param f_fcn: Function that evaluates the objective function, its gradients
  94                   and Hessian for a given value of M{x}. If there are
  95                   nonlinear constraints, the Hessian information is provided
  96                   by the 'hess_fcn' argument and is not required here.
  97     @type f_fcn: callable
  98     @param x0: Starting value of optimization vector M{x}.
  99     @type x0: array
 100     @param A: Optional linear constraints.
 101     @type A: csr_matrix
 102     @param l: Optional linear constraints. Default values are M{-Inf}.
 103     @type l: array
 104     @param u: Optional linear constraints. Default values are M{Inf}.
 105     @type u: array
 106     @param xmin: Optional lower bounds on the M{x} variables, defaults are
 107                  M{-Inf}.
 108     @type xmin: array
 109     @param xmax: Optional upper bounds on the M{x} variables, defaults are
 110                  M{Inf}.
 111     @type xmax: array
 112     @param gh_fcn: Function that evaluates the optional nonlinear constraints
 113                    and their gradients for a given value of M{x}.
 114     @type gh_fcn: callable
 115     @param hess_fcn: Handle to function that computes the Hessian of the
 116                      Lagrangian for given values of M{x}, M{lambda} and M{mu},
 117                      where M{lambda} and M{mu} are the multipliers on the
 118                      equality and inequality constraints, M{g} and M{h},
 119                      respectively.
 120     @type hess_fcn: callable
 121     @param opt: optional options dictionary with the following keys, all of
 122                 which are also optional (default values shown in parentheses)
 123                   - C{verbose} (False) - Controls level of progress output
 124                     displayed
 125                   - C{feastol} (1e-6) - termination tolerance for feasibility
 126                     condition
 127                   - C{gradtol} (1e-6) - termination tolerance for gradient
 128                     condition
 129                   - C{comptol} (1e-6) - termination tolerance for
 130                     complementarity condition
 131                   - C{costtol} (1e-6) - termination tolerance for cost
 132                     condition
 133                   - C{max_it} (150) - maximum number of iterations
 134                   - C{step_control} (False) - set to True to enable step-size
 135                     control
 136                   - C{max_red} (20) - maximum number of step-size reductions if
 137                     step-control is on
 138                   - C{cost_mult} (1.0) - cost multiplier used to scale the
 139                     objective function for improved conditioning. Note: This
 140                     value is also passed as the 3rd argument to the Hessian
 141                     evaluation function so that it can appropriately scale the
 142                     objective function term in the Hessian of the Lagrangian.
 143     @type opt: dict
 144
 145     @rtype: dict
 146     @return: The solution dictionary has the following keys:
 147                - C{x} - solution vector
 148                - C{f} - final objective function value
 149                - C{converged} - exit status
 150                    - True = first order optimality conditions satisfied
 151                    - False = maximum number of iterations reached
 152                    - None = numerically failed
 153                - C{output} - output dictionary with keys:
 154                    - C{iterations} - number of iterations performed
 155                    - C{hist} - list of arrays with trajectories of the
 156                      following: feascond, gradcond, compcond, costcond, gamma,
 157                      stepsize, obj, alphap, alphad
 158                    - C{message} - exit message
 159                - C{lmbda} - dictionary containing the Langrange and Kuhn-Tucker
 160                  multipliers on the constraints, with keys:
 161                    - C{eqnonlin} - nonlinear equality constraints
 162                    - C{ineqnonlin} - nonlinear inequality constraints
 163                    - C{mu_l} - lower (left-hand) limit on linear constraints
 164                    - C{mu_u} - upper (right-hand) limit on linear constraints
 165                    - C{lower} - lower bound on optimization variables
 166                    - C{upper} - upper bound on optimization variables
 167
 168     @see: U{http://www.pserc.cornell.edu/matpower/}
 169
 170     @author: Ray Zimmerman (PSERC Cornell)
 171     """
 172     if isinstance(f_fcn, dict):  ## problem dict
 173         p = f_fcn
 174         f_fcn = p['f_fcn']
 175         x0 = p['x0']
 176         if 'opt' in p: opt = p['opt']
 177         if 'hess_fcn' in p: hess_fcn = p['hess_fcn']
 178         if 'gh_fcn' in p: gh_fcn = p['gh_fcn']
 179         if 'xmax' in p: xmax = p['xmax']
 180         if 'xmin' in p: xmin = p['xmin']
 181         if 'u' in p: u = p['u']
 182         if 'l' in p: l = p['l']
 183         if 'A' in p: A = p['A']
 184
 185     nx = x0.shape[0]                        # number of variables
 186     nA = A.shape[0] if A is not None else 0 # number of original linear constr
 187
 188     # default argument values
 189     if l is None or len(l) == 0: l = -Inf * ones(nA)
 190     if u is None or len(u) == 0: u =  Inf * ones(nA)
 191     if xmin is None or len(xmin) == 0: xmin = -Inf * ones(x0.shape[0])
 192     if xmax is None or len(xmax) == 0: xmax =  Inf * ones(x0.shape[0])
 193     if gh_fcn is None:
 194         nonlinear = False
 195         gn = array([])
 196         hn = array([])
 197     else:
 198         nonlinear = True
 199
 200     if opt is None: opt = {}
 201     # options
 202     if "feastol" not in opt:
 203         opt["feastol"] = 1e-06
 204     if "gradtol" not in opt:
 205         opt["gradtol"] = 1e-06
 206     if "comptol" not in opt:
 207         opt["comptol"] = 1e-06
 208     if "costtol" not in opt:
 209         opt["costtol"] = 1e-06
 210     if "max_it" not in opt:
 211         opt["max_it"] = 150
 212     if "max_red" not in opt:
 213         opt["max_red"] = 20
 214     if "step_control" not in opt:
 215         opt["step_control"] = False
 216     if "cost_mult" not in opt:
 217         opt["cost_mult"] = 1
 218     if "verbose" not in opt:
 219         opt["verbose"] = 0
 220
 221     # initialize history
 222     hist = []
 223
 224     # constants
 225     xi = 0.99995
 226     sigma = 0.1
 227     z0 = 1
 228     alpha_min = 1e-8
 229     rho_min = 0.95
 230     rho_max = 1.05
 231     mu_threshold = 1e-5
 232
 233     # initialize
 234     i = 0                       # iteration counter
 235     converged = False           # flag
 236     eflag = False               # exit flag
 237
 238     # add var limits to linear constraints
 239     eyex = eye(nx, nx, format="csr")
 240     AA = eyex if A is None else vstack([eyex, A], "csr")
 241     ll = r_[xmin, l]
 242     uu = r_[xmax, u]
 243
 244     # split up linear constraints
 245     ieq = find( absolute(uu - ll) <= EPS )
 246     igt = find( (uu >=  1e10) & (ll > -1e10) )
 247     ilt = find( (ll <= -1e10) & (uu <  1e10) )
 248     ibx = find( (absolute(uu - ll) > EPS) & (uu < 1e10) & (ll > -1e10) )
 249     # zero-sized sparse matrices unsupported
 250     Ae = AA[ieq, :] if len(ieq) else None
 251     if len(ilt) or len(igt) or len(ibx):
 252         idxs = [(1, ilt), (-1, igt), (1, ibx), (-1, ibx)]
 253         Ai = vstack([sig * AA[idx, :] for sig, idx in idxs if len(idx)], 'csr')
 254     else:
 255         Ai = None
 256     be = uu[ieq]
 257     bi = r_[uu[ilt], -ll[igt], uu[ibx], -ll[ibx]]
 258
 259     # evaluate cost f(x0) and constraints g(x0), h(x0)
 260     x = x0
 261     f, df = f_fcn(x)                 # cost
 262     f = f * opt["cost_mult"]
 263     df = df * opt["cost_mult"]
 264     if nonlinear:
 265         hn, gn, dhn, dgn = gh_fcn(x)        # nonlinear constraints
 266         h = hn if Ai is None else r_[hn, Ai * x - bi] # inequality constraints
 267         g = gn if Ae is None else r_[gn, Ae * x - be] # equality constraints
 268
 269         if (dhn is None) and (Ai is None):
 270             dh = None
 271         elif dhn is None:
 272             dh = Ai.T
 273         elif Ai is None:
 274             dh = dhn
 275         else:
 276             dh = hstack([dhn, Ai.T])
 277
 278         if (dgn is None) and (Ae is None):
 279             dg = None
 280         elif dgn is None:
 281             dg = Ae.T
 282         elif Ae is None:
 283             dg = dgn
 284         else:
 285             dg = hstack([dgn, Ae.T])
 286     else:
 287         h = -bi if Ai is None else Ai * x - bi        # inequality constraints
 288         g = -be if Ae is None else Ae * x - be        # equality constraints
 289         dh = None if Ai is None else Ai.T     # 1st derivative of inequalities
 290         dg = None if Ae is None else Ae.T     # 1st derivative of equalities
 291
 292     # some dimensions
 293     neq = g.shape[0]           # number of equality constraints
 294     niq = h.shape[0]           # number of inequality constraints
 295     neqnln = gn.shape[0]       # number of nonlinear equality constraints
 296     niqnln = hn.shape[0]       # number of nonlinear inequality constraints
 297     nlt = len(ilt)             # number of upper bounded linear inequalities
 298     ngt = len(igt)             # number of lower bounded linear inequalities
 299     nbx = len(ibx)             # number of doubly bounded linear inequalities
 300
 301     # initialize gamma, lam, mu, z, e
 302     gamma = 1                  # barrier coefficient
 303     lam = zeros(neq)
 304     z = z0 * ones(niq)
 305     mu = z0 * ones(niq)
 306     k = find(h < -z0)
 307     z[k] = -h[k]
 308     k = find((gamma / z) > z0)
 309     mu[k] = gamma / z[k]
 310     e = ones(niq)
 311
 312     # check tolerance
 313     f0 = f
 314     if opt["step_control"]:
 315         L = f + dot(lam, g) + dot(mu, h + z) - gamma * sum(log(z))
 316
 317     Lx = df.copy()
 318     Lx = Lx + dg * lam if dg is not None else Lx
 319     Lx = Lx + dh * mu  if dh is not None else Lx
 320
 321     maxh = zeros(1) if len(h) == 0 else max(h)
 322
 323     gnorm = norm(g, Inf) if len(g) else 0.0
 324     lam_norm = norm(lam, Inf) if len(lam) else 0.0
 325     mu_norm = norm(mu, Inf) if len(mu) else 0.0
 326     znorm = norm(z, Inf) if len(z) else 0.0
 327     feascond = \
 328         max([gnorm, maxh]) / (1 + max([norm(x, Inf), znorm]))
 329     gradcond = \
 330         norm(Lx, Inf) / (1 + max([lam_norm, mu_norm]))
 331     compcond = dot(z, mu) / (1 + norm(x, Inf))
 332     costcond = absolute(f - f0) / (1 + absolute(f0))
 333
 334     # save history
 335     hist.append({'feascond': feascond, 'gradcond': gradcond,
 336         'compcond': compcond, 'costcond': costcond, 'gamma': gamma,
 337         'stepsize': 0, 'obj': f / opt["cost_mult"], 'alphap': 0, 'alphad': 0})
 338
 339     if opt["verbose"]:
 340         s = '-sc' if opt["step_control"] else ''
 341         v = pipsver('all')
 342         print('Python Interior Point Solver - PIPS%s, Version %s, %s' %
 343                     (s, v['Version'], v['Date']))
 344         if opt['verbose'] > 1:
 345             print(" it    objective   step size   feascond     gradcond     "
 346                   "compcond     costcond  ")
 347             print("----  ------------ --------- ------------ ------------ "
 348                   "------------ ------------")
 349             print("%3d  %12.8g %10s %12g %12g %12g %12g" %
 350                 (i, (f / opt["cost_mult"]), "",
 351                  feascond, gradcond, compcond, costcond))
 352
 353     if feascond < opt["feastol"] and gradcond < opt["gradtol"] and \
 354         compcond < opt["comptol"] and costcond < opt["costtol"]:
 355         converged = True
 356         if opt["verbose"]:
 357             print("Converged!")
 358
 359     # do Newton iterations
 360     while (not converged) and (i < opt["max_it"]):
 361         # update iteration counter
 362         i += 1
 363
 364         # compute update step
 365         lmbda = {"eqnonlin": lam[range(neqnln)],
 366                  "ineqnonlin": mu[range(niqnln)]}
 367         if nonlinear:
 368             if hess_fcn is None:
 369                 print("pips: Hessian evaluation via finite differences "
 370                       "not yet implemented.\nPlease provide "
 371                       "your own hessian evaluation function.")
 372             Lxx = hess_fcn(x, lmbda, opt["cost_mult"])
 373         else:
 374             _, _, d2f = f_fcn(x, True)      # cost
 375             Lxx = d2f * opt["cost_mult"]
 376         rz = range(len(z))
 377         zinvdiag = sparse((1.0 / z, (rz, rz))) if len(z) else None
 378         rmu = range(len(mu))
 379         mudiag = sparse((mu, (rmu, rmu))) if len(mu) else None
 380         dh_zinv = None if dh is None else dh * zinvdiag
 381         M = Lxx if dh is None else Lxx + dh_zinv * mudiag * dh.T
 382         N = Lx if dh is None else Lx + dh_zinv * (mudiag * h + gamma * e)
 383
 384         Ab = sparse(M) if dg is None else vstack([
 385             hstack([M, dg]),
 386             hstack([dg.T, sparse((neq, neq))])
 387         ])
 388         bb = r_[-N, -g]
 389
 390         dxdlam = pplinsolve(Ab.tocsr(), bb)
 391
 392         if any(isnan(dxdlam)):
 393             if opt["verbose"]:
 394                 print('\nNumerically Failed\n')
 395             eflag = -1
 396             break
 397
 398         dx = dxdlam[:nx]
 399         dlam = dxdlam[nx:nx + neq]
 400         dz = -h - z if dh is None else -h - z - dh.T * dx
 401         dmu = -mu if dh is None else -mu + zinvdiag * (gamma * e - mudiag * dz)
 402
 403         # optional step-size control
 404         sc = False
 405         if opt["step_control"]:
 406             x1 = x + dx
 407
 408             # evaluate cost, constraints, derivatives at x1
 409             f1, df1 = f_fcn(x1)          # cost
 410             f1 = f1 * opt["cost_mult"]
 411             df1 = df1 * opt["cost_mult"]
 412             if nonlinear:
 413                 hn1, gn1, dhn1, dgn1 = gh_fcn(x1) # nonlinear constraints
 414
 415                 h1 = hn1 if Ai is None else r_[hn1, Ai * x1 - bi] # ieq constraints
 416                 g1 = gn1 if Ae is None else r_[gn1, Ae * x1 - be] # eq constraints
 417
 418                 # 1st der of ieq
 419                 if (dhn1 is None) and (Ai is None):
 420                     dh1 = None
 421                 elif dhn1 is None:
 422                     dh1 = Ai.T
 423                 elif Ai is None:
 424                     dh1 = dhn1
 425                 else:
 426                     dh1 = hstack([dhn1, Ai.T])
 427
 428                 # 1st der of eqs
 429                 if (dgn1 is None) and (Ae is None):
 430                     dg1 = None
 431                 elif dgn is None:
 432                     dg1 = Ae.T
 433                 elif Ae is None:
 434                     dg1 = dgn1
 435                 else:
 436                     dg1 = hstack([dgn1, Ae.T])
 437             else:
 438                 h1 = -bi if Ai is None else Ai * x1 - bi    # inequality constraints
 439                 g1 = -be if Ae is None else Ae * x1 - be    # equality constraints
 440
 441                 dh1 = dh                       ## 1st derivative of inequalities
 442                 dg1 = dg                       ## 1st derivative of equalities
 443
 444             # check tolerance
 445             Lx1 = df1
 446             Lx1 = Lx1 + dg1 * lam if dg1 is not None else Lx1
 447             Lx1 = Lx1 + dh1 * mu  if dh1 is not None else Lx1
 448
 449             maxh1 = zeros(1) if len(h1) == 0 else max(h1)
 450
 451             g1norm = norm(g1, Inf) if len(g1) else 0.0
 452             lam1_norm = norm(lam, Inf) if len(lam) else 0.0
 453             mu1_norm = norm(mu, Inf) if len(mu) else 0.0
 454             z1norm = norm(z, Inf) if len(z) else 0.0
 455
 456             feascond1 = max([ g1norm, maxh1 ]) / \
 457                 (1 + max([ norm(x1, Inf), z1norm ]))
 458             gradcond1 = norm(Lx1, Inf) / (1 + max([ lam1_norm, mu1_norm ]))
 459
 460             if (feascond1 > feascond) and (gradcond1 > gradcond):
 461                 sc = True
 462         if sc:
 463             alpha = 1.0
 464             for j in range(opt["max_red"]):
 465                 dx1 = alpha * dx
 466                 x1 = x + dx1
 467                 f1, _ = f_fcn(x1)             # cost
 468                 f1 = f1 * opt["cost_mult"]
 469                 if nonlinear:
 470                     hn1, gn1, _, _ = gh_fcn(x1)              # nonlinear constraints
 471                     h1 = hn1 if Ai is None else r_[hn1, Ai * x1 - bi]         # inequality constraints
 472                     g1 = gn1 if Ae is None else r_[gn1, Ae * x1 - be]         # equality constraints
 473                 else:
 474                     h1 = -bi if Ai is None else Ai * x1 - bi    # inequality constraints
 475                     g1 = -be if Ae is None else Ae * x1 - be    # equality constraints
 476
 477                 L1 = f1 + dot(lam, g1) + dot(mu, h1 + z) - gamma * sum(log(z))
 478
 479                 if opt["verbose"] > 2:
 480                     print("   %3d            %10.5f" % (-j, norm(dx1)))
 481
 482                 rho = (L1 - L) / (dot(Lx, dx1) + 0.5 * dot(dx1, Lxx * dx1))
 483
 484                 if (rho > rho_min) and (rho < rho_max):
 485                     break
 486                 else:
 487                     alpha = alpha / 2.0
 488             dx = alpha * dx
 489             dz = alpha * dz
 490             dlam = alpha * dlam
 491             dmu = alpha * dmu
 492
 493         # do the update
 494         k = find(dz < 0.0)
 495         alphap = min([xi * min(z[k] / -dz[k]), 1]) if len(k) else 1.0
 496         k = find(dmu < 0.0)
 497         alphad = min([xi * min(mu[k] / -dmu[k]), 1]) if len(k) else 1.0
 498         x = x + alphap * dx
 499         z = z + alphap * dz
 500         lam = lam + alphad * dlam
 501         mu = mu + alphad * dmu
 502         if niq > 0:
 503             gamma = sigma * dot(z, mu) / niq
 504
 505         # evaluate cost, constraints, derivatives
 506         f, df = f_fcn(x)             # cost
 507         f = f * opt["cost_mult"]
 508         df = df * opt["cost_mult"]
 509         if nonlinear:
 510             hn, gn, dhn, dgn = gh_fcn(x)                   # nln constraints
 511 #            g = gn if Ai is None else r_[gn, Ai * x - bi] # ieq constraints
 512 #            h = hn if Ae is None else r_[hn, Ae * x - be] # eq constraints
 513             h = hn if Ai is None else r_[hn, Ai * x - bi] # ieq constr
 514             g = gn if Ae is None else r_[gn, Ae * x - be]  # eq constr
 515
 516             if (dhn is None) and (Ai is None):
 517                 dh = None
 518             elif dhn is None:
 519                 dh = Ai.T
 520             elif Ai is None:
 521                 dh = dhn
 522             else:
 523                 dh = hstack([dhn, Ai.T])
 524
 525             if (dgn is None) and (Ae is None):
 526                 dg = None
 527             elif dgn is None:
 528                 dg = Ae.T
 529             elif Ae is None:
 530                 dg = dgn
 531             else:
 532                 dg = hstack([dgn, Ae.T])
 533         else:
 534             h = -bi if Ai is None else Ai * x - bi    # inequality constraints
 535             g = -be if Ae is None else Ae * x - be    # equality constraints
 536             # 1st derivatives are constant, still dh = Ai.T, dg = Ae.T
 537
 538         Lx = df
 539         Lx = Lx + dg * lam if dg is not None else Lx
 540         Lx = Lx + dh * mu  if dh is not None else Lx
 541
 542         if len(h) == 0:
 543             maxh = zeros(1)
 544         else:
 545             maxh = max(h)
 546
 547         gnorm = norm(g, Inf) if len(g) else 0.0
 548         lam_norm = norm(lam, Inf) if len(lam) else 0.0
 549         mu_norm = norm(mu, Inf) if len(mu) else 0.0
 550         znorm = norm(z, Inf) if len(z) else 0.0
 551         feascond = \
 552             max([gnorm, maxh]) / (1 + max([norm(x, Inf), znorm]))
 553         gradcond = \
 554             norm(Lx, Inf) / (1 + max([lam_norm, mu_norm]))
 555         compcond = dot(z, mu) / (1 + norm(x, Inf))
 556         costcond = float(absolute(f - f0) / (1 + absolute(f0)))
 557
 558         hist.append({'feascond': feascond, 'gradcond': gradcond,
 559             'compcond': compcond, 'costcond': costcond, 'gamma': gamma,
 560             'stepsize': norm(dx), 'obj': f / opt["cost_mult"],
 561             'alphap': alphap, 'alphad': alphad})
 562
 563         if opt["verbose"] > 1:
 564             print("%3d  %12.8g %10.5g %12g %12g %12g %12g" %
 565                 (i, (f / opt["cost_mult"]), norm(dx), feascond, gradcond,
 566                  compcond, costcond))
 567
 568         if feascond < opt["feastol"] and gradcond < opt["gradtol"] and \
 569             compcond < opt["comptol"] and costcond < opt["costtol"]:
 570             converged = True
 571             if opt["verbose"]:
 572                 print("Converged!")
 573         else:
 574             if any(isnan(x)) or (alphap < alpha_min) or \
 575                 (alphad < alpha_min) or (gamma < EPS) or (gamma > 1.0 / EPS):
 576                 if opt["verbose"]:
 577                     print("Numerically failed.")
 578                 eflag = -1
 579                 break
 580             f0 = f
 581
 582             if opt["step_control"]:
 583                 L = f + dot(lam, g) + dot(mu, (h + z)) - gamma * sum(log(z))
 584
 585     if opt["verbose"]:
 586         if not converged:
 587             print("Did not converge in %d iterations." % i)
 588
 589     # package results
 590     if eflag != -1:
 591         eflag = converged
 592
 593     if eflag == 0:
 594         message = 'Did not converge'
 595     elif eflag == 1:
 596         message = 'Converged'
 597     elif eflag == -1:
 598         message = 'Numerically failed'
 599     else:
 600         raise
 601
 602     output = {"iterations": i, "hist": hist, "message": message}
 603
 604     # zero out multipliers on non-binding constraints
 605     mu[find( (h < -opt["feastol"]) & (mu < mu_threshold) )] = 0.0
 606
 607     # un-scale cost and prices
 608     f = f / opt["cost_mult"]
 609     lam = lam / opt["cost_mult"]
 610     mu = mu / opt["cost_mult"]
 611
 612     # re-package multipliers into struct
 613     lam_lin = lam[neqnln:neq]           # lambda for linear constraints
 614     mu_lin = mu[niqnln:niq]             # mu for linear constraints
 615     kl = find(lam_lin < 0.0)     # lower bound binding
 616     ku = find(lam_lin > 0.0)     # upper bound binding
 617
 618     mu_l = zeros(nx + nA)
 619     mu_l[ieq[kl]] = -lam_lin[kl]
 620     mu_l[igt] = mu_lin[nlt:nlt + ngt]
 621     mu_l[ibx] = mu_lin[nlt + ngt + nbx:nlt + ngt + nbx + nbx]
 622
 623     mu_u = zeros(nx + nA)
 624     mu_u[ieq[ku]] = lam_lin[ku]
 625     mu_u[ilt] = mu_lin[:nlt]
 626     mu_u[ibx] = mu_lin[nlt + ngt:nlt + ngt + nbx]
 627
 628     lmbda = {'mu_l': mu_l[nx:], 'mu_u': mu_u[nx:],
 629              'lower': mu_l[:nx], 'upper': mu_u[:nx]}
 630
 631     if niqnln > 0:
 632         lmbda['ineqnonlin'] = mu[:niqnln]
 633     if neqnln > 0:
 634         lmbda['eqnonlin'] = lam[:neqnln]
 635
 636 #    lmbda = {"eqnonlin": lam[:neqnln], 'ineqnonlin': mu[:niqnln],
 637 #             "mu_l": mu_l[nx:], "mu_u": mu_u[nx:],
 638 #             "lower": mu_l[:nx], "upper": mu_u[:nx]}
 639
 640     solution =  {"x": x, "f": f, "eflag": converged,
 641                  "output": output, "lmbda": lmbda}
 642
 643     return solution
 644
 645 if __name__ == "__main__":
 646     import doctest
 647     doctest.testmod()