pylon/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 IPOPT OPF solver.
  18 """
  19
  20 #------------------------------------------------------------------------------
  21 #  Imports:
  22 #------------------------------------------------------------------------------
  23
  24 import pyipopt # http://github.com/rwl/pyipopt
  25
  26 from numpy import Inf, ones, r_, zeros
  27
  28 from scipy.sparse import vstack, tril
  29
  30 from pylon.solver import PIPSSolver
  31
  32 #------------------------------------------------------------------------------
  33 #  "IPOPFSolver" class:
  34 #------------------------------------------------------------------------------
  35
  36 class IPOPFSolver(PIPSSolver):
  37     """ Solves AC optimal power flow using IPOPT.
  38     """
  39
  40     #--------------------------------------------------------------------------
  41     #  PIPSSolver interface:
  42     #--------------------------------------------------------------------------
  43
  44     def _solve(self, x0, A, l, u, xmin, xmax):
  45         """ Solves using the Interior Point OPTimizer.
  46         """
  47         # Indexes of constrained lines.
  48         il = [i for i,ln in enumerate(self._ln) if 0.0 < ln.rate_a < 1e10]
  49         nl2 = len(il)
  50
  51         neqnln = 2 * self._nb # no. of non-linear equality constraints
  52         niqnln = 2 * len(il)  # no. of lines with constraints
  53
  54         user_data = {"A": A, "neqnln": neqnln, "niqnln": niqnln}
  55
  56         self._f(x0)
  57         Jdata = self._dg(x0, False, user_data)
  58 #        Hdata = self._h(x0, ones(neqnln + niqnln), None, False, user_data)
  59
  60         lmbda = {"eqnonlin": ones(neqnln),
  61                  "ineqnonlin": ones(niqnln)}
  62         H = tril(self._hessfcn(x0, lmbda), format="coo")
  63         self._Hrow, self._Hcol = H.row, H.col
  64
  65         n = len(x0) # the number of variables
  66         xl = xmin
  67         xu = xmax
  68         gl = r_[zeros(2 * self._nb), -Inf * ones(2 * nl2), l]
  69         gu = r_[zeros(2 * self._nb),       zeros(2 * nl2), u]
  70         m = len(gl) # the number of constraints
  71         nnzj = len(Jdata) # the number of nonzeros in Jacobian matrix
  72         nnzh = 0#len(H.data) # the number of non-zeros in Hessian matrix
  73
  74         f_fcn, df_fcn, g_fcn, dg_fcn, h_fcn = \
  75             self._f, self._df, self._g, self._dg, self._h
  76
  77         nlp = pyipopt.create(n, xl, xu, m, gl, gu, nnzj, nnzh,
  78                              f_fcn, df_fcn, g_fcn, dg_fcn)#, h_fcn)
  79
  80 #        print dir(nlp)
  81 #        nlp.str_option("print_options_documentation", "yes")
  82 #        nlp.int_option("max_iter", 10)
  83
  84 #        x, zl, zu, obj = nlp.solve(x0)
  85         success = nlp.solve(x0, user_data)
  86         nlp.close()
  87
  88     #--------------------------------------------------------------------------
  89     #  IPOPFSolver interface:
  90     #--------------------------------------------------------------------------
  91
  92     def _g(self, x, user_data):
  93         A = user_data["A"]
  94         h, g = self._gh(x)
  95         if A is None:
  96             b = r_[g, h]
  97         else:
  98             b = r_[g, h, A * x]
  99         return b
 100
 101
 102     def _dg(self, x, flag, user_data):
 103         A = user_data["A"]
 104         dh, dg = self._dgh(x)
 105         if A is None:
 106             J = vstack([dg.T, dh.T], "coo")
 107         else:
 108             J = vstack([dg.T, dh.T, A], "coo")
 109         if flag:
 110 #            return (J.row, J.col)
 111             return (J.col, J.row)
 112         else:
 113             return J.data
 114
 115
 116     def _h(self, x, lagrange, obj_factor, flag, user_data=None):
 117         if flag:
 118 #            return (self._Hrow, self._Hcol)
 119             return (self._Hcol, self._Hrow)
 120         else:
 121             neqnln = user_data["neqnln"]
 122             niqnln = user_data["niqnln"]
 123             lmbda = {"eqnonlin": lagrange[:neqnln],
 124                      "ineqnonlin": lagrange[neqnln:neqnln + niqnln]}
 125             H = tril(self._hessfcn(x, lmbda), format="coo")
 126             return H.data
 127
 128 # EOF -------------------------------------------------------------------------