e2nIEE / pandapower
1
# Copyright 1996-2015 PSERC. All rights reserved.
2
# Use of this source code is governed by a BSD-style
3
# license that can be found in the LICENSE file.
4

5 1
"""Solves a DC optimal power flow.
6
"""
7

8 1
from sys import stderr
9

10 1
from copy import deepcopy
11

12 1
from numpy import \
13
    array, zeros, ones, any, diag, r_, pi, Inf, isnan, arange, c_, dot
14

15 1
from numpy import flatnonzero as find
16

17 1
from scipy.sparse import vstack, hstack, csr_matrix as sparse
18

19 1
from pandapower.pypower.idx_bus import BUS_TYPE, REF, VA, LAM_P, LAM_Q, MU_VMAX, MU_VMIN
20 1
from pandapower.pypower.idx_gen import PG, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN
21 1
from pandapower.pypower.idx_brch import PF, PT, QF, QT, RATE_A, MU_SF, MU_ST
22 1
from pandapower.pypower.idx_cost import MODEL, POLYNOMIAL, PW_LINEAR, NCOST, COST
23

24 1
from pandapower.pypower.util import sub2ind, have_fcn
25
#from pandapower.pypower.ipopt_options import ipopt_options
26
#from pandapower.pypower.cplex_options import cplex_options
27
#from pandapower.pypower.mosek_options import mosek_options
28
#from pandapower.pypower.gurobi_options import gurobi_options
29 1
from pandapower.pypower.qps_pypower import qps_pypower
30

31

32 1
def dcopf_solver(om, ppopt, out_opt=None):
33
    """Solves a DC optimal power flow.
34

35
    Inputs are an OPF model object, a PYPOWER options dict and
36
    a dict containing fields (can be empty) for each of the desired
37
    optional output fields.
38

39
    Outputs are a C{results} dict, C{success} flag and C{raw} output dict.
40

41
    C{results} is a PYPOWER case dict (ppc) with the usual baseMVA, bus
42
    branch, gen, gencost fields, along with the following additional
43
    fields:
44
        - C{order}      see 'help ext2int' for details of this field
45
        - C{x}          final value of optimization variables (internal order)
46
        - C{f}          final objective function value
47
        - C{mu}         shadow prices on ...
48
            - C{var}
49
                - C{l}  lower bounds on variables
50
                - C{u}  upper bounds on variables
51
            - C{lin}
52
                - C{l}  lower bounds on linear constraints
53
                - C{u}  upper bounds on linear constraints
54
        - C{g}          (optional) constraint values
55
        - C{dg}         (optional) constraint 1st derivatives
56
        - C{df}         (optional) obj fun 1st derivatives (not yet implemented)
57
        - C{d2f}        (optional) obj fun 2nd derivatives (not yet implemented)
58

59
    C{success} is C{True} if solver converged successfully, C{False} otherwise.
60

61
    C{raw} is a raw output dict in form returned by MINOS
62
        - C{xr}     final value of optimization variables
63
        - C{pimul}  constraint multipliers
64
        - C{info}   solver specific termination code
65
        - C{output} solver specific output information
66

67
    @see: L{opf}, L{qps_pypower}
68

69
    @author: Ray Zimmerman (PSERC Cornell)
70
    @author: Carlos E. Murillo-Sanchez (PSERC Cornell & Universidad
71
    Autonoma de Manizales)
72
    @author: Richard Lincoln
73
    """
74 1
    if out_opt is None:
75 1
        out_opt = {}
76

77
    ## options
78 1
    verbose = ppopt['VERBOSE']
79 1
    alg     = ppopt['OPF_ALG_DC']
80

81 1
    if alg == 0:
82 1
        if have_fcn('cplex'):        ## use CPLEX by default, if available
83 0
            alg = 500
84
        # MOSEK is currently not supported
85
#        elif have_fcn('mosek'):      ## if not, then MOSEK, if available
86
#            alg = 600
87 1
        elif have_fcn('gurobi'):     ## if not, then Gurobi, if available
88
            # Error in Gurobi pypower solver -> Issue with pypower 5.1.4. Gurobi won't work. Using alg 200 instead
