Package pypower :: Module qps_pypower
[hide private]
[frames] | no frames]

Source Code for Module pypower.qps_pypower

  1  # Copyright (C) 2010-2011 Power System Engineering Research Center (PSERC) 
  2  # Copyright (C) 2011 Richard Lincoln 
  3  # 
  4  # PYPOWER is free software: you can redistribute it and/or modify 
  5  # it under the terms of the GNU General Public License as published 
  6  # by the Free Software Foundation, either version 3 of the License, 
  7  # or (at your option) any later version. 
  8  # 
  9  # PYPOWER 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 PYPOWER. If not, see <http://www.gnu.org/licenses/>. 
 16   
 17  """Quadratic Program Solver for PYPOWER. 
 18  """ 
 19   
 20  import sys 
 21   
 22  from pypower.qps_pips import qps_pips 
 23  from pypower.qps_ipopt import qps_ipopt 
 24  from pypower.qps_cplex import qps_cplex 
 25  from pypower.qps_mosek import qps_mosek 
 26   
 27   
28 -def qps_pypower(H, c=None, A=None, l=None, u=None, xmin=None, xmax=None, 29 x0=None, opt=None):
30 """Quadratic Program Solver for PYPOWER. 31 32 A common wrapper function for various QP solvers. 33 Solves the following QP (quadratic programming) problem:: 34 35 min 1/2 x'*H*x + c'*x 36 x 37 38 subject to:: 39 40 l <= A*x <= u (linear constraints) 41 xmin <= x <= xmax (variable bounds) 42 43 Inputs (all optional except C{H}, C{c}, C{A} and C{l}): 44 - C{H} : matrix (possibly sparse) of quadratic cost coefficients 45 - C{c} : vector of linear cost coefficients 46 - C{A, l, u} : define the optional linear constraints. Default 47 values for the elements of C{l} and C{u} are -Inf and Inf, 48 respectively. 49 - C{xmin}, C{xmax} : optional lower and upper bounds on the 50 C{x} variables, defaults are -Inf and Inf, respectively. 51 - C{x0} : optional starting value of optimization vector C{x} 52 - C{opt} : optional options structure with the following fields, 53 all of which are also optional (default values shown in parentheses) 54 - C{alg} (0) - determines which solver to use 55 - 0 = automatic, first available of BPMPD_MEX, CPLEX, MIPS 56 - 100 = BPMPD_MEX 57 - 200 = MIPS, MATLAB Interior Point Solver 58 pure MATLAB implementation of a primal-dual 59 interior point method 60 - 250 = MIPS-sc, a step controlled variant of MIPS 61 - 300 = Optimization Toolbox, QUADPROG or LINPROG 62 - 400 = IPOPT 63 - 500 = CPLEX 64 - 600 = MOSEK 65 - C{verbose} (0) - controls level of progress output displayed 66 - 0 = no progress output 67 - 1 = some progress output 68 - 2 = verbose progress output 69 - C{max_it} (0) - maximum number of iterations allowed 70 - 0 = use algorithm default 71 - C{bp_opt} - options vector for BP 72 - C{cplex_opt} - options dict for CPLEX 73 - C{ipopt_opt} - options dict for IPOPT 74 - C{pips_opt} - options dict for L{qps_pips} 75 - C{mosek_opt} - options dict for MOSEK 76 - C{ot_opt} - options dict for QUADPROG/LINPROG 77 - C{problem} : The inputs can alternatively be supplied in a single 78 C{problem} dict with fields corresponding to the input arguments 79 described above: C{H, c, A, l, u, xmin, xmax, x0, opt} 80 81 Outputs: 82 - C{x} : solution vector 83 - C{f} : final objective function value 84 - C{exitflag} : exit flag 85 - 1 = converged 86 - 0 or negative values = algorithm specific failure codes 87 - C{output} : output struct with the following fields: 88 - C{alg} - algorithm code of solver used 89 - (others) - algorithm specific fields 90 - C{lmbda} : dict containing the Langrange and Kuhn-Tucker 91 multipliers on the constraints, with fields: 92 - C{mu_l} - lower (left-hand) limit on linear constraints 93 - C{mu_u} - upper (right-hand) limit on linear constraints 94 - C{lower} - lower bound on optimization variables 95 - C{upper} - upper bound on optimization variables 96 97 98 Example from U{http://www.uc.edu/sashtml/iml/chap8/sect12.htm}: 99 100 >>> from numpy import array, zeros, Inf 101 >>> from scipy.sparse import csr_matrix 102 >>> H = csr_matrix(array([[1003.1, 4.3, 6.3, 5.9], 103 ... [4.3, 2.2, 2.1, 3.9], 104 ... [6.3, 2.1, 3.5, 4.8], 105 ... [5.9, 3.9, 4.8, 10 ]])) 106 >>> c = zeros(4) 107 >>> A = csr_matrix(array([[1, 1, 1, 1 ], 108 ... [0.17, 0.11, 0.10, 0.18]])) 109 >>> l = array([1, 0.10]) 110 >>> u = array([1, Inf]) 111 >>> xmin = zeros(4) 112 >>> xmax = None 113 >>> x0 = array([1, 0, 0, 1]) 114 >>> solution = qps_pips(H, c, A, l, u, xmin, xmax, x0) 115 >>> round(solution["f"], 11) == 1.09666678128 116 True 117 >>> solution["converged"] 118 True 119 >>> solution["output"]["iterations"] 120 10 121 122 @author: Ray Zimmerman (PSERC Cornell) 123 @author: Richard Lincoln 124 """ 125 ##----- input argument handling ----- 126 ## gather inputs 127 if isinstance(H, dict): ## problem struct 128 p = H 129 if 'opt' in p: opt = p['opt'] 130 if 'x0' in p: x0 = p['x0'] 131 if 'xmax' in p: xmax = p['xmax'] 132 if 'xmin' in p: xmin = p['xmin'] 133 if 'u' in p: u = p['u'] 134 if 'l' in p: l = p['l'] 135 if 'A' in p: A = p['A'] 136 if 'c' in p: c = p['c'] 137 if 'H' in p: H = p['H'] 138 else: ## individual args 139 # assert H is not None zero dimensional sparse matrices not supported 140 assert c is not None 141 # assert A is not None zero dimensional sparse matrices not supported 142 # assert l is not None no lower bounds indicated by None 143 144 if opt is None: 145 opt = {} 146 # if x0 is None: 147 # x0 = array([]) 148 # if xmax is None: 149 # xmax = array([]) 150 # if xmin is None: 151 # xmin = array([]) 152 153 ## default options 154 if 'alg' in opt: 155 alg = opt['alg'] 156 else: 157 alg = 0 158 159 if 'verbose' in opt: 160 verbose = opt['verbose'] 161 else: 162 verbose = 0 163 164 if alg == 0: 165 try: 166 import pymosek #@UnusedImport 167 alg = 600 ## use MOSEK by default if available 168 except ImportError: 169 try: 170 import cplex #@UnusedImport 171 alg = 500 ## if not, then CPLEX if available 172 except ImportError: 173 alg = 200 ## otherwise PIPS 174 175 176 ##----- call the appropriate solver ----- 177 if alg == 200 or alg == 250: ## use MIPS or sc-MIPS 178 ## set up options 179 if 'pips_opt' in opt: 180 pips_opt = opt['pips_opt'] 181 else: 182 pips_opt = {} 183 184 if 'max_it' in opt: 185 pips_opt['max_it'] = opt['max_it'] 186 187 if alg == 200: 188 pips_opt['step_control'] = False 189 else: 190 pips_opt['step_control'] = True 191 192 pips_opt['verbose'] = verbose 193 194 ## call solver 195 x, f, eflag, output, lmbda = \ 196 qps_pips(H, c, A, l, u, xmin, xmax, x0, pips_opt) 197 elif alg == 400: ## use IPOPT 198 x, f, eflag, output, lmbda = \ 199 qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt) 200 elif alg == 500: ## use CPLEX 201 x, f, eflag, output, lmbda = \ 202 qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt) 203 elif alg == 600: ## use MOSEK 204 x, f, eflag, output, lmbda = \ 205 qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt) 206 else: 207 sys.stderr.write('qps_pypower: %d is not a valid algorithm code\n', alg) 208 209 if 'alg' not in output: 210 output['alg'] = alg 211 212 return x, f, eflag, output, lmbda
213 214 215 if __name__ == "__main__": 216 import doctest 217 doctest.testmod() 218