pypower.pipsopf

36 """Solves AC optimal power flow using PIPS. 37 38 Inputs are an OPF model object, a PYPOWER options vector and 39 a dict containing keys (can be empty) for each of the desired 40 optional output fields. 41 42 outputs are a C{results} dict, C{success} flag and C{raw} output dict. 43 44 C{results} is a PYPOWER case dict (ppc) with the usual baseMVA, bus 45 branch, gen, gencost fields, along with the following additional 46 fields: 47 - C{order} see 'help ext2int' for details of this field 48 - C{x} final value of optimization variables (internal order) 49 - C{f} final objective function value 50 - C{mu} shadow prices on ... 51 - C{var} 52 - C{l} lower bounds on variables 53 - C{u} upper bounds on variables 54 - C{nln} 55 - C{l} lower bounds on nonlinear constraints 56 - C{u} upper bounds on nonlinear constraints 57 - C{lin} 58 - C{l} lower bounds on linear constraints 59 - C{u} upper bounds on linear constraints 60 61 C{success} is C{True} if solver converged successfully, C{False} otherwise 62 63 C{raw} is a raw output dict in form returned by MINOS 64 - xr final value of optimization variables 65 - pimul constraint multipliers 66 - info solver specific termination code 67 - output solver specific output information 68 69 @see: L{opf}, L{pips} 70 71 @author: Ray Zimmerman (PSERC Cornell) 72 @author: Carlos E. Murillo-Sanchez (PSERC Cornell & Universidad 73 Autonoma de Manizales) 74 @author: Richard Lincoln 75 """ 76 ##----- initialization ----- 77 ## optional output 78 if out_opt is None: 79 out_opt = {} 80 81 ## options 82 verbose = ppopt['VERBOSE'] 83 feastol = ppopt['PDIPM_FEASTOL'] 84 gradtol = ppopt['PDIPM_GRADTOL'] 85 comptol = ppopt['PDIPM_COMPTOL'] 86 costtol = ppopt['PDIPM_COSTTOL'] 87 max_it = ppopt['PDIPM_MAX_IT'] 88 max_red = ppopt['SCPDIPM_RED_IT'] 89 step_control = (ppopt['OPF_ALG'] == 565) ## OPF_ALG == 565, PIPS-sc 90 if feastol == 0: 91 feastol = ppopt['OPF_VIOLATION'] 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 = bus.shape[0] ## number of buses 110 nl = branch.shape[0] ## number of branches 111 ny = om.getN('var', 'y') ## number of piece-wise linear costs 112 113 ## linear constraints 114 A, l, u = om.linear_constraints() 115 116 ## bounds on optimization vars 117 _, xmin, xmax = om.getv() 118 119 ## build admittance matrices 120 Ybus, Yf, Yt = makeYbus(baseMVA, bus, branch) 121 122 ## try to select an interior initial point 123 ll, uu = xmin.copy(), xmax.copy() 124 ll[xmin == -Inf] = -1e10 ## replace Inf with numerical proxies 125 uu[xmax == Inf] = 1e10 126 x0 = (ll + uu) / 2 127 Varefs = bus[bus[:, BUS_TYPE] == REF, VA] * (pi / 180) 128 ## angles set to first reference angle 129 x0[vv["i1"]["Va"]:vv["iN"]["Va"]] = Varefs[0] 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 c = gencost.flatten('F')[sub2ind(gencost.shape, ipwl, NCOST+2*gencost[ipwl, NCOST])] ## largest y-value in CCV data 135 x0[vv["i1"]["y"]:vv["iN"]["y"]] = max(c) + 0.1 * abs(max(c)) 136 # x0[vv["i1"]["y"]:vv["iN"]["y"]] = c + 0.1 * abs(c) 137 138 ## find branches with flow limits 139 il = find((branch[:, RATE_A] != 0) & (branch[:, RATE_A] < 1e10)) 140 nl2 = len(il) ## number of constrained lines 141 142 ##----- run opf ----- 143 f_fcn = lambda x, return_hessian=False: opf_costfcn(x, om, return_hessian) 144 gh_fcn = lambda x: opf_consfcn(x, om, Ybus, Yf[il, :], Yt[il,:], ppopt, il) 145 hess_fcn = lambda x, lmbda, cost_mult: opf_hessfcn(x, lmbda, om, Ybus, Yf[il, :], Yt[il, :], ppopt, il, cost_mult) 146 147 solution = pips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt) 148 x, f, info, lmbda, output = solution["x"], solution["f"], \ 149 solution["eflag"], solution["lmbda"], solution["output"] 150 151 success = (info > 0) 152 153 ## update solution data 154 Va = x[vv["i1"]["Va"]:vv["iN"]["Va"]] 155 Vm = x[vv["i1"]["Vm"]:vv["iN"]["Vm"]] 156 Pg = x[vv["i1"]["Pg"]:vv["iN"]["Pg"]] 157 Qg = x[vv["i1"]["Qg"]:vv["iN"]["Qg"]] 158 159 V = Vm * exp(1j * Va) 160 161 ##----- calculate return values ----- 162 ## update voltages & generator outputs 163 bus[:, VA] = Va * 180 / pi 164 bus[:, VM] = Vm 165 gen[:, PG] = Pg * baseMVA 166 gen[:, QG] = Qg * baseMVA 167 gen[:, VG] = Vm[ gen[:, GEN_BUS].astype(int) ] 168 169 ## compute branch flows 170 Sf = V[ branch[:, F_BUS].astype(int) ] * conj(Yf * V) ## cplx pwr at "from" bus, p["u"]. 171 St = V[ branch[:, T_BUS].astype(int) ] * conj(Yt * V) ## cplx pwr at "to" bus, p["u"]. 172 branch[:, PF] = Sf.real * baseMVA 173 branch[:, QF] = Sf.imag * baseMVA 174 branch[:, PT] = St.real * baseMVA 175 branch[:, QT] = St.imag * baseMVA 176 177 ## line constraint is actually on square of limit 178 ## so we must fix multipliers 179 muSf = zeros(nl) 180 muSt = zeros(nl) 181 if len(il) > 0: 182 muSf[il] = \ 183 2 * lmbda["ineqnonlin"][:nl2] * branch[il, RATE_A] / baseMVA 184 muSt[il] = \ 185 2 * lmbda["ineqnonlin"][nl2:nl2+nl2] * branch[il, RATE_A] / baseMVA 186 187 ## update Lagrange multipliers 188 bus[:, MU_VMAX] = lmbda["upper"][vv["i1"]["Vm"]:vv["iN"]["Vm"]] 189 bus[:, MU_VMIN] = lmbda["lower"][vv["i1"]["Vm"]:vv["iN"]["Vm"]] 190 gen[:, MU_PMAX] = lmbda["upper"][vv["i1"]["Pg"]:vv["iN"]["Pg"]] / baseMVA 191 gen[:, MU_PMIN] = lmbda["lower"][vv["i1"]["Pg"]:vv["iN"]["Pg"]] / baseMVA 192 gen[:, MU_QMAX] = lmbda["upper"][vv["i1"]["Qg"]:vv["iN"]["Qg"]] / baseMVA 193 gen[:, MU_QMIN] = lmbda["lower"][vv["i1"]["Qg"]:vv["iN"]["Qg"]] / baseMVA 194 195 bus[:, LAM_P] = \ 196 lmbda["eqnonlin"][nn["i1"]["Pmis"]:nn["iN"]["Pmis"]] / baseMVA 197 bus[:, LAM_Q] = \ 198 lmbda["eqnonlin"][nn["i1"]["Qmis"]:nn["iN"]["Qmis"]] / baseMVA 199 branch[:, MU_SF] = muSf / baseMVA 200 branch[:, MU_ST] = muSt / baseMVA 201 202 ## package up results 203 nlnN = om.getN('nln') 204 205 ## extract multipliers for nonlinear constraints 206 kl = find(lmbda["eqnonlin"] < 0) 207 ku = find(lmbda["eqnonlin"] > 0) 208 nl_mu_l = zeros(nlnN) 209 nl_mu_u = r_[zeros(2*nb), muSf, muSt] 210 nl_mu_l[kl] = -lmbda["eqnonlin"][kl] 211 nl_mu_u[ku] = lmbda["eqnonlin"][ku] 212 213 mu = { 214 'var': {'l': lmbda["lower"], 'u': lmbda["upper"]}, 215 'nln': {'l': nl_mu_l, 'u': nl_mu_u}, 216 'lin': {'l': lmbda["mu_l"], 'u': lmbda["mu_u"]} } 217 218 results = ppc 219 results["bus"], results["branch"], results["gen"], \ 220 results["om"], results["x"], results["mu"], results["f"] = \ 221 bus, branch, gen, om, x, mu, f 222 223 pimul = r_[ 224 results["mu"]["nln"]["l"] - results["mu"]["nln"]["u"], 225 results["mu"]["lin"]["l"] - results["mu"]["lin"]["u"], 226 -ones(ny > 0), 227 results["mu"]["var"]["l"] - results["mu"]["var"]["u"], 228 ] 229 raw = {'xr': x, 'pimul': pimul, 'info': info, 'output': output} 230 231 return results, success, raw

Source Code for Module pypower.pipsopf_solver