contrib/cvxopf.py

   1 #------------------------------------------------------------------------------
   2 # Copyright (C) 2010 Richard Lincoln
   3 #
   4 # This program is free software: you can redistribute it and/or modify
   5 # it under the terms of the GNU General Public License as published by
   6 # the Free Software Foundation, either version 3 of the License, or
   7 # (at your option) any later version.
   8 #
   9 # This program is distributed in the hope that it will be useful,
  10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12 # GNU General Public License for more details.
  13 #
  14 # You should have received a copy of the GNU General Public License
  15 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
  16 #------------------------------------------------------------------------------
  17
  18 """ Defines an OPF solver using CVXOPT [1].
  19
  20     Note: This module is licensed under GNU GPL version 3 due to the
  21     CVXOPT import.
  22
  23     [1] Ray Zimmerman, "mipsopf_solver.m", MATPOWER, PSERC Cornell,
  24         version 4.0b1, http://www.pserc.cornell.edu/matpower/, Dec 2009
  25 """
  26
  27 #------------------------------------------------------------------------------
  28 #  Imports:
  29 #------------------------------------------------------------------------------
  30
  31 import logging
  32
  33 from numpy import polyval, polyder, finfo, pi
  34
  35 from cvxopt import matrix, spmatrix, mul, sparse, exp, solvers, div, spdiag
  36
  37 from pylon.solver import _Solver, DCOPFSolver, SFLOW, PFLOW, IFLOW
  38
  39 #------------------------------------------------------------------------------
  40 #  Logging:
  41 #------------------------------------------------------------------------------
  42
  43 logger = logging.getLogger(__name__)
  44
  45 #------------------------------------------------------------------------------
  46 #  Constants:
  47 #------------------------------------------------------------------------------
  48
  49 EPS = finfo(float).eps
  50
  51 #------------------------------------------------------------------------------
  52 #  "DCCVXOPTSolver" class:
  53 #------------------------------------------------------------------------------
  54
  55 class DCCVXOPTSolver(DCOPFSolver):
  56     """ Solves DC optimal power flow using CVXOPT.
  57     """
  58
  59     #--------------------------------------------------------------------------
  60     #  "object" interface:
  61     #--------------------------------------------------------------------------
  62
  63     def __init__(self, om, opt=None, solver=None):
  64         """ Initialises the new DCOPF instance.
  65         """
  66         super(DCOPFSolver, self).__init__(om, opt)
  67
  68         # Choice of solver (May be None or "mosek" (or "glpk" for linear
  69         # formulation)). Specify None to use the Python solver from CVXOPT.
  70         self.solver = solver
  71
  72         if opt.has_key("verbose"):
  73             solvers.options["show_progress"] = opt["verbose"]
  74         if opt.has_key("max_it"):
  75             solvers.options["maxiters"] = opt["max_it"]
  76         if opt.has_key("feastol"):
  77             solvers.options["feastol"] = opt["feastol"]
  78         if opt.has_key("gradtol"):
  79             raise NotImplementedError
  80         if opt.has_key("comptol"):
  81             raise NotImplementedError
  82         if opt.has_key("costtol"):
  83             raise NotImplementedError
  84         if opt.has_key("max_red"):
  85             raise NotImplementedError
  86         if opt.has_key("step_control"):
  87             raise NotImplementedError
  88         if opt.has_key("cost_mult"):
  89             raise NotImplementedError
  90
  91
  92     def _run_opf(self, P, q, AA, ll, uu, xmin, xmax, x0, opt):
  93         """ Solves the either quadratic or linear program.
  94         """
  95         AAcvx = tocvx(AA)
  96         nx = x0.shape[0] # number of variables
  97         # add var limits to linear constraints
  98         eyex = spmatrix(1.0, range(nx), range(nx))
  99         A = sparse([eyex, AAcvx])
 100         l = matrix([-xmin, ll])
 101         u = matrix([ xmax, uu])
 102
 103         Ae, be, Ai, bi = split_linear_constraints(A, l, u)
 104
 105
 106         if len(P.V) > 0:
 107             cvx_sol = solvers.qp(P, q, Ai, bi, Ae, be, self.solver, {"x": x0})
 108         else:
 109             cvx_sol = solvers.lp(q, Ai, bi, Ae, be, self.solver, {"x": x0})
 110
 111         return cvx_sol
 112
 113
 114     def _update_case(self, bs, ln, gn, base_mva, Bf, Pfinj, Va, Pg, lmbda):
 115         """ Calculates the result attribute values.
 116         """
 117         for i, bus in enumerate(bs):
 118             bus.v_angle = Va[i] * 180.0 / pi
 119
 120 #------------------------------------------------------------------------------
 121 #  "CVXOPTSolver" class:
 122 #------------------------------------------------------------------------------
 123
 124 class CVXOPTSolver(_Solver):
 125     """ Solves AC optimal power flow using convex optimization.
 126     """
 127
 128     def __init__(self, om, opt=None, flow_lim=SFLOW):
 129         """ Initialises a new CVXOPTSolver instance.
 130         """
 131         super(CVXOPTSolver, self).__init__(om)
 132
 133         # Quantity to limit for branch flow constraints ("S", "P" or "I").
 134         self.flow_lim = flow_lim
 135
 136         self.opt = {} if opt is None else opt
 137
 138         if self.opt.has_key("verbose"):
 139             solvers.options["show_progress"] = self.opt["verbose"]
 140         if self.opt.has_key("max_it"):
 141             solvers.options["maxiters"] = self.opt["max_it"]
 142         if self.opt.has_key("feastol"):
 143             solvers.options["feastol"] = self.opt["feastol"]
 144         if self.opt.has_key("gradtol"):
 145             raise NotImplementedError
 146         if self.opt.has_key("comptol"):
 147             raise NotImplementedError
 148         if self.opt.has_key("costtol"):
 149             raise NotImplementedError
 150         if self.opt.has_key("max_red"):
 151             raise NotImplementedError
 152         if self.opt.has_key("step_control"):
 153             raise NotImplementedError
 154         if self.opt.has_key("cost_mult"):
 155             raise NotImplementedError
 156
 157
 158     def solve(self):
 159         case = self.om.case
 160         base_mva = case.base_mva
 161         self.opt["cost_mult"] = 1e-4
 162         # Unpack the OPF model.
 163         bs, ln, gn, _ = self._unpack_model(self.om)
 164         # Compute problem dimensions.
 165         ng = len(gn)
 166         ipol, _, nb, nl, _, ny, nxyz = self._dimension_data(bs, ln, gn)
 167         # The number of non-linear equality constraints.
 168 #        neq = 2 * nb
 169         # The number of control variables.
 170 #        nc = 2 * nb + 2 * ng
 171         # Indexes of constrained lines.
 172         il = matrix([i for i,l in enumerate(ln) if 0.0 < l.rate_a < 1e10])
 173 #        nl2 = len(il)
 174
 175         neqnln = 2 * nb # no. of non-linear equality constraints
 176         niqnln = 2 * len(il) # no. of lines with constraints
 177
 178         # Linear constraints (l <= A*x <= u).
 179         A, l, u = self.om.linear_constraints()
 180         if A is None:
 181             AA = spmatrix([], [], [], (0, self.om.var_N))
 182         else:
 183             AA = tocvx(A)
 184         Ae, be, Ai, bi = split_linear_constraints(AA, matrix(l), matrix(u))
 185
 186         _, xmin, xmax = self._var_bounds()
 187
 188         # Select an interior initial point for interior point solver.
 189         x0 = matrix(self._initial_interior_point(bs, gn, xmin, xmax, ny))
 190
 191         # Build admittance matrices.
 192         YYbus, YYf, YYt = case.Y
 193         Ybus = tocvx(YYbus)
 194         Yf = tocvx(YYf)
 195         Yt = tocvx(YYt)
 196
 197         # Optimisation variables.
 198         Va = self.om.get_var("Va")
 199         Vm = self.om.get_var("Vm")
 200         Pg = self.om.get_var("Pg")
 201         Qg = self.om.get_var("Qg")
 202
 203         def F(x=None, z=None):
 204             """ Evaluates the objective and nonlinear constraint functions.
 205             """
 206             if x is None:
 207                 # Include power mismatch constraint twice to force equality.
 208                 nln_N = self.om.nln_N + neqnln
 209                 return nln_N, x0
 210
 211             # Evaluate objective function -------------------------------------
 212
 213             Pgen = x[Pg.i1:Pg.iN + 1] # Active generation in p.u.
 214             Qgen = x[Qg.i1:Qg.iN + 1] # Reactive generation in p.u.
 215
 216             xx = matrix([Pgen, Qgen]) * base_mva
 217
 218             # Evaluate the objective function value.
 219             if len(ipol) > 0:
 220                 # FIXME: Implement reactive power costs.
 221                 f0 = sum([g.total_cost(xx[i]) for i, g in enumerate(gn)])
 222             else:
 223                 f0 = 0
 224
 225             # Piecewise linear cost of P and Q.
 226             if ny:
 227                 y = self.om.get_var("y")
 228                 ccost = spmatrix(matrix(1.0, (ny, 1)),
 229                                  range(y.i1, y.iN + 1),
 230                                  matrix(0.0, (ny, 1)), size=(nxyz, 1)).T
 231                 f0 = f0 + ccost * x
 232             else:
 233                 ccost = matrix(0.0, (1, nxyz))
 234             # TODO: Generalised cost term.
 235
 236             # Evaluate cost gradient ------------------------------------------
 237
 238             iPg = range(Pg.i1, Pg.iN + 1)
 239             iQg = range(Qg.i1, Qg.iN + 1)
 240
 241             # Polynomial cost of P and Q.
 242             df_dPgQg = matrix(0.0, (2 * ng, 1))        # w.r.t p.u. Pg and Qg
 243             for i in ipol:
 244                 df_dPgQg[i] = \
 245                     base_mva * polyval(polyder(list(gn[i].p_cost)), xx[i])
 246
 247             df0 = matrix(0.0, (nxyz, 1))
 248             df0[iPg] = df_dPgQg[:ng]
 249             df0[iQg] = df_dPgQg[ng:ng + ng]
 250
 251             # Piecewise linear cost of P and Q.
 252             df0 = df0 + ccost.T
 253             # TODO: Generalised cost term.
 254
 255             # Evaluate cost Hessian -------------------------------------------
 256
 257 #            d2f = None
 258
 259             # Evaluate nonlinear equality constraints -------------------------
 260
 261             for i, g in enumerate(gn):
 262                 g.p = Pgen[i] * base_mva # active generation in MW
 263                 g.q = Qgen[i] * base_mva # reactive generation in MVAr
 264
 265             # Rebuild the net complex bus power injection vector in p.u.
 266             Sbus = matrix(case.getSbus(bs))
 267
 268             Vang = x[Va.i1:Va.iN + 1]
 269             Vmag = x[Vm.i1:Vm.iN + 1]
 270             V = mul(Vmag, exp(1j * Vang))
 271
 272             # Evaluate the power flow equations.
 273             mis = mul(V, conj(Ybus * V)) - Sbus
 274
 275             # Equality constraints (power flow).
 276             g = matrix([mis.real(),  # active power mismatch for all buses
 277                         mis.imag()]) # reactive power mismatch for all buses
 278
 279             # Evaluate nonlinear inequality constraints -----------------------
 280
 281             # Inequality constraints (branch flow limits).
 282             # (line constraint is actually on square of limit)
 283             flow_max = matrix([(l.rate_a / base_mva)**2 for l in ln])
 284             # FIXME: There must be a more elegant way to do this.
 285             for i, rate in enumerate(flow_max):
 286                 if rate == 0.0:
 287                     flow_max[i] = 1e5
 288
 289             if self.flow_lim == IFLOW:
 290                 If = Yf * V
 291                 It = Yt * V
 292                 # Branch current limits.
 293                 h = matrix([(If * conj(If)) - flow_max,
 294                             (If * conj(It)) - flow_max])
 295             else:
 296                 i_fbus = [e.from_bus._i for e in ln]
 297                 i_tbus = [e.to_bus._i for e in ln]
 298                 # Complex power injected at "from" bus (p.u.).
 299                 Sf = mul(V[i_fbus], conj(Yf * V))
 300                 # Complex power injected at "to" bus (p.u.).
 301                 St = mul(V[i_tbus], conj(Yt * V))
 302                 if self.flow_lim == PFLOW: # active power limit, P (Pan Wei)
 303                     # Branch real power limits.
 304                     h = matrix([Sf.real()**2 - flow_max,
 305                                 St.real()**2 - flow_max])
 306                 elif self.flow_lim == SFLOW: # apparent power limit, |S|
 307                     # Branch apparent power limits.
 308                     h = matrix([mul(Sf, conj(Sf)) - flow_max,
 309                                 mul(St, conj(St)) - flow_max]).real()
 310                 else:
 311                     raise ValueError
 312
 313             # Evaluate partial derivatives of constraints ---------------------
 314
 315             iVa = matrix(range(Va.i1, Va.iN + 1))
 316             iVm = matrix(range(Vm.i1, Vm.iN + 1))
 317             iPg = matrix(range(Pg.i1, Pg.iN + 1))
 318             iQg = matrix(range(Qg.i1, Qg.iN + 1))
 319             iVaVmPgQg = matrix([iVa, iVm, iPg, iQg]).T
 320
 321             # Compute partials of injected bus powers.
 322             dSbus_dVm, dSbus_dVa = dSbus_dV(Ybus, V)
 323
 324             i_gbus = [gen.bus._i for gen in gn]
 325             neg_Cg = spmatrix(matrix(-1.0, (ng, 1)),
 326                                      i_gbus,
 327                                      range(ng), (nb, ng))
 328
 329             # Transposed Jacobian of the power balance equality constraints.
 330             dg = spmatrix([], [], [], (nxyz, 2 * nb))
 331
 332             blank = spmatrix([], [], [], (nb, ng))
 333
 334             dg[iVaVmPgQg, :] = sparse([
 335                 [dSbus_dVa.real(), dSbus_dVa.imag()],
 336                 [dSbus_dVm.real(), dSbus_dVm.imag()],
 337                 [neg_Cg, blank],
 338                 [blank, neg_Cg]
 339             ]).T
 340
 341             # Compute partials of flows w.r.t V.
 342             if self.flow_lim == IFLOW:
 343                 dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft = \
 344                     dIbr_dV(Yf, Yt, V)
 345             else:
 346                 dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft = \
 347                     dSbr_dV(Yf, Yt, V, bs, ln)
 348             if self.flow_lim == PFLOW:
 349                 dFf_dVa = dFf_dVa.real
 350                 dFf_dVm = dFf_dVm.real
 351                 dFt_dVa = dFt_dVa.real
 352                 dFt_dVm = dFt_dVm.real
 353                 Ff = Ff.real
 354                 Ft = Ft.real
 355
 356             # Squared magnitude of flow (complex power, current or real power).
 357             df_dVa, df_dVm, dt_dVa, dt_dVm = \
 358                 dAbr_dV(dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft)
 359
 360             # Construct Jacobian of inequality constraints (branch limits) and
 361             # transpose it.
 362             dh = spmatrix([], [], [], (nxyz, 2 * nl))
 363             dh[matrix([iVa, iVm]).T, :] = sparse([[df_dVa, dt_dVa],
 364                                                   [df_dVm, dt_dVm]]).T
 365
 366             f = matrix([f0, g, h, -g])
 367             df = matrix([[df0], [dg], [dh], [-dg]]).T
 368
 369             if z is None:
 370                 return f, df
 371
 372             # Evaluate cost Hessian -------------------------------------------
 373
 374             nxtra = nxyz - (2 * nb)
 375
 376             # Evaluate d2f ----------------------------------------------------
 377
 378             d2f_dPg2 = matrix(0.0, (ng, 1)) # w.r.t p.u. Pg
 379             d2f_dQg2 = matrix(0.0, (ng, 1)) # w.r.t p.u. Qg
 380
 381             for i in ipol:
 382                 d2f_dPg2[i, 0] = polyval(polyder(list(gn[i].p_cost), 2),
 383                                          Pg.v0[i] * base_mva) * base_mva**2
 384             # TODO: Reactive power costs.
 385 #            for i in ipol:
 386 #                d2f_dQg2[i] = polyval(polyder(list(gn[i].q_cost), 2),
 387 #                                      Qg.v0[i] * base_mva) * base_mva**2
 388
 389             i = matrix([range(Pg.i1, Pg.iN + 1), range(Qg.i1, Qg.iN + 1)])
 390
 391             d2f = spmatrix(matrix([d2f_dPg2, d2f_dQg2]), i, i, (nxyz, nxyz))
 392
 393             # TODO: Generalised cost model.
 394
 395             d2f = d2f * self.opt["cost_mult"]
 396
 397             # Evaluate Hessian of power balance constraints -------------------
 398
 399             eqnonlin = z[1:neqnln + 1]
 400             nlam = len(eqnonlin) / 2
 401             lamP = eqnonlin[:nlam]
 402             lamQ = eqnonlin[nlam:nlam + nlam]
 403
 404             Gpaa, Gpav, Gpva, Gpvv = d2Sbus_dV2(Ybus, V, lamP)
 405             Gqaa, Gqav, Gqva, Gqvv = d2Sbus_dV2(Ybus, V, lamQ)
 406
 407             d2G_1 = sparse([[sparse([[Gpaa, Gpva], [Gpav, Gpvv]]).real() +
 408                            sparse([[Gqaa, Gqva], [Gqav, Gqvv]]).imag()],
 409                            [spmatrix([], [], [], (2 * nb, nxtra))]])
 410             d2G_2 = spmatrix([], [], [], (nxtra, 2 * nb + nxtra))
 411             d2G = sparse([d2G_1, d2G_2])
 412
 413             #  Evaluate Hessian of flow constraints ---------------------------
 414
 415             ineqnonlin = z[1 + neqnln:1 + neqnln + niqnln]
 416
 417             nmu = len(ineqnonlin) / 2
 418             muF = ineqnonlin[:nmu]
 419             muT = ineqnonlin[nmu:nmu + nmu]
 420
 421             if self.flow_lim == IFLOW:
 422                 dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It = \
 423                     dIbr_dV(Yf, Yt, V)
 424                 Hfaa, Hfav, Hfva, Hfvv = \
 425                     d2AIbr_dV2(dIf_dVa, dIf_dVm, If, Yf, V, muF)
 426                 Htaa, Htav, Htva, Htvv = \
 427                     d2AIbr_dV2(dIt_dVa, dIt_dVm, It, Yt, V, muT)
 428             else:
 429                 fr = [e.from_bus._i for e in ln]
 430                 to = [e.to_bus._i for e in ln]
 431                 # Line-bus connection matrices.
 432                 Cf = spmatrix(1.0, range(nl), fr, (nl, nb))
 433                 Ct = spmatrix(1.0, range(nl), to, (nl, nb))
 434                 dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St = \
 435                     dSbr_dV(Yf, Yt, V, bs, ln)
 436                 if self.flow_lim == PFLOW:
 437                     Hfaa, Hfav, Hfva, Hfvv = \
 438                         d2ASbr_dV2(dSf_dVa.real(), dSf_dVm.real(),
 439                                         Sf.real(), Cf, Yf, V, muF)
 440                     Htaa, Htav, Htva, Htvv = \
 441                         d2ASbr_dV2(dSt_dVa.real(), dSt_dVm.real(),
 442                                         St.real(), Ct, Yt, V, muT)
 443                 elif self.flow_lim == SFLOW:
 444                     Hfaa, Hfav, Hfva, Hfvv = \
 445                         d2ASbr_dV2(dSf_dVa, dSf_dVm, Sf, Cf, Yf, V, muF)
 446                     Htaa, Htav, Htva, Htvv = \
 447                         d2ASbr_dV2(dSt_dVa, dSt_dVm, St, Ct, Yt, V, muT)
 448                 else:
 449                     raise ValueError
 450
 451             d2H_1 = sparse([[sparse([[Hfaa, Hfva], [Hfav, Hfvv]]) +
 452                              sparse([[Htaa, Htva], [Htav, Htvv]])],
 453                             [spmatrix([], [], [], (2 * nb, nxtra))]])
 454             d2H_2 = spmatrix([], [], [], (nxtra, 2 * nb + nxtra))
 455             d2H = sparse([d2H_1, d2H_2])
 456
 457             Lxx = d2f + d2G + d2H
 458
 459             return f, df, Lxx
 460
 461         # cp(F, G=None, h=None, dims=None, A=None, b=None, kktsolver=None)
 462         #
 463         #     minimize    f0(x)
 464         #     subject to  fk(x) <= 0, k = 1, ..., mnl
 465         #                 G*x   <= h
 466         #                 A*x   =  b.
 467         solution = solvers.cp(F, G=Ai, h=bi, dims=None, A=Ae, b=be)
 468
 469         return solution
 470
 471
 472 def split_linear_constraints(A, l, u):
 473     """ Returns the linear equality and inequality constraints.
 474     """
 475     ieq = []
 476     igt = []
 477     ilt = []
 478     ibx = []
 479     for i in range(len(l)):
 480         if abs(u[i] - l[i]) <= EPS:
 481             ieq.append(i)
 482         elif (u[i] > 1e10) and (l[i] > -1e10):
 483             igt.append(i)
 484         elif (l[i] <= -1e10) and (u[i] < 1e10):
 485             ilt.append(i)
 486         elif (abs(u[i] - l[i]) > EPS) and (u[i] < 1e10) and (l[i] > -1e10):
 487             ibx.append(i)
 488         else:
 489             raise ValueError
 490
 491     Ae = A[ieq, :]
 492     Ai = sparse([A[ilt, :], -A[igt, :], A[ibx, :], -A[ibx, :]])
 493     be = u[ieq, :]
 494     bi = matrix([u[ilt], -l[igt], u[ibx], -l[ibx]])
 495
 496     return Ae, be, Ai, bi
 497
 498 #--------------------------------------------------------------------------
 499 #  Partial derivative of power injection w.r.t. voltage:
 500 #--------------------------------------------------------------------------
 501
 502 def dSbus_dV(Y, V):
 503     """ Computes the partial derivative of power injection w.r.t. voltage.
 504
 505         References:
 506             Ray Zimmerman, "dSbus_dV.m", MATPOWER, version 3.2,
 507             PSERC (Cornell), http://www.pserc.cornell.edu/matpower/
 508     """
 509     I = Y * V
 510
 511     diagV = spdiag(V)
 512     diagIbus = spdiag(I)
 513     diagVnorm = spdiag(div(V, abs(V))) # Element-wise division.
 514
 515     dS_dVm = diagV * conj(Y * diagVnorm) + conj(diagIbus) * diagVnorm
 516     dS_dVa = 1j * diagV * conj(diagIbus - Y * diagV)
 517
 518     return dS_dVm, dS_dVa
 519
 520 #--------------------------------------------------------------------------
 521 #  Partial derivatives of branch currents w.r.t. voltage.
 522 #--------------------------------------------------------------------------
 523
 524 def dIbr_dV(Yf, Yt, V):
 525     """ Computes partial derivatives of branch currents w.r.t. voltage.
 526
 527         Ray Zimmerman, "dIbr_dV.m", MATPOWER, version 4.0b1,
 528         PSERC (Cornell), http://www.pserc.cornell.edu/matpower/
 529     """
 530 #        nb = len(V)
 531
 532     Vnorm = div(V, abs(V))
 533     diagV = spdiag(V)
 534     diagVnorm = spdiag(Vnorm)
 535     dIf_dVa = Yf * 1j * diagV
 536     dIf_dVm = Yf * diagVnorm
 537     dIt_dVa = Yt * 1j * diagV
 538     dIt_dVm = Yt * diagVnorm
 539
 540     # Compute currents.
 541     If = Yf * V
 542     It = Yt * V
 543
 544     return dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It
 545
 546 #--------------------------------------------------------------------------
 547 #  Partial derivative of branch power flow w.r.t voltage:
 548 #--------------------------------------------------------------------------
 549
 550 def dSbr_dV(Yf, Yt, V, buses, branches):
 551     """ Computes the branch power flow vector and the partial derivative of
 552         branch power flow w.r.t voltage.
 553     """
 554     nl = len(branches)
 555     nb = len(V)
 556
 557     f = matrix([l.from_bus._i for l in branches])
 558     t = matrix([l.to_bus._i for l in branches])
 559
 560     # Compute currents.
 561     If = Yf * V
 562     It = Yt * V
 563
 564     Vnorm = div(V, abs(V))
 565
 566     diagVf = spdiag(V[f])
 567     diagIf = spdiag(If)
 568     diagVt = spdiag(V[t])
 569     diagIt = spdiag(It)
 570     diagV = spdiag(V)
 571     diagVnorm = spdiag(Vnorm)
 572
 573     ibr = range(nl)
 574     size = (nl, nb)
 575     # Partial derivative of S w.r.t voltage phase angle.
 576     dSf_dVa = 1j * (conj(diagIf) *
 577         spmatrix(V[f], ibr, f, size) - diagVf * conj(Yf * diagV))
 578
 579     dSt_dVa = 1j * (conj(diagIt) *
 580         spmatrix(V[t], ibr, t, size) - diagVt * conj(Yt * diagV))
 581
 582     # Partial derivative of S w.r.t. voltage amplitude.
 583     dSf_dVm = diagVf * conj(Yf * diagVnorm) + conj(diagIf) * \
 584         spmatrix(Vnorm[f], ibr, f, size)
 585
 586     dSt_dVm = diagVt * conj(Yt * diagVnorm) + conj(diagIt) * \
 587         spmatrix(Vnorm[t], ibr, t, size)
 588
 589     # Compute power flow vectors.
 590     Sf = mul(V[f], conj(If))
 591     St = mul(V[t], conj(It))
 592
 593     return dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St
 594
 595 #--------------------------------------------------------------------------
 596 #  Partial derivative of apparent power flow w.r.t voltage:
 597 #--------------------------------------------------------------------------
 598
 599 def dAbr_dV(dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St):
 600     """ Partial derivatives of squared flow magnitudes w.r.t voltage.
 601
 602         Computes partial derivatives of apparent power w.r.t active and
 603         reactive power flows.  Partial derivative must equal 1 for lines
 604         with zero flow to avoid division by zero errors (1 comes from
 605         L'Hopital).
 606     """
 607     dAf_dPf = spdiag(2 * Sf.real())
 608     dAf_dQf = spdiag(2 * Sf.imag())
 609     dAt_dPt = spdiag(2 * St.real())
 610     dAt_dQt = spdiag(2 * St.imag())
 611
 612     # Partial derivative of apparent power magnitude w.r.t voltage
 613     # phase angle.
 614     dAf_dVa = dAf_dPf * dSf_dVa.real() + dAf_dQf * dSf_dVa.imag()
 615     dAt_dVa = dAt_dPt * dSt_dVa.real() + dAt_dQt * dSt_dVa.imag()
 616     # Partial derivative of apparent power magnitude w.r.t. voltage
 617     # amplitude.
 618     dAf_dVm = dAf_dPf * dSf_dVm.real() + dAf_dQf * dSf_dVm.imag()
 619     dAt_dVm = dAt_dPt * dSt_dVm.real() + dAt_dQt * dSt_dVm.imag()
 620
 621     return dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm
 622
 623 #--------------------------------------------------------------------------
 624 #  Second derivative of power injection w.r.t voltage:
 625 #--------------------------------------------------------------------------
 626
 627 def d2Sbus_dV2(Ybus, V, lam):
 628     """ Computes 2nd derivatives of power injection w.r.t. voltage.
 629     """
 630     n = len(V)
 631     Ibus = Ybus * V
 632     diaglam = spdiag(lam)
 633     diagV = spdiag(V)
 634
 635     A = spmatrix(mul(lam, V), range(n), range(n))
 636     B = Ybus * diagV
 637     C = A * conj(B)
 638     D = Ybus.H * diagV
 639     E = conj(diagV) * (D * diaglam - spmatrix(D*lam, range(n), range(n)))
 640     F = C - A * spmatrix(conj(Ibus), range(n), range(n))
 641     G = spmatrix(div(matrix(1.0, (n, 1)), abs(V)), range(n), range(n))
 642
 643     Gaa = E + F
 644     Gva = 1j * G * (E - F)
 645     Gav = Gva.T
 646     Gvv = G * (C + C.T) * G
 647
 648     return Gaa, Gav, Gva, Gvv
 649
 650 #--------------------------------------------------------------------------
 651 #  Second derivative of complex branch current w.r.t. voltage:
 652 #--------------------------------------------------------------------------
 653
 654 def d2Ibr_dV2(Ybr, V, lam):
 655     """ Computes 2nd derivatives of complex branch current w.r.t. voltage.
 656     """
 657     nb = len(V)
 658     diaginvVm = spdiag(div(matrix(1.0, (nb, 1)), abs(V)))
 659
 660     Haa = spdiag(mul(-(Ybr.T * lam), V))
 661     Hva = -1j * Haa * diaginvVm
 662     Hav = Hva
 663     Hvv = spmatrix([], [], [], (nb, nb))
 664
 665     return Haa, Hav, Hva, Hvv
 666
 667 #--------------------------------------------------------------------------
 668 #  Second derivative of complex power flow w.r.t. voltage:
 669 #--------------------------------------------------------------------------
 670
 671 def d2Sbr_dV2(Cbr, Ybr, V, lam):
 672     """ Computes 2nd derivatives of complex power flow w.r.t. voltage.
 673     """
 674     nb = len(V)
 675
 676     diaglam = spdiag(lam)
 677     diagV = spdiag(V)
 678
 679     A = Ybr.H * diaglam * Cbr
 680     B = conj(diagV) * A * diagV
 681     D = spdiag(mul((A*V), conj(V)))
 682     E = spdiag(mul((A.T * conj(V)), V))
 683     F = B + B.T
 684     G = spdiag(div(matrix(1.0, (nb, 1)), abs(V)))
 685
 686     Haa = F - D - E
 687     Hva = 1j * G * (B - B.T - D + E)
 688     Hav = Hva.T
 689     Hvv = G * F * G
 690
 691     return Haa, Hav, Hva, Hvv
 692
 693 #--------------------------------------------------------------------------
 694 #  Second derivative of |complex power flow|**2 w.r.t. voltage:
 695 #--------------------------------------------------------------------------
 696
 697 def d2ASbr_dV2(dSbr_dVa, dSbr_dVm, Sbr, Cbr, Ybr, V, lam):
 698     """ Computes 2nd derivatives of |complex power flow|**2 w.r.t. V.
 699     """
 700     diaglam = spdiag(lam)
 701     diagSbr_conj = spdiag(conj(Sbr))
 702
 703     Saa, Sav, Sva, Svv = d2Sbr_dV2(Cbr, Ybr, V, diagSbr_conj * lam)
 704
 705     Haa = 2 * ( Saa + dSbr_dVa.T * diaglam * conj(dSbr_dVa) ).real()
 706     Hva = 2 * ( Sva + dSbr_dVm.T * diaglam * conj(dSbr_dVa) ).real()
 707     Hav = 2 * ( Sav + dSbr_dVa.T * diaglam * conj(dSbr_dVm) ).real()
 708     Hvv = 2 * ( Svv + dSbr_dVm.T * diaglam * conj(dSbr_dVm) ).real()
 709
 710     return Haa, Hav, Hva, Hvv
 711
 712 #--------------------------------------------------------------------------
 713 #  Second derivative of |complex current|**2 w.r.t. voltage:
 714 #--------------------------------------------------------------------------
 715
 716 def d2AIbr_dV2(dIbr_dVa, dIbr_dVm, Ibr, Ybr, V, lam):
 717     """ Computes 2nd derivatives of |complex current|**2 w.r.t. V.
 718     """
 719     diaglam = spdiag(lam)
 720     diagIbr_conj = spdiag(conj(Ibr))
 721
 722     Iaa, Iav, Iva, Ivv = d2Ibr_dV2(Ybr, V, diagIbr_conj * lam)
 723
 724     Haa = 2 * ( Iaa + dIbr_dVa.T * diaglam * conj(dIbr_dVa) ).real()
 725     Hva = 2 * ( Iva + dIbr_dVm.T * diaglam * conj(dIbr_dVa) ).real()
 726     Hav = 2 * ( Iav + dIbr_dVa.T * diaglam * conj(dIbr_dVm) ).real()
 727     Hvv = 2 * ( Ivv + dIbr_dVm.T * diaglam * conj(dIbr_dVm) ).real()
 728
 729     return Haa, Hav, Hva, Hvv
 730
 731 #------------------------------------------------------------------------------
 732 # Complex conjugate:
 733 #------------------------------------------------------------------------------
 734
 735 def conj(A):
 736     """ Returns the complex conjugate of A as a new matrix.
 737     """
 738     return A.ctrans().trans()
 739
 740
 741 def tocvx(B):
 742     """ Converts a sparse SciPy matrix into a sparse CVXOPT matrix.
 743     """
 744     Bcoo = B.tocoo()
 745     return spmatrix(Bcoo.data, Bcoo.row.tolist(), Bcoo.col.tolist())
 746
 747 # EOF -------------------------------------------------------------------------