pypower.ipoptopf

Source Code for Module pypower.ipoptopf_solver

1 # Copyright (C) 2000-2011 Power System Engineering Research Center (PSERC) 2 # Copyright (C) 2011 Richard Lincoln 3 # 4 # PYPOWER is free software: you can redistribute it and/or modify 5 # it under the terms of the GNU General Public License as published 6 # by the Free Software Foundation, either version 3 of the License, 7 # or (at your option) any later version. 8 # 9 # PYPOWER is distributed in the hope that it will be useful, 10 # but WITHOUT ANY WARRANTY], without even the implied warranty of 11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 # GNU General Public License for more details. 13 # 14 # You should have received a copy of the GNU General Public License 15 # along with PYPOWER. If not, see <http://www.gnu.org/licenses/>. 16 17 """Solves AC optimal power flow using IPOPT. 18 """ 19 20 from numpy import ones, zeros, shape, Inf, pi, exp, conj, r_, arange 21 from numpy import flatnonzero as find 22 23 from scipy.sparse import issparse, tril, vstack, hstack, csr_matrix as sparse 24 from scipy.sparse import eye as speye 25 26 from pypower.idx_bus import BUS_TYPE, REF, VM, VA, MU_VMAX, MU_VMIN, LAM_P, LAM_Q 27 from pypower.idx_brch import F_BUS, T_BUS, RATE_A, PF, QF, PT, QT, MU_SF, MU_ST 28 from pypower.idx_gen import GEN_BUS, PG, QG, VG, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN 29 from pypower.idx_cost import MODEL, PW_LINEAR, NCOST 30 31 from pypower.makeYbus import makeYbus 32 from pypower.opf_costfcn import opf_costfcn 33 from pypower.opf_consfcn import opf_consfcn 34 from pypower.opf_hessfcn import opf_hessfcn 35 from pypower.util import sub2ind 36 from pypower.ipopt_options import ipopt_options 37 38

39 -def ipoptopf_solver(om, ppopt):