89
            # Reason for failure: In qps_gurobi of pypower len(H) raises Error:
90
            # TypeError: sparse matrix length is ambiguous; use getnnz() or shape[0]
91
            # Todo: Fix Gurobi and activate 700 again. ATM: Fallback on 200
92
            # alg = 700
93 0
            alg = 200
94 0
            UserWarning("Gurobi not working with pypower 5.1.4")
95
        else:                        ## otherwise PIPS
96 1
            alg = 200
97

98
    ## unpack data
99 1
    ppc = om.get_ppc()
100 1
    baseMVA, bus, gen, branch, gencost = \
101
        ppc["baseMVA"], ppc["bus"], ppc["gen"], ppc["branch"], ppc["gencost"]
102 1
    cp = om.get_cost_params()
103 1
    N, H, Cw = cp["N"], cp["H"], cp["Cw"]
104 1
    fparm = array(c_[cp["dd"], cp["rh"], cp["kk"], cp["mm"]])
105 1
    Bf = om.userdata('Bf')
106 1
    Pfinj = om.userdata('Pfinj')
107 1
    vv, ll, _, _ = om.get_idx()
108

109
    ## problem dimensions
110 1
    ipol = find(gencost[:, MODEL] == POLYNOMIAL) ## polynomial costs
111 1
    ipwl = find(gencost[:, MODEL] == PW_LINEAR)  ## piece-wise linear costs
112 1
    nb = bus.shape[0]              ## number of buses
113 1
    nl = branch.shape[0]           ## number of branches
114 1
    nw = N.shape[0]                ## number of general cost vars, w
115 1
    ny = om.getN('var', 'y')       ## number of piece-wise linear costs
116 1
    nxyz = om.getN('var')          ## total number of control vars of all types
117

118
    ## linear constraints & variable bounds
119 1
    A, l, u = om.linear_constraints()
120 1
    x0, xmin, xmax = om.getv()
121

122
    ## set up objective function of the form: f = 1/2 * X'*HH*X + CC'*X
123
    ## where X = [x;y;z]. First set up as quadratic function of w,
124
    ## f = 1/2 * w'*HHw*w + CCw'*w, where w = diag(M) * (N*X - Rhat). We
125
    ## will be building on the (optionally present) user supplied parameters.
126

127
    ## piece-wise linear costs
128 1
    any_pwl = int(ny > 0)
129 1
    if any_pwl:
130
        # Sum of y vars.
131 1
        Npwl = sparse((ones(ny), (zeros(ny), arange(vv["i1"]["y"], vv["iN"]["y"]))), (1, nxyz))
132 1
        Hpwl = sparse((1, 1))
133 1
        Cpwl = array([1])
134 1
        fparm_pwl = array([[1, 0, 0, 1]])
135
    else:
136 1
        Npwl = None#zeros((0, nxyz))
137 1
        Hpwl = None#array([])
138 1
        Cpwl = array([])
139 1
        fparm_pwl = zeros((0, 4))
140

141
    ## quadratic costs
142 1
    npol = len(ipol)
143 1
    if any(len(gencost[ipol, NCOST] > 3)) and sum(gencost[find(gencost[ipol, NCOST] > 3)][:][NCOST+1:]):
144 0
        stderr.write('DC opf cannot handle polynomial costs with higher '
145
                     'than quadratic order.\n')
146 1
    iqdr = find(gencost[ipol, NCOST] == 3)
147 1
    ilin = find(gencost[ipol, NCOST] == 2)
148 1
    polycf = zeros((npol, 3))         ## quadratic coeffs for Pg
149 1
    if len(iqdr) > 0:
150 0
        polycf[iqdr, :] = gencost[ipol[iqdr], COST:COST + 3]
151 1
    if npol:
152 1
        polycf[ilin, 1:3] = gencost[ipol[ilin], COST:COST + 2]
153 1
    polycf = dot(polycf, diag([ baseMVA**2, baseMVA, 1]))     ## convert to p.u.
154 1
    if npol:
155 1
        Npol = sparse((ones(npol), (arange(npol), vv["i1"]["Pg"] + ipol)),
156
                      (npol, nxyz))  # Pg vars
