pypower.dcopf

44 """Solves a DC optimal power flow. 45 46 Inputs are an OPF model object, a PYPOWER options dict and 47 a dict containing fields (can be empty) for each of the desired 48 optional output fields. 49 50 Outputs are a C{results} dict, C{success} flag and C{raw} output dict. 51 52 C{results} is a PYPOWER case dict (ppc) with the usual baseMVA, bus 53 branch, gen, gencost fields, along with the following additional 54 fields: 55 - C{order} see 'help ext2int' for details of this field 56 - C{x} final value of optimization variables (internal order) 57 - C{f} final objective function value 58 - C{mu} shadow prices on ... 59 - C{var} 60 - C{l} lower bounds on variables 61 - C{u} upper bounds on variables 62 - C{lin} 63 - C{l} lower bounds on linear constraints 64 - C{u} upper bounds on linear constraints 65 - C{g} (optional) constraint values 66 - C{dg} (optional) constraint 1st derivatives 67 - C{df} (optional) obj fun 1st derivatives (not yet implemented) 68 - C{d2f} (optional) obj fun 2nd derivatives (not yet implemented) 69 70 C{success} is C{True} if solver converged successfully, C{False} otherwise. 71 72 C{raw} is a raw output dict in form returned by MINOS 73 - C{xr} final value of optimization variables 74 - C{pimul} constraint multipliers 75 - C{info} solver specific termination code 76 - C{output} solver specific output information 77 78 @see: L{opf}, L{qps_pypower} 79 80 @author: Ray Zimmerman (PSERC Cornell) 81 @author: Carlos E. Murillo-Sanchez (PSERC Cornell & Universidad 82 Autonoma de Manizales) 83 @author: Richard Lincoln 84 """ 85 if out_opt is None: 86 out_opt = {} 87 88 ## options 89 verbose = ppopt['VERBOSE'] 90 alg = ppopt['OPF_ALG_DC'] 91 92 if alg == 0: 93 alg = 200 # Use PIPS 94 95 ## unpack data 96 ppc = om.get_ppc() 97 baseMVA, bus, gen, branch, gencost = \ 98 ppc["baseMVA"], ppc["bus"], ppc["gen"], ppc["branch"], ppc["gencost"] 99 cp = om.get_cost_params() 100 N, H, Cw = cp["N"], cp["H"], cp["Cw"] 101 fparm = array(c_[cp["dd"], cp["rh"], cp["kk"], cp["mm"]]) 102 Bf = om.userdata('Bf') 103 Pfinj = om.userdata('Pfinj') 104 vv, ll, _, _ = om.get_idx() 105 106 ## problem dimensions 107 ipol = find(gencost[:, MODEL] == POLYNOMIAL) ## polynomial costs 108 ipwl = find(gencost[:, MODEL] == PW_LINEAR) ## piece-wise linear costs 109 nb = bus.shape[0] ## number of buses 110 nl = branch.shape[0] ## number of branches 111 nw = N.shape[0] ## number of general cost vars, w 112 ny = om.getN('var', 'y') ## number of piece-wise linear costs 113 nxyz = om.getN('var') ## total number of control vars of all types 114 115 ## linear constraints & variable bounds 116 A, l, u = om.linear_constraints() 117 x0, xmin, xmax = om.getv() 118 119 ## set up objective function of the form: f = 1/2 * X'*HH*X + CC'*X 120 ## where X = [x;y;z]. First set up as quadratic function of w, 121 ## f = 1/2 * w'*HHw*w + CCw'*w, where w = diag(M) * (N*X - Rhat). We 122 ## will be building on the (optionally present) user supplied parameters. 123 124 ## piece-wise linear costs 125 any_pwl = int(ny > 0) 126 if any_pwl: 127 # Sum of y vars. 128 Npwl = sparse((ones(ny), (zeros(ny), arange(vv["i1"]["y"], vv["iN"]["y"]))), (1, nxyz)) 129 Hpwl = sparse((1, 1)) 130 Cpwl = array([1]) 131 fparm_pwl = array([[1, 0, 0, 1]]) 132 else: 133 Npwl = None#zeros((0, nxyz)) 134 Hpwl = None#array([]) 135 Cpwl = array([]) 136 fparm_pwl = zeros((0, 4)) 137 138 ## quadratic costs 139 npol = len(ipol) 140 if any(find(gencost[ipol, NCOST] > 3)): 141 stderr.write('DC opf cannot handle polynomial costs with higher ' 142 'than quadratic order.\n') 143 iqdr = find(gencost[ipol, NCOST] == 3) 144 ilin = find(gencost[ipol, NCOST] == 2) 145 polycf = zeros((npol, 3)) ## quadratic coeffs for Pg 146 if len(iqdr) > 0: 147 polycf[iqdr, :] = gencost[ipol[iqdr], COST:COST + 3] 148 if npol: 149 polycf[ilin, 1:3] = gencost[ipol[ilin], COST:COST + 2] 150 polycf = dot(polycf, diag([ baseMVA**2, baseMVA, 1])) ## convert to p.u. 151 if npol: 152 Npol = sparse((ones(npol), (arange(npol), vv["i1"]["Pg"] + ipol)), 153 (npol, nxyz)) # Pg vars 154 Hpol = sparse((2 * polycf[:, 0], (arange(npol), arange(npol))), 155 (npol, npol)) 156 else: 157 Npol = None 158 Hpol = None 159 Cpol = polycf[:, 1] 160 fparm_pol = ones((npol, 1)) * array([[1, 0, 0, 1]]) 161 162 ## combine with user costs 163 NN = vstack([n for n in [Npwl, Npol, N] if n is not None and n.shape[0] > 0], "csr") 164 # FIXME: Zero dimension sparse matrices. 165 if (Hpwl is not None) and any_pwl and (npol + nw): 166 Hpwl = hstack([Hpwl, sparse((any_pwl, npol + nw))]) 167 if Hpol is not None: 168 if any_pwl and npol: 169 Hpol = hstack([sparse((npol, any_pwl)), Hpol]) 170 if npol and nw: 171 Hpol = hstack([Hpol, sparse((npol, nw))]) 172 if (H is not None) and nw and (any_pwl + npol): 173 H = hstack([sparse((nw, any_pwl + npol)), H]) 174 HHw = vstack([h for h in [Hpwl, Hpol, H] if h is not None and h.shape[0] > 0], "csr") 175 CCw = r_[Cpwl, Cpol, Cw] 176 ffparm = r_[fparm_pwl, fparm_pol, fparm] 177 178 ## transform quadratic coefficients for w into coefficients for X 179 nnw = any_pwl + npol + nw 180 M = sparse((ffparm[:, 3], (range(nnw), range(nnw)))) 181 MR = M * ffparm[:, 1] 182 HMR = HHw * MR 183 MN = M * NN 184 HH = MN.T * HHw * MN 185 CC = MN.T * (CCw - HMR) 186 C0 = 0.5 * dot(MR, HMR) + sum(polycf[:, 2]) # Constant term of cost. 187 188 ## set up input for QP solver 189 opt = {'alg': alg, 'verbose': verbose} 190 if (alg == 200) or (alg == 250): 191 ## try to select an interior initial point 192 Varefs = bus[bus[:, BUS_TYPE] == REF, VA] * (pi / 180.0) 193 194 lb, ub = xmin.copy(), xmax.copy() 195 lb[xmin == -Inf] = -1e10 ## replace Inf with numerical proxies 196 ub[xmax == Inf] = 1e10 197 x0 = (lb + ub) / 2; 198 # angles set to first reference angle 199 x0[vv["i1"]["Va"]:vv["iN"]["Va"]] = Varefs[0] 200 if ny > 0: 201 ipwl = find(gencost[:, MODEL] == PW_LINEAR) 202 # largest y-value in CCV data 203 c = gencost.flatten('F')[sub2ind(gencost.shape, ipwl, 204 NCOST + 2 * gencost[ipwl, NCOST])] 205 x0[vv["i1"]["y"]:vv["iN"]["y"]] = max(c) + 0.1 * abs(max(c)) 206 207 ## set up options 208 feastol = ppopt['PDIPM_FEASTOL'] 209 gradtol = ppopt['PDIPM_GRADTOL'] 210 comptol = ppopt['PDIPM_COMPTOL'] 211 costtol = ppopt['PDIPM_COSTTOL'] 212 max_it = ppopt['PDIPM_MAX_IT'] 213 max_red = ppopt['SCPDIPM_RED_IT'] 214 if feastol == 0: 215 feastol = ppopt['OPF_VIOLATION'] ## = OPF_VIOLATION by default 216 opt["pips_opt"] = { 'feastol': feastol, 217 'gradtol': gradtol, 218 'comptol': comptol, 219 'costtol': costtol, 220 'max_it': max_it, 221 'max_red': max_red, 222 'cost_mult': 1 } 223 elif alg == 400: 224 opt['ipopt_opt'] = ipopt_options([], ppopt) 225 elif alg == 500: 226 opt['cplex_opt'] = cplex_options([], ppopt) 227 elif alg == 600: 228 opt['mosek_opt'] = mosek_options([], ppopt) 229 else: 230 raise ValueError, "Unrecognised solver [%d]." % alg 231 232 ##----- run opf ----- 233 x, f, info, output, lmbda = \ 234 qps_pypower(HH, CC, A, l, u, xmin, xmax, x0, opt) 235 success = (info == 1) 236 237 ##----- calculate return values ----- 238 if not any(isnan(x)): 239 ## update solution data 240 Va = x[vv["i1"]["Va"]:vv["iN"]["Va"]] 241 Pg = x[vv["i1"]["Pg"]:vv["iN"]["Pg"]] 242 f = f + C0 243 244 ## update voltages & generator outputs 245 bus[:, VA] = Va * 180 / pi 246 gen[:, PG] = Pg * baseMVA 247 248 ## compute branch flows 249 branch[:, [QF, QT]] = zeros((nl, 2)) 250 branch[:, PF] = (Bf * Va + Pfinj) * baseMVA 251 branch[:, PT] = -branch[:, PF] 252 253 ## package up results 254 mu_l = lmbda["mu_l"] 255 mu_u = lmbda["mu_u"] 256 muLB = lmbda["lower"] 257 muUB = lmbda["upper"] 258 259 ## update Lagrange multipliers 260 il = find((branch[:, RATE_A] != 0) & (branch[:, RATE_A] < 1e10)) 261 bus[:, [LAM_P, LAM_Q, MU_VMIN, MU_VMAX]] = zeros((nb, 4)) 262 gen[:, [MU_PMIN, MU_PMAX, MU_QMIN, MU_QMAX]] = zeros((gen.shape[0], 4)) 263 branch[:, [MU_SF, MU_ST]] = zeros((nl, 2)) 264 bus[:, LAM_P] = (mu_u[ll["i1"]["Pmis"]:ll["iN"]["Pmis"]] - 265 mu_l[ll["i1"]["Pmis"]:ll["iN"]["Pmis"]]) / baseMVA 266 branch[il, MU_SF] = mu_u[ll["i1"]["Pf"]:ll["iN"]["Pf"]] / baseMVA 267 branch[il, MU_ST] = mu_u[ll["i1"]["Pt"]:ll["iN"]["Pt"]] / baseMVA 268 gen[:, MU_PMIN] = muLB[vv["i1"]["Pg"]:vv["iN"]["Pg"]] / baseMVA 269 gen[:, MU_PMAX] = muUB[vv["i1"]["Pg"]:vv["iN"]["Pg"]] / baseMVA 270 271 pimul = r_[ 272 mu_l - mu_u, 273 -ones((ny > 0)), ## dummy entry corresponding to linear cost row in A 274 muLB - muUB 275 ] 276 277 mu = { 'var': {'l': muLB, 'u': muUB}, 278 'lin': {'l': mu_l, 'u': mu_u} } 279 280 results = deepcopy(ppc) 281 results["bus"], results["branch"], results["gen"], \ 282 results["om"], results["x"], results["mu"], results["f"] = \ 283 bus, branch, gen, om, x, mu, f 284 285 raw = {'xr': x, 'pimul': pimul, 'info': info, 'output': output} 286 287 return results, success, raw

Source Code for Module pypower.dcopf_solver