40 """Solves AC optimal power flow using IPOPT. 41 42 Inputs are an OPF model object and a PYPOWER options vector. 43 44 Outputs are a C{results} dict, C{success} flag and C{raw} output dict. 45 46 C{results} is a PYPOWER case dict (ppc) with the usual C{baseMVA}, C{bus} 47 C{branch}, C{gen}, C{gencost} fields, along with the following additional 48 fields: 49 - C{order} see 'help ext2int' for details of this field 50 - C{x} final value of optimization variables (internal order) 51 - C{f} final objective function value 52 - C{mu} shadow prices on ... 53 - C{var} 54 - C{l} lower bounds on variables 55 - C{u} upper bounds on variables 56 - C{nln} 57 - C{l} lower bounds on nonlinear constraints 58 - C{u} upper bounds on nonlinear constraints 59 - C{lin} 60 - C{l} lower bounds on linear constraints 61 - C{u} upper bounds on linear constraints 62 63 C{success} is C{True} if solver converged successfully, C{False} otherwise 64 65 C{raw} is a raw output dict in form returned by MINOS 66 - C{xr} final value of optimization variables 67 - C{pimul} constraint multipliers 68 - C{info} solver specific termination code 69 - C{output} solver specific output information 70 71 @see: L{opf}, L{pips} 72 73 @author: Ray Zimmerman (PSERC Cornell) 74 @author: Carlos E. Murillo-Sanchez (PSERC Cornell & Universidad 75 Autonoma de Manizales) 76 @author: Richard Lincoln 77 """ 78 import pyipopt 79 80 ## options 81 verbose = ppopt['VERBOSE'] 82 feastol = ppopt['PDIPM_FEASTOL'] 83 gradtol = ppopt['PDIPM_GRADTOL'] 84 comptol = ppopt['PDIPM_COMPTOL'] 85 costtol = ppopt['PDIPM_COSTTOL'] 86 max_it = ppopt['PDIPM_MAX_IT'] 87 max_red = ppopt['SCPDIPM_RED_IT'] 88 step_control = ppopt['OPF_ALG'] == 565 # PIPS-sc 89 if feastol == 0: 90 feastol = ppopt['OPF_VIOLATION'] ## = OPF_VIOLATION by default 91 92 opt = { 'feastol': feastol, \ 93 'gradtol': gradtol, \ 94 'comptol': comptol, \ 95 'costtol': costtol, \ 96 'max_it': max_it, \ 97 'max_red': max_red, \ 98 'step_control': step_control, \ 99 'cost_mult': 1e-4, \ 100 'verbose': verbose } 101 102 ## unpack data 103 ppc = om.get_ppc() 104 baseMVA, bus, gen, branch, gencost = \ 105 ppc['baseMVA'], ppc['bus'], ppc['gen'], ppc['branch'], ppc['gencost'] 106 vv, _, nn, _ = om.get_idx() 107 108 ## problem dimensions 109 nb = shape(bus)[0] ## number of buses 110 ng = shape(gen)[0] ## number of buses 111 nl = shape(branch)[0] ## number of branches 112 ny = om.getN('var', 'y') ## number of piece-wise linear costs 113 114 ## linear constraints 115 A, l, u = om.linear_constraints() 116 117 ## bounds on optimization vars 118 _, xmin, xmax = om.getv() 119 120 ## build admittance matrices 121 Ybus, Yf, Yt = makeYbus(baseMVA, bus, branch) 122 123 ## try to select an interior initial point 124 ll = xmin.copy(); uu = xmax.copy() 125 ll[xmin == -Inf] = -2e19 ## replace Inf with numerical proxies 126 uu[xmax == Inf] = 2e19 127 x0 = (ll + uu) / 2 128 Varefs = bus[bus[:, BUS_TYPE] == REF, VA] * (pi / 180) 129 x0[vv['i1']['Va']:vv['iN']['Va']] = Varefs[0] ## angles set to first reference angle 130 if ny > 0: 131 ipwl = find(gencost[:, MODEL] == PW_LINEAR) 132 # PQ = r_[gen[:, PMAX], gen[:, QMAX]] 133 # c = totcost(gencost[ipwl, :], PQ[ipwl]) 134 ## largest y-value in CCV data 135 c = gencost.flatten('F')[sub2ind(shape(gencost), ipwl, NCOST + 2 * gencost[ipwl, NCOST])] 136 x0[vv['i1']['y']:vv['iN']['y']] = max(c) + 0.1 * abs(max(c)) 137 # x0[vv['i1']['y']:vv['iN']['y']) = c + 0.1 * abs(c) 138 139 ## find branches with flow limits 140 il = find((branch[:, RATE_A] != 0) & (branch[:, RATE_A] < 1e10)) 141 nl2 = len(il) ## number of constrained lines 142 143 ##----- run opf ----- 144 ## build Jacobian and Hessian structure 145 if A is not None and issparse(A): 146 nA = A.shape[0] ## number of original linear constraints 147 else: 148 nA = 0 149 nx = len(x0) 150 f = branch[:, F_BUS] ## list of "from" buses 151 t = branch[:, T_BUS] ## list of "to" buses 152 Cf = sparse((ones(nl), (arange(nl), f)), (nl, nb)) ## connection matrix for line & from buses 153 Ct = sparse((ones(nl), (arange(nl), t)), (nl, nb)) ## connection matrix for line & to buses 154 Cl = Cf + Ct 155 Cb = Cl.T * Cl + speye(nb, nb) 156 Cl2 = Cl[il, :] 157 Cg = sparse((ones(ng), (gen[:, GEN_BUS], arange(ng))), (nb, ng)) 158 nz = nx - 2 * (nb + ng) 159 nxtra = nx - 2 * nb 160 if nz > 0: 161 Js = vstack([ 162 hstack([Cb, Cb, Cg, sparse((nb, ng)), sparse((nb, nz))]), 163 hstack([Cb, Cb, sparse((nb, ng)), Cg, sparse((nb, nz))]), 164 hstack([Cl2, Cl2, sparse((nl2, 2 * ng)), sparse((nl2, nz))]), 165 hstack([Cl2, Cl2, sparse((nl2, 2 * ng)), sparse((nl2, nz))]) 166 ], 'coo') 167 else: 168 Js = vstack([ 169 hstack([Cb, Cb, Cg, sparse((nb, ng))]), 170 hstack([Cb, Cb, sparse((nb, ng)), Cg, ]), 171 hstack([Cl2, Cl2, sparse((nl2, 2 * ng)), ]), 172 hstack([Cl2, Cl2, sparse((nl2, 2 * ng)), ]) 173 ], 'coo') 174 175 if A is not None and issparse(A): 176 Js = vstack([Js, A], 'coo') 177 178 f, _, d2f = opf_costfcn(x0, om, True) 179 Hs = tril(d2f + vstack([ 180 hstack([Cb, Cb, sparse((nb, nxtra))]), 181 hstack([Cb, Cb, sparse((nb, nxtra))]), 182 sparse((nxtra, nx)) 183 ]), format='coo') 184 185 ## set options struct for IPOPT 186 # options = {} 187 # options['ipopt'] = ipopt_options([], ppopt) 188 189 ## extra data to pass to functions 190 userdata = { 191 'om': om, 192 'Ybus': Ybus, 193 'Yf': Yf[il, :], 194 'Yt': Yt[il, :], 195 'ppopt': ppopt, 196 'il': il, 197 'A': A, 198 'nA': nA, 199 'neqnln': 2 * nb, 200 'niqnln': 2 * nl2, 201 'Js': Js, 202 'Hs': Hs 203 } 204 205 ## check Jacobian and Hessian structure 206 #xr = rand(x0.shape) 207 #lmbda = rand( 2 * nb + 2 * nl2) 208 #Js1 = eval_jac_g(x, flag, userdata) #(xr, options.auxdata) 209 #Hs1 = eval_h(xr, 1, lmbda, userdata) 210 #i1, j1, s = find(Js) 211 #i2, j2, s = find(Js1) 212 #if (len(i1) != len(i2)) | (norm(i1 - i2) != 0) | (norm(j1 - j2) != 0): 213 # raise ValueError, 'something''s wrong with the Jacobian structure' 214 # 215 #i1, j1, s = find(Hs) 216 #i2, j2, s = find(Hs1) 217 #if (len(i1) != len(i2)) | (norm(i1 - i2) != 0) | (norm(j1 - j2) != 0): 218 # raise ValueError, 'something''s wrong with the Hessian structure' 219 220 ## define variable and constraint bounds 221 # n is the number of variables 222 n = x0.shape[0] 223 # xl is the lower bound of x as bounded constraints 224 xl = xmin 225 # xu is the upper bound of x as bounded constraints 226 xu = xmax 227 228 neqnln = 2 * nb 229 niqnln = 2 * nl2 230 231 # number of constraints 232 m = neqnln + niqnln + nA 233 # lower bound of constraint 234 gl = r_[zeros(neqnln), -Inf * ones(niqnln), l] 235 # upper bound of constraints 236 gu = r_[zeros(neqnln), zeros(niqnln), u] 237 238 # number of nonzeros in Jacobi matrix 239 nnzj = Js.nnz 240 # number of non-zeros in Hessian matrix, you can set it to 0 241 nnzh = Hs.nnz 242 243 eval_hessian = True 244 if eval_hessian: 245 hessian = lambda x, lagrange, obj_factor, flag, user_data=None: \ 246 eval_h(x, lagrange, obj_factor, flag, userdata) 247 248 nlp = pyipopt.create(n, xl, xu, m, gl, gu, nnzj, nnzh, 249 eval_f, eval_grad_f, eval_g, eval_jac_g, hessian) 250 else: 251 nnzh = 0 252 nlp = pyipopt.create(n, xl, xu, m, gl, gu, nnzj, nnzh, 253 eval_f, eval_grad_f, eval_g, eval_jac_g) 254 255 nlp.int_option('print_level', 5) 256 nlp.num_option('tol', 1.0000e-12) 257 nlp.int_option('max_iter', 250) 258 nlp.num_option('dual_inf_tol', 0.10000) 259 nlp.num_option('constr_viol_tol', 1.0000e-06) 260 nlp.num_option('compl_inf_tol', 1.0000e-05) 261 nlp.num_option('acceptable_tol', 1.0000e-08) 262 nlp.num_option('acceptable_constr_viol_tol', 1.0000e-04) 263 nlp.num_option('acceptable_compl_inf_tol', 0.0010000) 264 nlp.str_option('mu_strategy', 'adaptive') 265 266 iter = 0 267 def intermediate_callback(algmod, iter_count, obj_value, inf_pr, inf_du, 268 mu, d_norm, regularization_size, alpha_du, alpha_pr, ls_trials, 269 user_data=None): 270 iter = iter_count 271 return True

