e2nIEE / pandapower
1
# Copyright (c) 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
"""Evaluates nonlinear constraints and their Jacobian for OPF.
6
"""
7

8 1
from numpy import zeros, ones, conj, exp, r_, Inf, arange
9

10 1
from scipy.sparse import lil_matrix, vstack, hstack, csr_matrix as sparse
11

12 1
from pandapower.pypower.idx_gen import GEN_BUS, PG, QG
13 1
from pandapower.pypower.idx_brch import F_BUS, T_BUS, RATE_A
14

15 1
from pandapower.pypower.makeSbus import makeSbus
16 1
from pandapower.pypower.dSbus_dV import dSbus_dV
17 1
from pandapower.pypower.dIbr_dV import dIbr_dV
18 1
from pandapower.pypower.dSbr_dV import dSbr_dV
19 1
from pandapower.pypower.dAbr_dV import dAbr_dV
20

21

22 1
def opf_consfcn(x, om, Ybus, Yf, Yt, ppopt, il=None, *args):
23
    """Evaluates nonlinear constraints and their Jacobian for OPF.
24

25
    Constraint evaluation function for AC optimal power flow, suitable
26
    for use with L{pips}. Computes constraint vectors and their gradients.
27

28
    @param x: optimization vector
29
    @param om: OPF model object
30
    @param Ybus: bus admittance matrix
31
    @param Yf: admittance matrix for "from" end of constrained branches
32
    @param Yt: admittance matrix for "to" end of constrained branches
33
    @param ppopt: PYPOWER options vector
34
    @param il: (optional) vector of branch indices corresponding to
35
    branches with flow limits (all others are assumed to be
36
    unconstrained). The default is C{range(nl)} (all branches).
37
    C{Yf} and C{Yt} contain only the rows corresponding to C{il}.
38

39
    @return: C{h} - vector of inequality constraint values (flow limits)
40
    limit^2 - flow^2, where the flow can be apparent power real power or
41
    current, depending on value of C{OPF_FLOW_LIM} in C{ppopt} (only for
42
    constrained lines). C{g} - vector of equality constraint values (power
43
    balances). C{dh} - (optional) inequality constraint gradients, column
44
    j is gradient of h(j). C{dg} - (optional) equality constraint gradients.
45

46
    @see: L{opf_costfcn}, L{opf_hessfcn}
47

48
    @author: Carlos E. Murillo-Sanchez (PSERC Cornell & Universidad
49
    Autonoma de Manizales)
50
    @author: Ray Zimmerman (PSERC Cornell)
51
    """
52
    ##----- initialize -----
53

54
    ## unpack data
55 1
    ppc = om.get_ppc()
56 1
    baseMVA, bus, gen, branch = \
57
        ppc["baseMVA"], ppc["bus"], ppc["gen"], ppc["branch"]
58 1
    vv, _, _, _ = om.get_idx()
59

60
    ## problem dimensions
61 1
    nb = bus.shape[0]          ## number of buses
62 1
    nl = branch.shape[0]       ## number of branches
63 1
    ng = gen.shape[0]          ## number of dispatchable injections
64 1
    nxyz = len(x)              ## total number of control vars of all types
65

66
    ## set default constrained lines
67 1
    if il is None:
68 0
        il = arange(nl)         ## all lines have limits by default
69 1
    nl2 = len(il)              ## number of constrained lines
70

71
    ## grab Pg & Qg
72 1
    Pg = x[vv["i1"]["Pg"]:vv["iN"]["Pg"]]  ## active generation in p.u.
73 1
    Qg = x[vv["i1"]["Qg"]:vv["iN"]["Qg"]]  ## reactive generation in p.u.
74

75
    ## put Pg & Qg back in gen
76 1
    gen[:, PG] = Pg * baseMVA  ## active generation in MW
77 1
    gen[:, QG] = Qg * baseMVA  ## reactive generation in MVAr
78

79
    ## rebuild Sbus
80 1
    Sbus = makeSbus(baseMVA, bus, gen) ## net injected power in p.u.
81

82
    ## ----- evaluate constraints -----
83
    ## reconstruct V
84 1
    Va = x[vv["i1"]["Va"]:vv["iN"]["Va"]]
85 1
    Vm = x[vv["i1"]["Vm"]:vv["iN"]["Vm"]]
86 1
    V = Vm * exp(1j * Va)
87

88
    ## evaluate power flow equations
89 1
    mis = V * conj(Ybus * V) - Sbus
90

91
    ##----- evaluate constraint function values -----
92
    ## first, the equality constraints (power flow)
93 1
    g = r_[ mis.real,            ## active power mismatch for all buses
94
            mis.imag ]           ## reactive power mismatch for all buses
95

96
    ## then, the inequality constraints (branch flow limits)
