e2nIEE / pandapower
1
# -*- coding: utf-8 -*-
2

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

7
# Copyright (c) 2016-2021 by University of Kassel and Fraunhofer Institute for Energy Economics
8
# and Energy System Technology (IEE), Kassel. All rights reserved.
9

10

11 1
"""Evaluates Hessian of Lagrangian for AC OPF.
12
"""
13

14 1
from numpy import array, zeros, ones, exp, arange, r_, flatnonzero as find
15 1
from scipy.sparse import vstack, hstack, issparse, csr_matrix as sparse
16

17 1
from pandapower.pypower.d2AIbr_dV2 import d2AIbr_dV2
18 1
from pandapower.pypower.d2ASbr_dV2 import d2ASbr_dV2
19 1
from pandapower.pypower.d2Sbus_dV2 import d2Sbus_dV2
20 1
from pandapower.pypower.dIbr_dV import dIbr_dV
21 1
from pandapower.pypower.dSbr_dV import dSbr_dV
22 1
from pandapower.pypower.idx_brch import F_BUS, T_BUS
23 1
from pandapower.pypower.idx_cost import MODEL, POLYNOMIAL
24 1
from pandapower.pypower.idx_gen import PG, QG
25 1
from pandapower.pypower.opf_consfcn import opf_consfcn
26 1
from pandapower.pypower.opf_costfcn import opf_costfcn
27 1
from pandapower.pypower.polycost import polycost
28

29

30 1
def opf_hessfcn(x, lmbda, om, Ybus, Yf, Yt, ppopt, il=None, cost_mult=1.0):
31
    """Evaluates Hessian of Lagrangian for AC OPF.
32

33
    Hessian evaluation function for AC optimal power flow, suitable
34
    for use with L{pips}.
35

36
    Examples::
37
        Lxx = opf_hessfcn(x, lmbda, om, Ybus, Yf, Yt, ppopt)
38
        Lxx = opf_hessfcn(x, lmbda, om, Ybus, Yf, Yt, ppopt, il)
39
        Lxx = opf_hessfcn(x, lmbda, om, Ybus, Yf, Yt, ppopt, il, cost_mult)
40

41
    @param x: optimization vector
42
    @param lmbda: C{eqnonlin} - Lagrange multipliers on power balance
43
    equations. C{ineqnonlin} - Kuhn-Tucker multipliers on constrained
44
    branch flows.
45
    @param om: OPF model object
46
    @param Ybus: bus admittance matrix
47
    @param Yf: admittance matrix for "from" end of constrained branches
48
    @param Yt: admittance matrix for "to" end of constrained branches
49
    @param ppopt: PYPOWER options vector
50
    @param il: (optional) vector of branch indices corresponding to
51
    branches with flow limits (all others are assumed to be unconstrained).
52
    The default is C{range(nl)} (all branches). C{Yf} and C{Yt} contain
53
    only the rows corresponding to C{il}.
54
    @param cost_mult: (optional) Scale factor to be applied to the cost
55
    (default = 1).
56

57
    @return: Hessian of the Lagrangian.
58

59
    @see: L{opf_costfcn}, L{opf_consfcn}
60

61
    @author: Ray Zimmerman (PSERC Cornell)
62
    @author: Carlos E. Murillo-Sanchez (PSERC Cornell & Universidad
63
    Autonoma de Manizales)
64
    @author: Richard Lincoln
65

66
    Modified by University of Kassel (Friederike Meier): Bugfix in line 173
67
    """
68
    ##----- initialize -----
69
    ## unpack data
70 1
    ppc = om.get_ppc()
71 1
    baseMVA, bus, gen, branch, gencost = \
72
        ppc["baseMVA"], ppc["bus"], ppc["gen"], ppc["branch"], ppc["gencost"]
73 1
    cp = om.get_cost_params()
74 1
    N, Cw, H, dd, rh, kk, mm = \
75
        cp["N"], cp["Cw"], cp["H"], cp["dd"], cp["rh"], cp["kk"], cp["mm"]
76 1
    vv, _, _, _ = om.get_idx()
77

78
    ## unpack needed parameters
79 1
    nb = bus.shape[0]          ## number of buses
80 1
    nl = branch.shape[0]       ## number of branches
81 1
    ng = gen.shape[0]          ## number of dispatchable injections
82 1
    nxyz = len(x)              ## total number of control vars of all types
83

84
    ## set default constrained lines
85 1
    if il is None:
86 0
        il = arange(nl)            ## all lines have limits by default
87 1
    nl2 = len(il)           ## number of constrained lines
88

89
    ## grab Pg & Qg
90 1
    Pg = x[vv["i1"]["Pg"]:vv["iN"]["Pg"]]  ## active generation in p.u.
91 1
    Qg = x[vv["i1"]["Qg"]:vv["iN"]["Qg"]]  ## reactive generation in p.u.
92