272 273 nlp.set_intermediate_callback(intermediate_callback) 274 275 ## run the optimization 276 # returns final solution x, upper and lower bound for multiplier, final 277 # objective function obj and the return status of ipopt 278 x, zl, zu, obj, status, zg = nlp.solve(x0, m, userdata) 279 280 info = {'x': x, 'zl': zl, 'zu': zu, 'obj': obj, 'status': status, 'lmbda': zg} 281 282 nlp.close() 283 284 success = (status == 0) | (status == 1) 285 286 output = {'iterations': iter} 287 288 f, _ = opf_costfcn(x, om) 289 290 ## update solution data 291 Va = x[vv['i1']['Va']:vv['iN']['Va']] 292 Vm = x[vv['i1']['Vm']:vv['iN']['Vm']] 293 Pg = x[vv['i1']['Pg']:vv['iN']['Pg']] 294 Qg = x[vv['i1']['Qg']:vv['iN']['Qg']] 295 V = Vm * exp(1j * Va) 296 297 ##----- calculate return values ----- 298 ## update voltages & generator outputs 299 bus[:, VA] = Va * 180 / pi 300 bus[:, VM] = Vm 301 gen[:, PG] = Pg * baseMVA 302 gen[:, QG] = Qg * baseMVA 303 gen[:, VG] = Vm[gen[:, GEN_BUS].astype(int)] 304 305 ## compute branch flows 306 f_br = branch[:, F_BUS].astype(int) 307 t_br = branch[:, T_BUS].astype(int) 308 Sf = V[f_br] * conj(Yf * V) ## cplx pwr at "from" bus, p.u. 309 St = V[t_br] * conj(Yt * V) ## cplx pwr at "to" bus, p.u. 310 branch[:, PF] = Sf.real * baseMVA 311 branch[:, QF] = Sf.imag * baseMVA 312 branch[:, PT] = St.real * baseMVA 313 branch[:, QT] = St.imag * baseMVA 314 315 ## line constraint is actually on square of limit 316 ## so we must fix multipliers 317 muSf = zeros(nl) 318 muSt = zeros(nl) 319 if len(il) > 0: 320 muSf[il] = 2 * info['lmbda'][2 * nb + arange(nl2)] * branch[il, RATE_A] / baseMVA 321 muSt[il] = 2 * info['lmbda'][2 * nb + nl2 + arange(nl2)] * branch[il, RATE_A] / baseMVA 322 323 ## update Lagrange multipliers 324 bus[:, MU_VMAX] = info['zu'][vv['i1']['Vm']:vv['iN']['Vm']] 325 bus[:, MU_VMIN] = info['zl'][vv['i1']['Vm']:vv['iN']['Vm']] 326 gen[:, MU_PMAX] = info['zu'][vv['i1']['Pg']:vv['iN']['Pg']] / baseMVA 327 gen[:, MU_PMIN] = info['zl'][vv['i1']['Pg']:vv['iN']['Pg']] / baseMVA 328 gen[:, MU_QMAX] = info['zu'][vv['i1']['Qg']:vv['iN']['Qg']] / baseMVA 329 gen[:, MU_QMIN] = info['zl'][vv['i1']['Qg']:vv['iN']['Qg']] / baseMVA 330 bus[:, LAM_P] = info['lmbda'][nn['i1']['Pmis']:nn['iN']['Pmis']] / baseMVA 331 bus[:, LAM_Q] = info['lmbda'][nn['i1']['Qmis']:nn['iN']['Qmis']] / baseMVA 332 branch[:, MU_SF] = muSf / baseMVA 333 branch[:, MU_ST] = muSt / baseMVA 334 335 ## package up results 336 nlnN = om.getN('nln') 337 338 ## extract multipliers for nonlinear constraints 339 kl = find(info['lmbda'][:2 * nb] < 0) 340 ku = find(info['lmbda'][:2 * nb] > 0) 341 nl_mu_l = zeros(nlnN) 342 nl_mu_u = r_[zeros(2 * nb), muSf, muSt] 343 nl_mu_l[kl] = -info['lmbda'][kl] 344 nl_mu_u[ku] = info['lmbda'][ku] 345 346 ## extract multipliers for linear constraints 347 lam_lin = info['lmbda'][2 * nb + 2 * nl2 + arange(nA)] ## lmbda for linear constraints 348 kl = find(lam_lin < 0) ## lower bound binding 349 ku = find(lam_lin > 0) ## upper bound binding 350 mu_l = zeros(nA) 351 mu_l[kl] = -lam_lin[kl] 352 mu_u = zeros(nA) 353 mu_u[ku] = lam_lin[ku] 354 355 mu = { 356 'var': {'l': info['zl'], 'u': info['zu']}, 357 'nln': {'l': nl_mu_l, 'u': nl_mu_u}, \ 358 'lin': {'l': mu_l, 'u': mu_u} 359 } 360 361 results = ppc 362 results['bus'], results['branch'], results['gen'], \ 363 results['om'], results['x'], results['mu'], results['f'] = \ 364 bus, branch, gen, om, x, mu, f 365 366 pimul = r_[ 367 results['mu']['nln']['l'] - results['mu']['nln']['u'], 368 results['mu']['lin']['l'] - results['mu']['lin']['u'], 369 -ones(ny > 0), 370 results['mu']['var']['l'] - results['mu']['var']['u'] 371 ] 372 raw = {'xr': x, 'pimul': pimul, 'info': info['status'], 'output': output} 373 374 return results, success, raw 375 376

