contrib/ipopf.py

   1 #------------------------------------------------------------------------------
   2 # Copyright (C) 2007-2010 Richard Lincoln
   3 #
   4 # Licensed under the Apache License, Version 2.0 (the "License");
   5 # you may not use this file except in compliance with the License.
   6 # You may obtain a copy of the License at
   7 #
   8 #     http://www.apache.org/licenses/LICENSE-2.0
   9 #
  10 # Unless required by applicable law or agreed to in writing, software
  11 # distributed under the License is distributed on an "AS IS" BASIS,
  12 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  13 # See the License for the specific language governing permissions and
  14 # limitations under the License.
  15 #------------------------------------------------------------------------------
  16
  17 """ Defines an OPF solver for Pylon using IPOPT.
  18 """
  19
  20 #------------------------------------------------------------------------------
  21 #  Imports:
  22 #------------------------------------------------------------------------------
  23
  24 import pyipopt
  25
  26 from numpy import \
  27     array, polyder, polyval, exp, conj, Inf, ones, r_, zeros, asarray
  28
  29 from scipy.sparse import lil_matrix, csr_matrix, hstack, vstack
  30
  31 from pylon import REFERENCE
  32 from pylon.solver import _Solver, SFLOW, IFLOW, PFLOW
  33
  34 #------------------------------------------------------------------------------
  35 #  "IPOPFSolver" class:
  36 #------------------------------------------------------------------------------
  37
  38 class IPOPFSolver(_Solver):
  39     """ Solves AC optimal power flow using IPOPT.
  40     """
  41
  42     def __init__(self, om, flow_lim=SFLOW, opt=None):
  43         """ Initialises a new IPOPFSolver instance.
  44         """
  45         super(IPOPFSolver, self).__init__(om)
  46
  47         #: Quantity to limit for branch flow constraints ("S", "P" or "I").
  48         self.flow_lim = flow_lim
  49
  50         #: Options for the PIPS.
  51         self.opt = {} if opt is None else opt
  52
  53
  54     def _ref_bus_angle_constraint(self, buses, Va, xmin, xmax):
  55         """ Adds a constraint on the reference bus angles.
  56         """
  57         refs = [bus._i for bus in buses if bus.type == REFERENCE]
  58         Varefs = array([b.v_angle for b in buses if b.type == REFERENCE])
  59
  60         xmin[Va.i1 - 1 + refs] = Varefs
  61         xmax[Va.iN - 1 + refs] = Varefs
  62
  63         return xmin, xmax
  64
  65
  66     def solve(self):
  67         """ Solves AC optimal power flow.
  68         """
  69         case = self.om.case
  70         base_mva = case.base_mva
  71         # TODO: Explain this value.
  72         self.opt["cost_mult"] = 1e-4
  73
  74         # Unpack the OPF model.
  75         bs, ln, gn, _ = self._unpack_model(self.om)
  76         # Compute problem dimensions.
  77         ipol, _, nb, nl, _, ny, nxyz = self._dimension_data(bs, ln, gn)
  78
  79         # Compute problem dimensions.
  80         ng = len(gn)
  81 #        gpol = [g for g in gn if g.pcost_model == POLYNOMIAL]
  82         # Indexes of constrained lines.
  83         il = array([i for i,l in enumerate(ln) if 0.0 < l.rate_a < 1e10])
  84         nl2 = len(il)
  85
  86         # Linear constraints (l <= A*x <= u).
  87         A, l, u = self.om.linear_constraints()
  88 #        AA, bb = self._linear_constraints(self.om)
  89
  90         _, xmin, xmax = self._var_bounds()
  91
  92         # Select an interior initial point for interior point solver.
  93         x0 = self._initial_interior_point(bs, gn, xmin, xmax, ny)
  94
  95         # Build admittance matrices.
  96         Ybus, Yf, Yt = case.Y
  97
  98         # Optimisation variables.
  99         Va = self.om.get_var("Va")
 100         Vm = self.om.get_var("Vm")
 101         Pg = self.om.get_var("Pg")
 102         Qg = self.om.get_var("Qg")
 103
 104         # Adds a constraint on the reference bus angles.
 105 #        xmin, xmax = self._ref_bus_angle_constraint(bs, Va, xmin, xmax)
 106
 107         def f_fcn(x, user_data=None):
 108             """ Evaluates the objective function.
 109             """
 110             p_gen = x[Pg.i1:Pg.iN + 1] # Active generation in p.u.
 111             q_gen = x[Qg.i1:Qg.iN + 1] # Reactive generation in p.u.
 112
 113             # Polynomial cost of P and Q.
 114             xx = r_[p_gen, q_gen] * base_mva
 115             if len(ipol) > 0:
 116                 f = sum([g.total_cost(xx[i]) for i,g in enumerate(gn)])
 117             else:
 118                 f = 0
 119
 120             # Piecewise linear cost of P and Q.
 121             if ny:
 122                 y = self.om.get_var("y")
 123                 ccost = csr_matrix((ones(ny),
 124                     (range(y.i1, y.iN + 1), zeros(ny))), shape=(nxyz, 1)).T
 125                 f = f + ccost * x
 126             else:
 127                 ccost = zeros((1, nxyz))
 128                 # TODO: Generalised cost term.
 129
 130             return f
 131
 132
 133         def df_fcn(x, usr_data=None):
 134             """ Calculates gradient of the objective function.
 135             """
 136             p_gen = x[Pg.i1:Pg.iN + 1] # Active generation in p.u.
 137             q_gen = x[Qg.i1:Qg.iN + 1] # Reactive generation in p.u.
 138
 139             xx = r_[p_gen, q_gen] * base_mva
 140
 141             if ny > 0:
 142                 y = self.om.get_var("y")
 143                 iy = range(y.i1, y.iN + 1)
 144                 ccost = \
 145                     csr_matrix((ones(ny), (iy, zeros(ny))), shape=(nxyz, 1)).T
 146             else:
 147                 ccost = zeros((1, nxyz))
 148                 # TODO: Generalised cost term.
 149
 150             iPg = range(Pg.i1, Pg.iN + 1)
 151             iQg = range(Qg.i1, Qg.iN + 1)
 152
 153             # Polynomial cost of P and Q.
 154             df_dPgQg = zeros((2 * ng, 1))        # w.r.t p.u. Pg and Qg
 155 #            df_dPgQg[ipol] = matrix([g.poly_cost(xx[i], 1) for g in gpol])
 156 #            for i, g in enumerate(gn):
 157 #                der = polyder(list(g.p_cost))
 158 #                df_dPgQg[i] = polyval(der, xx[i]) * base_mva
 159             for i in ipol:
 160                 df_dPgQg[i] = \
 161                     base_mva * polyval(polyder(list(gn[i].p_cost)), xx[i])
 162
 163             df = zeros((nxyz, 1))
 164             df[iPg] = df_dPgQg[:ng]
 165             df[iQg] = df_dPgQg[ng:ng + ng]
 166
 167             # Piecewise linear cost of P and Q.
 168             df = df + ccost.T
 169             # TODO: Generalised cost term.
 170
 171             return asarray(df).flatten()
 172
 173
 174         def g_fcn(x, usr_data=None):
 175             """ Evaluates the non-linear constraint values.
 176             """
 177             Pgen = x[Pg.i1:Pg.iN + 1] # Active generation in p.u.
 178             Qgen = x[Qg.i1:Qg.iN + 1] # Reactive generation in p.u.
 179
 180             for i, g in enumerate(gn):
 181                 g.p = Pgen[i] * base_mva # active generation in MW
 182                 g.q = Qgen[i] * base_mva # reactive generation in MVAr
 183
 184             # Rebuild the net complex bus power injection vector in p.u.
 185             Sbus = case.getSbus(bs)
 186
 187             Vang = x[Va.i1:Va.iN + 1]
 188             Vmag = x[Vm.i1:Vm.iN + 1]
 189             V = Vmag * exp(1j * Vang)
 190
 191             # Evaluate the power flow equations.
 192             mis = V * conj(Ybus * V) - Sbus
 193
 194             # Equality constraints (power flow).
 195             g = r_[mis.real,  # active power mismatch for all buses
 196                    mis.imag]  # reactive power mismatch for all buses
 197
 198             # Inequality constraints (branch flow limits).
 199             # (line constraint is actually on square of limit)
 200             flow_max = array([(l.rate_a / base_mva)**2 for l in ln])
 201             # FIXME: There must be a more elegant method for this.
 202             for i, v in enumerate(flow_max):
 203                 if v == 0.0:
 204                     flow_max[i] = Inf
 205
 206             if self.flow_lim == IFLOW:
 207                 If = Yf * V
 208                 It = Yt * V
 209                 # Branch current limits.
 210                 h = r_[(If * conj(If)) - flow_max,
 211                        (If * conj(It)) - flow_max]
 212             else:
 213                 i_fbus = [e.from_bus._i for e in ln]
 214                 i_tbus = [e.to_bus._i for e in ln]
 215                 # Complex power injected at "from" bus (p.u.).
 216                 Sf = V[i_fbus] * conj(Yf * V)
 217                 # Complex power injected at "to" bus (p.u.).
 218                 St = V[i_tbus] * conj(Yt * V)
 219                 if self.flow_lim == PFLOW: # active power limit, P (Pan Wei)
 220                     # Branch real power limits.
 221                     h = r_[Sf.real()**2 - flow_max,
 222                            St.real()**2 - flow_max]
 223                 elif self.flow_lim == SFLOW: # apparent power limit, |S|
 224                     # Branch apparent power limits.
 225                     h = r_[(Sf * conj(Sf)) - flow_max,
 226                            (St * conj(St)) - flow_max].real
 227                 else:
 228                     raise ValueError
 229
 230             return r_[g, h]
 231
 232
 233         def dg_fcn(x, flag, usr_data=None):
 234             """ Calculates the Jacobian matrix. It takes two arguments, the
 235                 first is the variable x and the second is a Boolean flag. If
 236                 the flag is true, the function returns a tuple of arrays
 237                 (row, col) to indicate the sparse structure of the Jacobian
 238                 matrix. If the flag is false the function returns the values
 239                 of the Jacobian matrix with length nnzj.
 240             """
 241             iVa = range(Va.i1, Va.iN + 1)
 242             iVm = range(Vm.i1, Vm.iN + 1)
 243             iPg = range(Pg.i1, Pg.iN + 1)
 244             iQg = range(Qg.i1, Qg.iN + 1)
 245             iVaVmPgQg = r_[iVa, iVm, iPg, iQg].T
 246
 247             Vang = x[Va.i1:Va.iN + 1]
 248             Vmag = x[Vm.i1:Vm.iN + 1]
 249             V = Vmag * exp(1j * Vang)
 250
 251             # Compute partials of injected bus powers.
 252             dSbus_dVm, dSbus_dVa = case.dSbus_dV(Ybus, V)
 253
 254             i_gbus = [gen.bus._i for gen in gn]
 255             neg_Cg = csr_matrix((-ones(ng), (i_gbus, range(ng))), (nb, ng))
 256
 257             # Transposed Jacobian of the power balance equality constraints.
 258             dg = lil_matrix((nxyz, 2 * nb))
 259
 260             blank = csr_matrix((nb, ng))
 261             dg[iVaVmPgQg, :] = vstack([
 262                 hstack([dSbus_dVa.real, dSbus_dVm.real, neg_Cg, blank]),
 263                 hstack([dSbus_dVa.imag, dSbus_dVm.imag, blank, neg_Cg])
 264             ], "csr").T
 265
 266             # Compute partials of flows w.r.t V.
 267             if self.flow_lim == IFLOW:
 268                 dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft = \
 269                     case.dIbr_dV(Yf, Yt, V)
 270             else:
 271                 dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft = \
 272                     case.dSbr_dV(Yf, Yt, V, bs, ln)
 273             if self.flow_lim == PFLOW:
 274                 dFf_dVa = dFf_dVa.real
 275                 dFf_dVm = dFf_dVm.real
 276                 dFt_dVa = dFt_dVa.real
 277                 dFt_dVm = dFt_dVm.real
 278                 Ff = Ff.real
 279                 Ft = Ft.real
 280
 281             # Squared magnitude of flow (complex power, current or real power).
 282             df_dVa, df_dVm, dt_dVa, dt_dVm = \
 283                 case.dAbr_dV(dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft)
 284
 285             # Construct Jacobian of inequality constraints (branch limits) and
 286             # transpose it.
 287             dh = lil_matrix((nxyz, 2 * nl))
 288             dh[r_[iVa, iVm].T, :] = vstack([hstack([df_dVa, df_dVm]),
 289                                             hstack([dt_dVa, dt_dVm])], "csr").T
 290
 291             J = vstack([dg, dh, A]).tocoo()
 292
 293             if flag:
 294                 return (J.row, J.col)
 295             else:
 296                 return J.data
 297
 298
 299         def h_fcn(x, lagrange, obj_factor, flag, usr_data=None):
 300             """ Evaluates the Hessian of the Lagrangian.
 301             """
 302             neqnln = 2 * nb
 303             niqnln = 2 * len(il) # no. of lines with constraints
 304
 305             Pgen = x[Pg.i1:Pg.iN + 1] # Active generation in p.u.
 306             Qgen = x[Qg.i1:Qg.iN + 1] # Reactive generation in p.u.
 307
 308             for i, g in enumerate(gn):
 309                 g.p = Pgen[i] * base_mva # active generation in MW
 310                 g.q = Qgen[i] * base_mva # reactive generation in MVAr
 311
 312             Vang = x[Va.i1:Va.iN + 1]
 313             Vmag = x[Vm.i1:Vm.iN + 1]
 314             V = Vmag * exp(1j * Vang)
 315             nxtra = nxyz - 2 * nb
 316
 317             #------------------------------------------------------------------
 318             #  Evaluate d2f.
 319             #------------------------------------------------------------------
 320
 321             d2f_dPg2 = lil_matrix((ng, 1)) # w.r.t p.u. Pg
 322             d2f_dQg2 = lil_matrix((ng, 1)) # w.r.t p.u. Qg]
 323
 324             for i in ipol:
 325                 d2f_dPg2[i, 0] = polyval(polyder(list(gn[i].p_cost), 2),
 326                                          Pg.v0[i] * base_mva) * base_mva**2
 327 #            for i in ipol:
 328 #                d2f_dQg2[i] = polyval(polyder(list(gn[i].p_cost), 2),
 329 #                                      Qg.v0[i] * base_mva) * base_mva**2
 330
 331             i = r_[range(Pg.i1, Pg.iN + 1), range(Qg.i1, Qg.iN + 1)]
 332
 333             d2f = csr_matrix((vstack([d2f_dPg2, d2f_dQg2]).toarray().flatten(),
 334                               (i, i)), shape=(nxyz, nxyz))
 335             # TODO: Generalised cost model.
 336             d2f = d2f * self.opt["cost_mult"]
 337
 338             #------------------------------------------------------------------
 339             #  Evaluate Hessian of power balance constraints.
 340             #------------------------------------------------------------------
 341
 342             eqnonlin = lagrange[:neqnln]
 343 #            nlam = len(lagrange["eqnonlin"]) / 2
 344             nlam = len(eqnonlin) / 2
 345             lamP = eqnonlin[:nlam]
 346             lamQ = eqnonlin[nlam:nlam + nlam]
 347             Gpaa, Gpav, Gpva, Gpvv = case.d2Sbus_dV2(Ybus, V, lamP)
 348             Gqaa, Gqav, Gqva, Gqvv = case.d2Sbus_dV2(Ybus, V, lamQ)
 349
 350             d2G = vstack([
 351                 hstack([
 352                     vstack([hstack([Gpaa, Gpav]),
 353                             hstack([Gpva, Gpvv])]).real +
 354                     vstack([hstack([Gqaa, Gqav]),
 355                             hstack([Gqva, Gqvv])]).imag,
 356                     csr_matrix((2 * nb, nxtra))]),
 357                 hstack([
 358                     csr_matrix((nxtra, 2 * nb)),
 359                     csr_matrix((nxtra, nxtra))
 360                 ])
 361             ], "csr")
 362
 363             #------------------------------------------------------------------
 364             #  Evaluate Hessian of flow constraints.
 365             #------------------------------------------------------------------
 366
 367             ineqnonlin = lagrange[neqnln:neqnln + niqnln]
 368             nmu = len(ineqnonlin) / 2
 369             muF = ineqnonlin[:nmu]
 370             muT = ineqnonlin[nmu:nmu + nmu]
 371             if self.flow_lim == "I":
 372                 dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It = \
 373                     case.dIbr_dV(Yf, Yt, V)
 374                 Hfaa, Hfav, Hfva, Hfvv = \
 375                     case.d2AIbr_dV2(dIf_dVa, dIf_dVm, If, Yf, V, muF)
 376                 Htaa, Htav, Htva, Htvv = \
 377                     case.d2AIbr_dV2(dIt_dVa, dIt_dVm, It, Yt, V, muT)
 378             else:
 379                 f = [e.from_bus._i for e in ln]
 380                 t = [e.to_bus._i for e in ln]
 381                 # Line-bus connection matrices.
 382                 Cf = csr_matrix((ones(nl), (range(nl), f)), (nl, nb))
 383                 Ct = csr_matrix((ones(nl), (range(nl), t)), (nl, nb))
 384                 dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St = \
 385                     case.dSbr_dV(Yf, Yt, V)
 386                 if self.flow_lim == PFLOW:
 387                     Hfaa, Hfav, Hfva, Hfvv = \
 388                         case.d2ASbr_dV2(dSf_dVa.real(), dSf_dVm.real(),
 389                                         Sf.real(), Cf, Yf, V, muF)
 390                     Htaa, Htav, Htva, Htvv = \
 391                         case.d2ASbr_dV2(dSt_dVa.real(), dSt_dVm.real(),
 392                                         St.real(), Ct, Yt, V, muT)
 393                 elif self.flow_lim == SFLOW:
 394                     Hfaa, Hfav, Hfva, Hfvv = \
 395                         case.d2ASbr_dV2(dSf_dVa, dSf_dVm, Sf, Cf, Yf, V, muF)
 396                     Htaa, Htav, Htva, Htvv = \
 397                         case.d2ASbr_dV2(dSt_dVa, dSt_dVm, St, Ct, Yt, V, muT)
 398                 else:
 399                     raise ValueError
 400
 401             d2H = vstack([
 402                 hstack([
 403                     vstack([hstack([Hfaa, Hfav]),
 404                             hstack([Hfva, Hfvv])]) +
 405                     vstack([hstack([Htaa, Htav]),
 406                             hstack([Htva, Htvv])]),
 407                     csr_matrix((2 * nb, nxtra))
 408                 ]),
 409                 hstack([
 410                     csr_matrix((nxtra, 2 * nb)),
 411                     csr_matrix((nxtra, nxtra))
 412                 ])
 413             ], "csr")
 414
 415             H = d2f + d2G + d2H
 416
 417             if flag:
 418                 return (H.row, H.col)
 419             else:
 420                 return H.data
 421
 422         n = len(x0) # the number of variables
 423         gl = r_[zeros(2 * nb), -Inf * ones(2 * nl2), l]
 424         gu = r_[zeros(2 * nb),       zeros(2 * nl2), u]
 425         m = len(gl) # the number of constraints
 426         nnzj = 0 # the number of nonzeros in Jacobian matrix
 427         nnzh = 0 # the number of non-zeros in Hessian matrix
 428
 429         nlp = pyipopt.create(n, xmin, xmax, m, gl, gu, nnzj, nnzh,
 430                              f_fcn, df_fcn, g_fcn, dg_fcn, h_fcn)
 431
 432 #        x, zl, zu, obj = nlp.solve(x0)
 433         success = nlp.solve(x0)
 434         nlp.close()
 435
 436         print "Success:", success
 437         print "Solution of the primal variables, x"
 438 #        print x
 439         print "Solution of the bound multipliers, z_L and z_U"
 440 #        print zl, zu
 441         print "Objective value"
 442 #        print "f(x*) =", obj
 443
 444
 445 if __name__ == "__main__":
 446     import os
 447     import pylon
 448     c = pylon.Case.load(os.path.join(os.path.dirname(pylon.__file__),
 449                                      "test", "data", "case6ww.pkl"))
 450     s = pylon.OPF(c, dc=False).solve(IPOPFSolver)
 451
 452 # EOF -------------------------------------------------------------------------