pypower.toggle

1 # Copyright (C) 2009-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 """Enable or disable fixed reserve requirements. 18 """ 19 20 from sys import stderr 21 22 from pprint import pprint 23 24 from numpy import zeros, ones, arange, Inf, any, flatnonzero as find 25 26 from scipy.sparse import eye as speye 27 from scipy.sparse import csr_matrix as sparse 28 from scipy.sparse import hstack 29 30 from pypower.add_userfcn import add_userfcn 31 from pypower.remove_userfcn import remove_userfcn 32 from pypower.ext2int import ext2int 33 from pypower.int2ext import int2ext 34 from pypower.idx_gen import RAMP_10, PMAX, GEN_STATUS, GEN_BUS 35 36

37 -def toggle_reserves(ppc, on_off):

38 """Enable or disable fixed reserve requirements. 39 40 Enables or disables a set of OPF userfcn callbacks to implement 41 co-optimization of reserves with fixed zonal reserve requirements. 42 43 These callbacks expect to find a 'reserves' field in the input C{ppc}, 44 where C{ppc['reserves']} is a dict with the following fields: 45 - C{zones} C{nrz x ng}, C{zone(i, j) = 1}, if gen C{j} belongs 46 to zone C{i} 0, otherwise 47 - C{req} C{nrz x 1}, zonal reserve requirement in MW 48 - C{cost} (C{ng} or C{ngr}) C{x 1}, cost of reserves in $/MW 49 - C{qty} (C{ng} or C{ngr}) C{x 1}, max quantity of reserves 50 in MW (optional) 51 where C{nrz} is the number of reserve zones and C{ngr} is the number of 52 generators belonging to at least one reserve zone and C{ng} is the total 53 number of generators. 54 55 The 'int2ext' callback also packages up results and stores them in 56 the following output fields of C{results['reserves']}: 57 - C{R} - C{ng x 1}, reserves provided by each gen in MW 58 - C{Rmin} - C{ng x 1}, lower limit on reserves provided by 59 each gen, (MW) 60 - C{Rmax} - C{ng x 1}, upper limit on reserves provided by 61 each gen, (MW) 62 - C{mu.l} - C{ng x 1}, shadow price on reserve lower limit, ($/MW) 63 - C{mu.u} - C{ng x 1}, shadow price on reserve upper limit, ($/MW) 64 - C{mu.Pmax} - C{ng x 1}, shadow price on C{Pg + R <= Pmax} 65 constraint, ($/MW) 66 - C{prc} - C{ng x 1}, reserve price for each gen equal to 67 maximum of the shadow prices on the zonal requirement constraint 68 for each zone the generator belongs to 69 70 @see: L{runopf_w_res}, L{add_userfcn}, L{remove_userfcn}, L{run_userfcn}, 71 L{t.t_case30_userfcns} 72 73 @author: Ray Zimmerman (PSERC Cornell) 74 @author: Richard Lincoln 75 """ 76 if on_off == 'on': 77 ## check for proper reserve inputs 78 if ('reserves' not in ppc) | (not isinstance(ppc['reserves'], dict)) | \ 79 ('zones' not in ppc['reserves']) | \ 80 ('req' not in ppc['reserves']) | \ 81 ('cost' not in ppc['reserves']): 82 stderr.write('toggle_reserves: case must contain a \'reserves\' field, a struct defining \'zones\', \'req\' and \'cost\'\n') 83 84 ## add callback functions 85 ## note: assumes all necessary data included in 1st arg (ppc, om, results) 86 ## so, no additional explicit args are needed 87 ppc = add_userfcn(ppc, 'ext2int', userfcn_reserves_ext2int) 88 ppc = add_userfcn(ppc, 'formulation', userfcn_reserves_formulation) 89 ppc = add_userfcn(ppc, 'int2ext', userfcn_reserves_int2ext) 90 ppc = add_userfcn(ppc, 'printpf', userfcn_reserves_printpf) 91 ppc = add_userfcn(ppc, 'savecase', userfcn_reserves_savecase) 92 elif on_off == 'off': 93 ppc = remove_userfcn(ppc, 'savecase', userfcn_reserves_savecase) 94 ppc = remove_userfcn(ppc, 'printpf', userfcn_reserves_printpf) 95 ppc = remove_userfcn(ppc, 'int2ext', userfcn_reserves_int2ext) 96 ppc = remove_userfcn(ppc, 'formulation', userfcn_reserves_formulation) 97 ppc = remove_userfcn(ppc, 'ext2int', userfcn_reserves_ext2int) 98 else: 99 stderr.write('toggle_reserves: 2nd argument must be either ''on'' or ''off''') 100 101 return ppc

