pypower.pips

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

Source Code for Module pypower.pips