pypower.qps

40 """Quadratic Program Solver based on CPLEX. 41 42 A wrapper function providing a PYPOWER standardized interface for using 43 C{cplexqp} or C{cplexlp} to solve the following QP (quadratic programming) 44 problem:: 45 46 min 1/2 X'*H*x + c'*x 47 x 48 49 subject to:: 50 51 l <= A*x <= u (linear constraints) 52 xmin <= x <= xmax (variable bounds) 53 54 Inputs (all optional except C{H}, C{c}, C{A} and C{l}): 55 - C{H} : matrix (possibly sparse) of quadratic cost coefficients 56 - C{c} : vector of linear cost coefficients 57 - C{A, l, u} : define the optional linear constraints. Default 58 values for the elements of L and U are -Inf and Inf, respectively. 59 - C{xmin, xmax} : optional lower and upper bounds on the 60 C{x} variables, defaults are -Inf and Inf, respectively. 61 - C{x0} : optional starting value of optimization vector C{x} 62 - C{opt} : optional options structure with the following fields, 63 all of which are also optional (default values shown in parentheses) 64 - C{verbose} (0) - controls level of progress output displayed 65 - 0 = no progress output 66 - 1 = some progress output 67 - 2 = verbose progress output 68 - C{cplex_opt} - options dict for CPLEX, value in 69 verbose overrides these options 70 - C{problem} : The inputs can alternatively be supplied in a single 71 C{problem} dict with fields corresponding to the input arguments 72 described above: C{H, c, A, l, u, xmin, xmax, x0, opt} 73 74 Outputs: 75 - C{x} : solution vector 76 - C{f} : final objective function value 77 - C{exitflag} : CPLEXQP/CPLEXLP exit flag 78 (see C{cplexqp} and C{cplexlp} documentation for details) 79 - C{output} : CPLEXQP/CPLEXLP output dict 80 (see C{cplexqp} and C{cplexlp} documentation for details) 81 - C{lmbda} : dict containing the Langrange and Kuhn-Tucker 82 multipliers on the constraints, with fields: 83 - mu_l - lower (left-hand) limit on linear constraints 84 - mu_u - upper (right-hand) limit on linear constraints 85 - lower - lower bound on optimization variables 86 - upper - upper bound on optimization variables 87 88 @author: Ray Zimmerman (PSERC Cornell) 89 @author: Richard Lincoln 90 """ 91 ##----- input argument handling ----- 92 ## gather inputs 93 if isinstance(H, dict): ## problem struct 94 p = H 95 if 'opt' in p: opt = p['opt'] 96 if 'x0' in p: x0 = p['x0'] 97 if 'xmax' in p: xmax = p['xmax'] 98 if 'xmin' in p: xmin = p['xmin'] 99 if 'u' in p: u = p['u'] 100 if 'l' in p: l = p['l'] 101 if 'A' in p: A = p['A'] 102 if 'c' in p: c = p['c'] 103 if 'H' in p: H = p['H'] 104 else: ## individual args 105 assert H is not None 106 assert c is not None 107 assert A is not None 108 assert l is not None 109 110 if opt is None: 111 opt = {} 112 # if x0 is None: 113 # x0 = array([]) 114 # if xmax is None: 115 # xmax = array([]) 116 # if xmin is None: 117 # xmin = array([]) 118 119 ## define nx, set default values for missing optional inputs 120 if len(H) == 0 or not any(any(H)): 121 if len(A) == 0 and len(xmin) == 0 and len(xmax) == 0: 122 stderr.write('qps_cplex: LP problem must include constraints or variable bounds\n') 123 else: 124 if len(A) > 0: 125 nx = shape(A)[1] 126 elif len(xmin) > 0: 127 nx = len(xmin) 128 else: # if len(xmax) > 0 129 nx = len(xmax) 130 else: 131 nx = shape(H)[0] 132 133 if len(c) == 0: 134 c = zeros(nx) 135 136 if len(A) > 0 and (len(l) == 0 or all(l == -Inf)) and \ 137 (len(u) == 0 or all(u == Inf)): 138 A = None ## no limits => no linear constraints 139 140 nA = shape(A)[0] ## number of original linear constraints 141 if len(u) == 0: ## By default, linear inequalities are ... 142 u = Inf * ones(nA) ## ... unbounded above and ... 143 144 if len(l) == 0: 145 l = -Inf * ones(nA) ## ... unbounded below. 146 147 if len(xmin) == 0: ## By default, optimization variables are ... 148 xmin = -Inf * ones(nx) ## ... unbounded below and ... 149 150 if len(xmax) == 0: 151 xmax = Inf * ones(nx) ## ... unbounded above. 152 153 if len(x0) == 0: 154 x0 = zeros(nx) 155 156 ## default options 157 if 'verbose' in opt: 158 verbose = opt['verbose'] 159 else: 160 verbose = 0 161 162 #if 'max_it' in opt: 163 # max_it = opt['max_it'] 164 #else: 165 # max_it = 0 166 167 ## split up linear constraints 168 ieq = find( abs(u-l) <= EPS ) ## equality 169 igt = find( u >= 1e10 & l > -1e10 ) ## greater than, unbounded above 170 ilt = find( l <= -1e10 & u < 1e10 ) ## less than, unbounded below 171 ibx = find( (abs(u-l) > EPS) & (u < 1e10) & (l > -1e10) ) 172 Ae = A[ieq, :] 173 be = u[ieq] 174 Ai = r_[ A[ilt, :], -A[igt, :], A[ibx, :] -A[ibx, :] ] 175 bi = r_[ u[ilt], -l[igt], u[ibx], -l[ibx] ] 176 177 ## grab some dimensions 178 nlt = len(ilt) ## number of upper bounded linear inequalities 179 ngt = len(igt) ## number of lower bounded linear inequalities 180 nbx = len(ibx) ## number of doubly bounded linear inequalities 181 182 ## set up options struct for CPLEX 183 if 'cplex_opt' in opt: 184 cplex_opt = cplex_options(opt['cplex_opt']) 185 else: 186 cplex_opt = cplex_options 187 188 189 vrb = max([0, verbose - 1]) 190 cplex_opt['barrier']['display'] = vrb 191 cplex_opt['conflict']['display'] = vrb 192 cplex_opt['mip']['display'] = vrb 193 cplex_opt['sifting']['display'] = vrb 194 cplex_opt['simplex']['display'] = vrb 195 cplex_opt['tune']['display'] = vrb 196 #if max_it: 197 # cplex_opt. ## not sure what to set here 198 199 if len(Ai) == 0 and len(Ae) == 0: 200 unconstrained = 1 201 Ae = sparse((1, nx)) 202 be = 0 203 else: 204 unconstrained = 0 205 206 ## call the solver 207 if verbose: 208 cplex = Cplex('null') 209 methods = [ 210 'default', 211 'primal simplex', 212 'dual simplex', 213 'network simplex', 214 'barrier', 215 'sifting', 216 'concurrent' 217 ] 218 219 if len(H) == 0 or not any(any(H)): 220 if verbose: 221 stdout.write('CPLEX Version %s -- %s LP solver\n' % 222 (cplex.getVersion(), methods[cplex_opt['lpmethod'] + 1])) 223 224 x, f, eflag, output, lam = \ 225 cplexlp(c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt) 226 else: 227 if verbose: 228 stdout.write('CPLEX Version %s -- %s QP solver\n' % 229 (cplex.getVersion(), methods[cplex_opt['qpmethod'] + 1])) 230 231 x, f, eflag, output, lam = \ 232 cplexqp(H, c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt) 233 234 235 ## check for empty results (in case optimization failed) 236 if len(x) == 0: 237 x = NaN * zeros(nx) 238 239 if len(f) == 0: 240 f = NaN 241 242 if len(lam) == 0: 243 lam['ineqlin'] = NaN * zeros(len(bi)) 244 lam['eqlin'] = NaN * zeros(len(be)) 245 lam['lower'] = NaN * zeros(nx) 246 lam['upper'] = NaN * zeros(nx) 247 mu_l = NaN * zeros(nA) 248 mu_u = NaN * zeros(nA) 249 else: 250 mu_l = zeros(nA) 251 mu_u = zeros(nA) 252 253 if unconstrained: 254 lam['eqlin'] = array([]) 255 256 ## repackage lambdas 257 kl = find(lam.eqlin > 0) ## lower bound binding 258 ku = find(lam.eqlin < 0) ## upper bound binding 259 260 mu_l[ieq[kl]] = lam['eqlin'][kl] 261 mu_l[igt] = -lam['ineqlin'][nlt + arange(ngt)] 262 mu_l[ibx] = -lam['ineqlin'][nlt + ngt + nbx + arange(nbx)] 263 264 mu_u[ieq[ku]] = -lam['eqlin'][ku] 265 mu_u[ilt] = -lam['ineqlin'][:nlt] 266 mu_u[ibx] = -lam['ineqlin'][nlt + ngt + arange(nbx)] 267 268 lmbda = { 269 'mu_l': mu_l, 270 'mu_u': mu_u, 271 'lower': lam.lower, 272 'upper': lam.upper 273 } 274 275 return x, f, eflag, output, lmbda

Source Code for Module pypower.qps_cplex