102 103

104 -def userfcn_reserves_ext2int(ppc, *args):

105 """This is the 'ext2int' stage userfcn callback that prepares the input 106 data for the formulation stage. It expects to find a 'reserves' field 107 in ppc as described above. The optional args are not currently used. 108 """ 109 ## initialize some things 110 r = ppc['reserves'] 111 o = ppc['order'] 112 ng0 = o['ext']['gen'].shape[0] ## number of original gens (+ disp loads) 113 nrz = r['req'].shape[0] ## number of reserve zones 114 if nrz > 1: 115 ppc['reserves']['rgens'] = any(r['zones'], 0) ## mask of gens available to provide reserves 116 else: 117 ppc['reserves']['rgens'] = r['zones'] 118 119 igr = find(ppc['reserves']['rgens']) ## indices of gens available to provide reserves 120 ngr = len(igr) ## number of gens available to provide reserves 121 122 ## check data for consistent dimensions 123 if r['zones'].shape[0] != nrz: 124 stderr.write('userfcn_reserves_ext2int: the number of rows in ppc[\'reserves\'][\'req\'] (%d) and ppc[\'reserves\'][\'zones\'] (%d) must match\n' % (nrz, r['zones'].shape[0])) 125 126 if (r['cost'].shape[0] != ng0) & (r['cost'].shape[0] != ngr): 127 stderr.write('userfcn_reserves_ext2int: the number of rows in ppc[\'reserves\'][\'cost\'] (%d) must equal the total number of generators (%d) or the number of generators able to provide reserves (%d)\n' % (r['cost'].shape[0], ng0, ngr)) 128 129 if 'qty' in r: 130 if r['qty'].shape[0] != r['cost'].shape[0]: 131 stderr.write('userfcn_reserves_ext2int: ppc[\'reserves\'][\'cost\'] (%d x 1) and ppc[\'reserves\'][\'qty\'] (%d x 1) must be the same dimension\n' % (r['cost'].shape[0], r['qty'].shape[0])) 132 133 134 ## convert both cost and qty from ngr x 1 to full ng x 1 vectors if necessary 135 if r['cost'].shape[0] < ng0: 136 if 'original' not in ppc['reserves']: 137 ppc['reserves']['original'] = {} 138 ppc['reserves']['original']['cost'] = r['cost'].copy() ## save original 139 cost = zeros(ng0) 140 cost[igr] = r['cost'] 141 ppc['reserves']['cost'] = cost 142 if 'qty' in r: 143 ppc['reserves']['original']['qty'] = r['qty'].copy() ## save original 144 qty = zeros(ng0) 145 qty[igr] = r['qty'] 146 ppc['reserves']['qty'] = qty 147 148 ##----- convert stuff to internal indexing ----- 149 ## convert all reserve parameters (zones, costs, qty, rgens) 150 if 'qty' in r: 151 ppc = ext2int(ppc, ['reserves', 'qty'], 'gen') 152 153 ppc = ext2int(ppc, ['reserves', 'cost'], 'gen') 154 ppc = ext2int(ppc, ['reserves', 'zones'], 'gen', 1) 155 ppc = ext2int(ppc, ['reserves', 'rgens'], 'gen', 1) 156 157 ## save indices of gens available to provide reserves 158 ppc['order']['ext']['reserves']['igr'] = igr ## external indexing 159 ppc['reserves']['igr'] = find(ppc['reserves']['rgens']) ## internal indexing 160 161 return ppc

162 163

164 -def userfcn_reserves_formulation(om, *args):