157 1
        Hpol = sparse((2 * polycf[:, 0], (arange(npol), arange(npol))),
158
                      (npol, npol))
159
    else:
160 1
        Npol = None
161 1
        Hpol = None
162 1
    Cpol = polycf[:, 1]
163 1
    fparm_pol = ones((npol, 1)) * array([[1, 0, 0, 1]])
164

165
    ## combine with user costs
166 1
    NN = vstack([n for n in [Npwl, Npol, N] if n is not None and n.shape[0] > 0], "csr")
167
    # FIXME: Zero dimension sparse matrices.
168 1
    if (Hpwl is not None) and any_pwl and (npol + nw):
169 0
        Hpwl = hstack([Hpwl, sparse((any_pwl, npol + nw))])
170 1
    if Hpol is not None:
171 1
        if any_pwl and npol:
172 0
            Hpol = hstack([sparse((npol, any_pwl)), Hpol])
173 1
        if npol and nw:
174 0
            Hpol = hstack([Hpol, sparse((npol, nw))])
175 1
    if (H is not None) and nw and (any_pwl + npol):
176 0
        H = hstack([sparse((nw, any_pwl + npol)), H])
177 1
    HHw = vstack([h for h in [Hpwl, Hpol, H] if h is not None and h.shape[0] > 0], "csr")
178 1
    CCw = r_[Cpwl, Cpol, Cw]
179 1
    ffparm = r_[fparm_pwl, fparm_pol, fparm]
180

181
    ## transform quadratic coefficients for w into coefficients for X
182 1
    nnw = any_pwl + npol + nw
183 1
    M = sparse((ffparm[:, 3], (range(nnw), range(nnw))))
184 1
    MR = M * ffparm[:, 1]
185 1
    HMR = HHw * MR
186 1
    MN = M * NN
187 1
    HH = MN.T * HHw * MN
188 1
    CC = MN.T * (CCw - HMR)
189 1
    C0 = 0.5 * dot(MR, HMR) + sum(polycf[:, 2])  # Constant term of cost.
190

191
    ## set up input for QP solver
192 1
    opt = {'alg': alg, 'verbose': verbose}
193 1
    if (alg == 200) or (alg == 250):
194
        ## try to select an interior initial point
195 1
        Varefs = bus[bus[:, BUS_TYPE] == REF, VA] * (pi / 180.0)
196

197 1
        lb, ub = xmin.copy(), xmax.copy()
198 1
        lb[xmin == -Inf] = -1e10   ## replace Inf with numerical proxies
199 1
        ub[xmax ==  Inf] =  1e10
200 1
        x0 = (lb + ub) / 2;
201
        # angles set to first reference angle
202 1
        x0[vv["i1"]["Va"]:vv["iN"]["Va"]] = Varefs[0]
203 1
        if ny > 0:
204 1
            ipwl = find(gencost[:, MODEL] == PW_LINEAR)
205
            # largest y-value in CCV data
206 1
            c = gencost.flatten('F')[sub2ind(gencost.shape, ipwl,
207
                                NCOST + 2 * gencost[ipwl, NCOST])]
208 1
            x0[vv["i1"]["y"]:vv["iN"]["y"]] = max(c) + 0.1 * abs(max(c))
209

210
        ## set up options
211 1
        feastol = ppopt['PDIPM_FEASTOL']
212 1
        gradtol = ppopt['PDIPM_GRADTOL']
213 1
        comptol = ppopt['PDIPM_COMPTOL']
214 1
        costtol = ppopt['PDIPM_COSTTOL']
215 1
        max_it  = ppopt['PDIPM_MAX_IT']
216 1
        max_red = ppopt['SCPDIPM_RED_IT']
217 1
        if feastol == 0:
218 1
            feastol = ppopt['OPF_VIOLATION']    ## = OPF_VIOLATION by default
219 1
        opt["pips_opt"] = {  'feastol': feastol,
220
                             'gradtol': gradtol,
221
                             'comptol': comptol,
222
                             'costtol': costtol,
223
                             'max_it':  max_it,
224
                             'max_red': max_red,
225
                             'cost_mult': 1  }