377 -def eval_f(x, user_data=None):

378 """Calculates the objective value. 379 380 @param x: input vector 381 """ 382 om = user_data['om'] 383 f, _ = opf_costfcn(x, om) 384 return f

385 386

387 -def eval_grad_f(x, user_data=None):

388 """Calculates gradient for objective function. 389 """ 390 om = user_data['om'] 391 _, df = opf_costfcn(x, om) 392 return df

393 394

395 -def eval_g(x, user_data=None):

396 """Calculates the constraint values and returns an array. 397 """ 398 om = user_data['om'] 399 Ybus = user_data['Ybus'] 400 Yf = user_data['Yf'] 401 Yt = user_data['Yt'] 402 ppopt = user_data['ppopt'] 403 il = user_data['il'] 404 A = user_data['A'] 405 406 hn, gn, _, _ = opf_consfcn(x, om, Ybus, Yf, Yt, ppopt, il) 407 408 if A is not None and issparse(A): 409 c = r_[gn, hn, A * x] 410 else: 411 c = r_[gn, hn] 412 return c

413 414

415 -def eval_jac_g(x, flag, user_data=None):

416 """Calculates the Jacobi matrix. 417 418 If the flag is true, returns a tuple (row, col) to indicate the 419 sparse Jacobi matrix's structure. 420 If the flag is false, returns the values of the Jacobi matrix 421 with length nnzj. 422 """ 423 Js = user_data['Js'] 424 if flag: 425 return (Js.row, Js.col) 426 else: 427 om = user_data['om'] 428 Ybus = user_data['Ybus'] 429 Yf = user_data['Yf'] 430 Yt = user_data['Yt'] 431 ppopt = user_data['ppopt'] 432 il = user_data['il'] 433 A = user_data['A'] 434 435 _, _, dhn, dgn = opf_consfcn(x, om, Ybus, Yf, Yt, ppopt, il) 436 437 if A is not None and issparse(A): 438 J = vstack([dgn.T, dhn.T, A], 'coo') 439 else: 440 J = vstack([dgn.T, dhn.T], 'coo') 441 442 ## FIXME: Extend PyIPOPT to handle changes in sparsity structure 443 nnzj = Js.nnz 444 Jd = zeros(nnzj) 445 Jc = J.tocsc() 446 for i in xrange(nnzj): 447 Jd[i] = Jc[Js.row[i], Js.col[i]] 448 449 return Jd