165 """This is the 'formulation' stage userfcn callback that defines the 166 user costs and constraints for fixed reserves. It expects to find 167 a 'reserves' field in the ppc stored in om, as described above. 168 By the time it is passed to this callback, ppc['reserves'] should 169 have two additional fields: 170 - C{igr} C{1 x ngr}, indices of generators available for reserves 171 - C{rgens} C{1 x ng}, 1 if gen avaiable for reserves, 0 otherwise 172 It is also assumed that if cost or qty were C{ngr x 1}, they have been 173 expanded to C{ng x 1} and that everything has been converted to 174 internal indexing, i.e. all gens are on-line (by the 'ext2int' 175 callback). The optional args are not currently used. 176 """ 177 ## initialize some things 178 ppc = om.get_ppc() 179 r = ppc['reserves'] 180 igr = r['igr'] ## indices of gens available to provide reserves 181 ngr = len(igr) ## number of gens available to provide reserves 182 ng = ppc['gen'].shape[0] ## number of on-line gens (+ disp loads) 183 184 ## variable bounds 185 Rmin = zeros(ngr) ## bound below by 0 186 Rmax = Inf * ones(ngr) ## bound above by ... 187 k = find(ppc['gen'][igr, RAMP_10]) 188 Rmax[k] = ppc['gen'][igr[k], RAMP_10] ## ... ramp rate and ... 189 if 'qty' in r: 190 k = find(r['qty'][igr] < Rmax) 191 Rmax[k] = r['qty'][igr[k]] ## ... stated max reserve qty 192 Rmax = Rmax / ppc['baseMVA'] 193 194 ## constraints 195 I = speye(ngr, ngr, format='csr') ## identity matrix 196 Ar = hstack([sparse((ones(ngr), (arange(ngr), igr)), (ngr, ng)), I], 'csr') 197 ur = ppc['gen'][igr, PMAX] / ppc['baseMVA'] 198 lreq = r['req'] / ppc['baseMVA'] 199 200 ## cost 201 Cw = r['cost'][igr] * ppc['baseMVA'] ## per unit cost coefficients 202 203 ## add them to the model 204 om.add_vars('R', ngr, [], Rmin, Rmax) 205 om.add_constraints('Pg_plus_R', Ar, [], ur, ['Pg', 'R']) 206 om.add_constraints('Rreq', sparse( r['zones'][:, igr] ), lreq, [], ['R']) 207 om.add_costs('Rcost', {'N': I, 'Cw': Cw}, ['R']) 208 209 return om

210 211

212 -def userfcn_reserves_int2ext(results, *args):

213 """This is the 'int2ext' stage userfcn callback that converts everything 214 back to external indexing and packages up the results. It expects to 215 find a 'reserves' field in the results struct as described for ppc 216 above, including the two additional fields 'igr' and 'rgens'. It also 217 expects the results to contain a variable 'R' and linear constraints 218 'Pg_plus_R' and 'Rreq' which are used to populate output fields in 219 results.reserves. The optional args are not currently used. 220 """ 221 ## initialize some things 222 r = results['reserves'] 223 224 ## grab some info in internal indexing order 225 igr = r['igr'] ## indices of gens available to provide reserves 226 ng = results['gen'].shape[0] ## number of on-line gens (+ disp loads) 227 228 ##----- convert stuff back to external indexing ----- 229 ## convert all reserve parameters (zones, costs, qty, rgens) 230 if 'qty' in r: 231 results = int2ext(results, ['reserves', 'qty'], ordering='gen') 232 233 results = int2ext(results, ['reserves', 'cost'], ordering='gen') 234 results = int2ext(results, ['reserves', 'zones'], ordering='gen', dim=1) 235 results = int2ext(results, ['reserves', 'rgens'], ordering='gen', dim=1) 236 results['order']['int']['reserves']['igr'] = results['reserves']['igr'] ## save internal version 237 results['reserves']['igr'] = results['order']['ext']['reserves']['igr'] ## use external version 238 r = results['reserves'] ## update 239 o = results['order'] ## update 240 241 ## grab same info in external indexing order 242 igr0 = r['igr'] ## indices of gens available to provide reserves 243 ng0 = o['ext']['gen'].shape[0] ## number of gens (+ disp loads) 244 245 ##----- results post-processing ----- 246 ## get the results (per gen reserves, multipliers) with internal gen indexing 247 ## and convert from p.u. to per MW units 248 _, Rl, Ru = results['om'].getv('R') 249 R = zeros(ng) 250 Rmin = zeros(ng) 251 Rmax = zeros(ng) 252 mu_l = zeros(ng) 253 mu_u = zeros(ng) 254 mu_Pmax = zeros(ng) 255 R[igr] = results['var']['val']['R'] * results['baseMVA'] 256 Rmin[igr] = Rl * results['baseMVA'] 257 Rmax[igr] = Ru * results['baseMVA'] 258 mu_l[igr] = results['var']['mu']['l']['R'] / results['baseMVA'] 259 mu_u[igr] = results['var']['mu']['u']['R'] / results['baseMVA'] 260 mu_Pmax[igr] = results['lin']['mu']['u']['Pg_plus_R'] / results['baseMVA'] 261 262 ## store in results in results struct 263 z = zeros(ng0) 264 results['reserves']['R'] = int2ext(results, R, z, 'gen') 265 results['reserves']['Rmin'] = int2ext(results, Rmin, z, 'gen') 266 results['reserves']['Rmax'] = int2ext(results, Rmax, z, 'gen') 267 if 'mu' not in results['reserves']: 268 results['reserves']['mu'] = {} 269 results['reserves']['mu']['l'] = int2ext(results, mu_l, z, 'gen') 270 results['reserves']['mu']['u'] = int2ext(results, mu_u, z, 'gen') 271 results['reserves']['mu']['Pmax'] = int2ext(results, mu_Pmax, z, 'gen') 272 results['reserves']['prc'] = z 273 for k in igr0: 274 iz = find(r['zones'][:, k]) 275 results['reserves']['prc'][k] = max(results['lin']['mu']['l']['Rreq'][iz]) / results['baseMVA'] 276 277 results['reserves']['totalcost'] = results['cost']['Rcost'] 278 279 ## replace ng x 1 cost, qty with ngr x 1 originals 280 if 'original' in r: 281 if 'qty' in r: 282 results['reserves']['qty'] = r['original']['qty'] 283 results['reserves']['cost'] = r['original']['cost'] 284 del results['reserves']['original'] 285 286 return results