93
    ## put Pg & Qg back in gen
94 1
    gen[:, PG] = Pg * baseMVA  ## active generation in MW
95 1
    gen[:, QG] = Qg * baseMVA  ## reactive generation in MVAr
96

97
    ## reconstruct V
98 1
    Va = x[vv["i1"]["Va"]:vv["iN"]["Va"]]
99 1
    Vm = x[vv["i1"]["Vm"]:vv["iN"]["Vm"]]
100 1
    V = Vm * exp(1j * Va)
101 1
    nxtra = nxyz - 2 * nb
102 1
    pcost = gencost[arange(ng), :]
103 1
    if gencost.shape[0] > ng:
104 1
        qcost = gencost[arange(ng, 2 * ng), :]
105
    else:
106 1
        qcost = array([])
107

108
    ## ----- evaluate d2f -----
109 1
    d2f_dPg2 = zeros(ng)#sparse((ng, 1))               ## w.r.t. p.u. Pg
110 1
    d2f_dQg2 = zeros(ng)#sparse((ng, 1))               ## w.r.t. p.u. Qg
111 1
    ipolp = find(pcost[:, MODEL] == POLYNOMIAL)
112 1
    if len(ipolp):
113 1
        d2f_dPg2[ipolp] = \
114
            baseMVA**2 * polycost(pcost[ipolp, :], Pg[ipolp] * baseMVA, 2)
115 1
    if qcost.any():          ## Qg is not free
116 1
        ipolq = find(qcost[:, MODEL] == POLYNOMIAL)
117 1
        d2f_dQg2[ipolq] = \
118
                baseMVA**2 * polycost(qcost[ipolq, :], Qg[ipolq] * baseMVA, 2)
119 1
    i = r_[arange(vv["i1"]["Pg"], vv["iN"]["Pg"]),
120
           arange(vv["i1"]["Qg"], vv["iN"]["Qg"])]
121
#    d2f = sparse((vstack([d2f_dPg2, d2f_dQg2]).toarray().flatten(),
122
#                  (i, i)), shape=(nxyz, nxyz))
123 1
    d2f = sparse((r_[d2f_dPg2, d2f_dQg2], (i, i)), (nxyz, nxyz))
124

125
    ## generalized cost
126
    if issparse(N) and N.nnz > 0: # pragma: no cover
127
        nw = N.shape[0]
128
        r = N * x - rh                    ## Nx - rhat
129
        iLT = find(r < -kk)               ## below dead zone
130
        iEQ = find((r == 0) & (kk == 0))  ## dead zone doesn't exist
131
        iGT = find(r > kk)                ## above dead zone
132
        iND = r_[iLT, iEQ, iGT]           ## rows that are Not in the Dead region
133
        iL = find(dd == 1)                ## rows using linear function
134
        iQ = find(dd == 2)                ## rows using quadratic function
135
        LL = sparse((ones(len(iL)), (iL, iL)), (nw, nw))
136
        QQ = sparse((ones(len(iQ)), (iQ, iQ)), (nw, nw))
137
        kbar = sparse((r_[ones(len(iLT)), zeros(len(iEQ)), -ones(len(iGT))],
138
                       (iND, iND)), (nw, nw)) * kk
139
        rr = r + kbar                  ## apply non-dead zone shift
140
        M = sparse((mm[iND], (iND, iND)), (nw, nw))  ## dead zone or scale
141
        diagrr = sparse((rr, (arange(nw), arange(nw))), (nw, nw))
142

143
        ## linear rows multiplied by rr(i), quadratic rows by rr(i)^2
144
        w = M * (LL + QQ * diagrr) * rr
145
        HwC = H * w + Cw
146
        AA = N.T * M * (LL + 2 * QQ * diagrr)
147

148
        d2f = d2f + AA * H * AA.T + 2 * N.T * M * QQ * \
149
                sparse((HwC, (arange(nw), arange(nw))), (nw, nw)) * N
150 1
    d2f = d2f * cost_mult
151

152
    ##----- evaluate Hessian of power balance constraints -----
153 1
    nlam = len(lmbda["eqnonlin"]) // 2
154 1
    lamP = lmbda["eqnonlin"][:nlam]
155 1
    lamQ = lmbda["eqnonlin"][nlam:nlam + nlam]
156 1
    Gpaa, Gpav, Gpva, Gpvv = d2Sbus_dV2(Ybus, V, lamP)
157 1
    Gqaa, Gqav, Gqva, Gqvv = d2Sbus_dV2(Ybus, V, lamQ)
158

159 1
    d2G = vstack([
160
            hstack([
161
                vstack([hstack([Gpaa, Gpav]),
162
                        hstack([Gpva, Gpvv])]).real +
163
                vstack([hstack([Gqaa, Gqav]),
164
                        hstack([Gqva, Gqvv])]).imag,
165
                sparse((2 * nb, nxtra))]),
166
            hstack([
167
                sparse((nxtra, 2 * nb)),
168
                sparse((nxtra, nxtra))
169
            ])
170
        ], "csr")