450 451

452 -def eval_h(x, lagrange, obj_factor, flag, user_data=None):

453 """Calculates the Hessian matrix (optional). 454 455 If omitted, set nnzh to 0 and Ipopt will use approximated Hessian 456 which will make the convergence slower. 457 """ 458 Hs = user_data['Hs'] 459 if flag: 460 return (Hs.row, Hs.col) 461 else: 462 neqnln = user_data['neqnln'] 463 niqnln = user_data['niqnln'] 464 om = user_data['om'] 465 Ybus = user_data['Ybus'] 466 Yf = user_data['Yf'] 467 Yt = user_data['Yt'] 468 ppopt = user_data['ppopt'] 469 il = user_data['il'] 470 471 lam = {} 472 lam['eqnonlin'] = lagrange[:neqnln] 473 lam['ineqnonlin'] = lagrange[arange(niqnln) + neqnln] 474 475 H = opf_hessfcn(x, lam, om, Ybus, Yf, Yt, ppopt, il, obj_factor) 476 477 Hl = tril(H, format='csc') 478 479 ## FIXME: Extend PyIPOPT to handle changes in sparsity structure 480 nnzh = Hs.nnz 481 Hd = zeros(nnzh) 482 for i in xrange(nnzh): 483 Hd[i] = Hl[Hs.row[i], Hs.col[i]] 484 485 return Hd

486