287 288

289 -def userfcn_reserves_printpf(results, fd, ppopt, *args):

290 """This is the 'printpf' stage userfcn callback that pretty-prints the 291 results. It expects a C{results} dict, a file descriptor and a PYPOWER 292 options vector. The optional args are not currently used. 293 """ 294 ##----- print results ----- 295 r = results['reserves'] 296 nrz = r['req'].shape[0] 297 OUT_ALL = ppopt['OUT_ALL'] 298 if OUT_ALL != 0: 299 fd.write('\n================================================================================') 300 fd.write('\n| Reserves |') 301 fd.write('\n================================================================================') 302 fd.write('\n Gen Bus Status Reserves Price') 303 fd.write('\n # # (MW) ($/MW) Included in Zones ...') 304 fd.write('\n---- ----- ------ -------- -------- ------------------------') 305 for k in r['igr']: 306 iz = find(r['zones'][:, k]) 307 fd.write('\n%3d %6d %2d ' % (k, results['gen'][k, GEN_BUS], results['gen'][k, GEN_STATUS])) 308 if (results['gen'][k, GEN_STATUS] > 0) & (abs(results['reserves']['R'][k]) > 1e-6): 309 fd.write('%10.2f' % results['reserves']['R'][k]) 310 else: 311 fd.write(' - ') 312 313 fd.write('%10.2f ' % results['reserves']['prc'][k]) 314 for i in range(len(iz)): 315 if i != 0: 316 fd.write(', ') 317 fd.write('%d' % iz[i]) 318 319 fd.write('\n --------') 320 fd.write('\n Total:%10.2f Total Cost: $%.2f' % 321 (sum(results['reserves']['R'][r['igr']]), results['reserves']['totalcost'])) 322 fd.write('\n') 323 324 fd.write('\nZone Reserves Price ') 325 fd.write('\n # (MW) ($/MW) ') 326 fd.write('\n---- -------- --------') 327 for k in range(nrz): 328 iz = find(r['zones'][k, :]) ## gens in zone k 329 fd.write('\n%3d%10.2f%10.2f' % (k, sum(results['reserves']['R'][iz]), 330 results['lin']['mu']['l']['Rreq'][k] / results['baseMVA'])) 331 fd.write('\n') 332 333 fd.write('\n================================================================================') 334 fd.write('\n| Reserve Limits |') 335 fd.write('\n================================================================================') 336 fd.write('\n Gen Bus Status Rmin mu Rmin Reserves Rmax Rmax mu Pmax mu ') 337 fd.write('\n # # ($/MW) (MW) (MW) (MW) ($/MW) ($/MW) ') 338 fd.write('\n---- ----- ------ -------- -------- -------- -------- -------- --------') 339 for k in r['igr']: 340 fd.write('\n%3d %6d %2d ' % (k, results['gen'][k, GEN_BUS], results['gen'][k, GEN_STATUS])) 341 if (results['gen'][k, GEN_STATUS] > 0) & (results['reserves']['mu']['l'][k] > 1e-6): 342 fd.write('%10.2f' % results['reserves']['mu']['l'][k]) 343 else: 344 fd.write(' - ') 345 346 fd.write('%10.2f' % results['reserves']['Rmin'][k]) 347 if (results['gen'][k, GEN_STATUS] > 0) & (abs(results['reserves']['R'][k]) > 1e-6): 348 fd.write('%10.2f' % results['reserves']['R'][k]) 349 else: 350 fd.write(' - ') 351 352 fd.write('%10.2f' % results['reserves']['Rmax'][k]) 353 if (results['gen'][k, GEN_STATUS] > 0) & (results['reserves']['mu']['u'][k] > 1e-6): 354 fd.write('%10.2f' % results['reserves']['mu']['u'][k]) 355 else: 356 fd.write(' - ') 357 358 if (results['gen'][k, GEN_STATUS] > 0) & (results['reserves']['mu']['Pmax'][k] > 1e-6): 359 fd.write('%10.2f' % results['reserves']['mu']['Pmax'][k]) 360 else: 361 fd.write(' - ') 362 363 fd.write('\n --------') 364 fd.write('\n Total:%10.2f' % sum(results['reserves']['R'][r['igr']])) 365 fd.write('\n') 366 367 return results