226
#    elif alg == 400:
227
#        opt['ipopt_opt'] = ipopt_options([], ppopt)
228
#    elif alg == 500:
229
#        opt['cplex_opt'] = cplex_options([], ppopt)
230
#    elif alg == 600:
231
#        opt['mosek_opt'] = mosek_options([], ppopt)
232
#    elif alg == 700:
233
#        ppopt['GRB_OPT'] = 0
234
#        ppopt['GRB_METHOD'] = "automatic"
235
#        ppopt['GRB_TIMELIMIT'] = Inf
236
#        ppopt['GRB_THREADS'] = 0
237
#        opt['GRB_OPT'] = gurobi_options(None, ppopt)
238
#    else:
239
#        raise ValueError("Unrecognised solver [%d]." % alg)
240

241
    ##-----  run opf  -----
242 1
    x, f, info, output, lmbda = \
243
            qps_pypower(HH, CC, A, l, u, xmin, xmax, x0, opt)
244 1
    success = (info == 1)
245

246
    ##-----  calculate return values  -----
247 1
    if not any(isnan(x)):
248
        ## update solution data
249 1
        Va = x[vv["i1"]["Va"]:vv["iN"]["Va"]]
250 1
        Pg = x[vv["i1"]["Pg"]:vv["iN"]["Pg"]]
251 1
        f = f + C0
252

253
        ## update voltages & generator outputs
254 1
        bus[:, VA] = Va * 180 / pi
255 1
        gen[:, PG] = Pg * baseMVA
256

257
        ## compute branch flows
258 1
        branch[:, [QF, QT]] = zeros((nl, 2))
259 1
        branch[:, PF] = (Bf * Va + Pfinj) * baseMVA
260 1
        branch[:, PT] = -branch[:, PF]
261

262
    ## package up results
263 1
    mu_l = lmbda["mu_l"]
264 1
    mu_u = lmbda["mu_u"]
265 1
    muLB = lmbda["lower"]
266 1
    muUB = lmbda["upper"]
267

268
    ## update Lagrange multipliers
269 1
    il = find((branch[:, RATE_A] != 0) & (branch[:, RATE_A] < 1e10))
270 1
    bus[:, [LAM_P, LAM_Q, MU_VMIN, MU_VMAX]] = zeros((nb, 4))
271 1
    gen[:, [MU_PMIN, MU_PMAX, MU_QMIN, MU_QMAX]] = zeros((gen.shape[0], 4))
272 1
    branch[:, [MU_SF, MU_ST]] = zeros((nl, 2))
273 1
    bus[:, LAM_P]       = (mu_u[ll["i1"]["Pmis"]:ll["iN"]["Pmis"]] -
274
                           mu_l[ll["i1"]["Pmis"]:ll["iN"]["Pmis"]]) / baseMVA
275 1
    branch[il, MU_SF]   = mu_u[ll["i1"]["Pf"]:ll["iN"]["Pf"]] / baseMVA
276 1
    branch[il, MU_ST]   = mu_u[ll["i1"]["Pt"]:ll["iN"]["Pt"]] / baseMVA
277 1
    gen[:, MU_PMIN]     = muLB[vv["i1"]["Pg"]:vv["iN"]["Pg"]] / baseMVA
278 1
    gen[:, MU_PMAX]     = muUB[vv["i1"]["Pg"]:vv["iN"]["Pg"]] / baseMVA
279

280 1
    pimul = r_[
281
      mu_l - mu_u,
282
     -ones((ny)), ## dummy entry corresponding to linear cost row in A
283
      muLB - muUB
284
    ]
285

286 1
    mu = { 'var': {'l': muLB, 'u': muUB},
287
           'lin': {'l': mu_l, 'u': mu_u} }
288

289 1
    results = deepcopy(ppc)
290 1
    results["bus"], results["branch"], results["gen"], \
291
        results["om"], results["x"], results["mu"], results["f"] = \
292
            bus, branch, gen, om, x, mu, f
293

294 1
    raw = {'xr': x, 'pimul': pimul, 'info': info, 'output': output}
295

296 1
    return results, success, raw

Read our documentation on viewing source code .

Loading