171

172
    ##----- evaluate Hessian of flow constraints -----
173 1
    nmu = len(lmbda["ineqnonlin"]) // 2
174 1
    muF = lmbda["ineqnonlin"][:nmu]
175 1
    muT = lmbda["ineqnonlin"][nmu:nmu + nmu]
176 1
    if ppopt['OPF_FLOW_LIM'] == 2:       ## current
177 1
        if Yf.size:
178 1
            dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It = dIbr_dV(branch, Yf, Yt, V)
179 1
            Hfaa, Hfav, Hfva, Hfvv = d2AIbr_dV2(dIf_dVa, dIf_dVm, If, Yf, V, muF)
180 1
            Htaa, Htav, Htva, Htvv = d2AIbr_dV2(dIt_dVa, dIt_dVm, It, Yt, V, muT)
181
        else:
182 1
            Hfaa= Hfav= Hfva= Hfvv= Htaa= Htav= Htva= Htvv = sparse(zeros((nb,nb)))
183
    else: # pragma: no cover
184
        f = branch[il, F_BUS].astype(int)    ## list of "from" buses
185
        t = branch[il, T_BUS].astype(int)    ## list of "to" buses
186
        ## connection matrix for line & from buses
187
        Cf = sparse((ones(nl2), (arange(nl2), f)), (nl2, nb))
188
        ## connection matrix for line & to buses
189
        Ct = sparse((ones(nl2), (arange(nl2), t)), (nl2, nb))
190
        dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St = \
191
                dSbr_dV(branch[il,:], Yf, Yt, V)
192
        if ppopt['OPF_FLOW_LIM'] == 1:     ## real power
193
            Hfaa, Hfav, Hfva, Hfvv = d2ASbr_dV2(dSf_dVa.real, dSf_dVm.real,
194
                                                Sf.real, Cf, Yf, V, muF)
195
            Htaa, Htav, Htva, Htvv = d2ASbr_dV2(dSt_dVa.real, dSt_dVm.real,
196
                                                St.real, Ct, Yt, V, muT)
197
        else:                  ## apparent power
198
            Hfaa, Hfav, Hfva, Hfvv = \
199
                    d2ASbr_dV2(dSf_dVa, dSf_dVm, Sf, Cf, Yf, V, muF)
200
            Htaa, Htav, Htva, Htvv = \
201
                    d2ASbr_dV2(dSt_dVa, dSt_dVm, St, Ct, Yt, V, muT)
202

203 1
    d2H = vstack([
204
            hstack([
205
                vstack([hstack([Hfaa, Hfav]),
206
				hstack([Hfva, Hfvv])]) +
207
				vstack([hstack([Htaa, Htav]),
208
                        hstack([Htva, Htvv])]),
209
                sparse((2 * nb, nxtra))
210
            ]),
211
            hstack([
212
                sparse((nxtra, 2 * nb)),
213
                sparse((nxtra, nxtra))
214
            ])
215
        ], "csr")
216

217
    ##-----  do numerical check using (central) finite differences  -----
218
    if 0:
219
        nx = len(x)
220
        step = 1e-5
221
        num_d2f = sparse((nx, nx))
222
        num_d2G = sparse((nx, nx))
223
        num_d2H = sparse((nx, nx))
224
        for i in range(nx):
225
            xp = x
226
            xm = x
227
            xp[i] = x[i] + step / 2
228
            xm[i] = x[i] - step / 2
229
            # evaluate cost & gradients
230
            _, dfp = opf_costfcn(xp, om)
231
            _, dfm = opf_costfcn(xm, om)
232
            # evaluate constraints & gradients
233
            _, _, dHp, dGp = opf_consfcn(xp, om, Ybus, Yf, Yt, ppopt, il)
234
            _, _, dHm, dGm = opf_consfcn(xm, om, Ybus, Yf, Yt, ppopt, il)
235
            num_d2f[:, i] = cost_mult * (dfp - dfm) / step
236
            num_d2G[:, i] = (dGp - dGm) * lmbda["eqnonlin"]   / step
237
            num_d2H[:, i] = (dHp - dHm) * lmbda["ineqnonlin"] / step
238
        d2f_err = max(max(abs(d2f - num_d2f)))
239
        d2G_err = max(max(abs(d2G - num_d2G)))
240
        d2H_err = max(max(abs(d2H - num_d2H)))
241
        if d2f_err > 1e-6:
242
            print('Max difference in d2f: %g' % d2f_err)
243
        if d2G_err > 1e-5:
244
            print('Max difference in d2G: %g' % d2G_err)
245
        if d2H_err > 1e-6:
246
            print('Max difference in d2H: %g' % d2H_err)
247

248 1
    return d2f + d2G + d2H

Read our documentation on viewing source code .

Loading