368 369

370 -def userfcn_reserves_savecase(ppc, fd, prefix, *args):

371 """This is the 'savecase' stage userfcn callback that prints the Python 372 file code to save the 'reserves' field in the case file. It expects a 373 PYPOWER case dict (ppc), a file descriptor and variable prefix 374 (usually 'ppc'). The optional args are not currently used. 375 """ 376 r = ppc['reserves'] 377 378 fd.write('\n####----- Reserve Data -----####\n') 379 fd.write('#### reserve zones, element i, j is 1 if gen j is in zone i, 0 otherwise\n') 380 fd.write('%sreserves.zones = [\n' % prefix) 381 template = '' 382 for _ in range(r['zones'].shape[1]): 383 template = template + '\t%d' 384 template = template + ';\n' 385 fd.write(template, r.zones.T) 386 fd.write('];\n') 387 388 fd.write('\n#### reserve requirements for each zone in MW\n') 389 fd.write('%sreserves.req = [\t%g' % (prefix, r['req'][0])) 390 if len(r['req']) > 1: 391 fd.write(';\t%g' % r['req'][1:]) 392 fd.write('\t];\n') 393 394 fd.write('\n#### reserve costs in $/MW for each gen that belongs to at least 1 zone\n') 395 fd.write('#### (same order as gens, but skipping any gen that does not belong to any zone)\n') 396 fd.write('%sreserves.cost = [\t%g' % (prefix, r['cost'][0])) 397 if len(r['cost']) > 1: 398 fd.write(';\t%g' % r['cost'][1:]) 399 fd.write('\t];\n') 400 401 if 'qty' in r: 402 fd.write('\n#### OPTIONAL max reserve quantities for each gen that belongs to at least 1 zone\n') 403 fd.write('#### (same order as gens, but skipping any gen that does not belong to any zone)\n') 404 fd.write('%sreserves.qty = [\t%g' % (prefix, r['qty'][0])) 405 if len(r['qty']) > 1: 406 fd.write(';\t%g' % r['qty'][1:]) 407 fd.write('\t];\n') 408 409 ## save output fields for solved case 410 if 'R' in r: 411 fd.write('\n#### solved values\n') 412 fd.write('%sreserves.R = %s\n' % (prefix, pprint(r['R']))) 413 fd.write('%sreserves.Rmin = %s\n' % (prefix, pprint(r['Rmin']))) 414 fd.write('%sreserves.Rmax = %s\n' % (prefix, pprint(r['Rmax']))) 415 fd.write('%sreserves.mu.l = %s\n' % (prefix, pprint(r['mu']['l']))) 416 fd.write('%sreserves.mu.u = %s\n' % (prefix, pprint(r['mu']['u']))) 417 fd.write('%sreserves.prc = %s\n' % (prefix, pprint(r['prc']))) 418 fd.write('%sreserves.totalcost = %s\n' % (prefix, pprint(r['totalcost']))) 419 420 return ppc

421

Source Code for Module pypower.toggle_reserves