97 1
    if nl2 > 0:
98 1
        flow_max = (branch[il, RATE_A] / baseMVA)**2
99 1
        flow_max[flow_max == 0] = Inf
100 1
        if ppopt['OPF_FLOW_LIM'] == 2:       ## current magnitude limit, |I|
101 1
            If = Yf * V
102 1
            It = Yt * V
103 1
            h = r_[ If * conj(If) - flow_max,     ## branch I limits (from bus)
104
                    It * conj(It) - flow_max ].real    ## branch I limits (to bus)
105
        else:
106
            ## compute branch power flows
107
            ## complex power injected at "from" bus (p.u.)
108 0
            Sf = V[ branch[il, F_BUS].astype(int) ] * conj(Yf * V)
109
            ## complex power injected at "to" bus (p.u.)
110 0
            St = V[ branch[il, T_BUS].astype(int) ] * conj(Yt * V)
111 0
            if ppopt['OPF_FLOW_LIM'] == 1:   ## active power limit, P (Pan Wei)
112 0
                h = r_[ Sf.real**2 - flow_max,   ## branch P limits (from bus)
113
                        St.real**2 - flow_max ]  ## branch P limits (to bus)
114
            else:                ## apparent power limit, |S|
115 0
                h = r_[ Sf * conj(Sf) - flow_max, ## branch S limits (from bus)
116
                        St * conj(St) - flow_max ].real  ## branch S limits (to bus)
117
    else:
118 1
        h = zeros((0,1))
119

120
    ##----- evaluate partials of constraints -----
121
    ## index ranges
122 1
    iVa = arange(vv["i1"]["Va"], vv["iN"]["Va"])
123 1
    iVm = arange(vv["i1"]["Vm"], vv["iN"]["Vm"])
124 1
    iPg = arange(vv["i1"]["Pg"], vv["iN"]["Pg"])
125 1
    iQg = arange(vv["i1"]["Qg"], vv["iN"]["Qg"])
126 1
    iVaVmPgQg = r_[iVa, iVm, iPg, iQg].T
127

128
    ## compute partials of injected bus powers
129 1
    dSbus_dVm, dSbus_dVa = dSbus_dV(Ybus, V)           ## w.r.t. V
130
    ## Pbus w.r.t. Pg, Qbus w.r.t. Qg
131 1
    neg_Cg = sparse((-ones(ng), (gen[:, GEN_BUS], range(ng))), (nb, ng))
132

133
    ## construct Jacobian of equality constraints (power flow) and transpose it
134 1
    dg = lil_matrix((2 * nb, nxyz))
135 1
    blank = sparse((nb, ng))
136 1
    dg[:, iVaVmPgQg] = vstack([
137
            ## P mismatch w.r.t Va, Vm, Pg, Qg
138
            hstack([dSbus_dVa.real, dSbus_dVm.real, neg_Cg, blank]),
139
            ## Q mismatch w.r.t Va, Vm, Pg, Qg
140
            hstack([dSbus_dVa.imag, dSbus_dVm.imag, blank, neg_Cg])
141
        ], "csr")
142 1
    dg = dg.T
143

144 1
    if nl2 > 0:
145
        ## compute partials of Flows w.r.t. V
146 1
        if ppopt['OPF_FLOW_LIM'] == 2:     ## current
147 1
            dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft = \
148
                    dIbr_dV(branch[il, :], Yf, Yt, V)
149
        else:                  ## power
150 0
            dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft = \
151
                    dSbr_dV(branch[il, :], Yf, Yt, V)
152 1
        if ppopt['OPF_FLOW_LIM'] == 1:     ## real part of flow (active power)
153 0
            dFf_dVa = dFf_dVa.real
154 0
            dFf_dVm = dFf_dVm.real
155 0
            dFt_dVa = dFt_dVa.real
156 0
            dFt_dVm = dFt_dVm.real
157 0
            Ff = Ff.real
158 0
            Ft = Ft.real
159

160
        ## squared magnitude of flow (of complex power or current, or real power)
161 1
        df_dVa, df_dVm, dt_dVa, dt_dVm = \
162
                dAbr_dV(dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft)
163

164
        ## construct Jacobian of inequality constraints (branch limits)
165
        ## and transpose it.
166 1
        dh = lil_matrix((2 * nl2, nxyz))
167 1
        dh[:, r_[iVa, iVm].T] = vstack([
168
                hstack([df_dVa, df_dVm]),    ## "from" flow limit
169
                hstack([dt_dVa, dt_dVm])     ## "to" flow limit
170
            ], "csr")
171 1
        dh = dh.T
172
    else:
173 1
        dh = None
174

175 1
    return h, g, dh, dg

Read our documentation on viewing source code .

Loading