pypower.qps

37 """Quadratic Program Solver based on IPOPT. 38 39 Uses IPOPT to solve the following QP (quadratic programming) problem:: 40 41 min 1/2 x'*H*x + c'*x 42 x 43 44 subject to:: 45 46 l <= A*x <= u (linear constraints) 47 xmin <= x <= xmax (variable bounds) 48 49 Inputs (all optional except C{H}, C{C}, C{A} and C{L}): 50 - C{H} : matrix (possibly sparse) of quadratic cost coefficients 51 - C{C} : vector of linear cost coefficients 52 - C{A, l, u} : define the optional linear constraints. Default 53 values for the elements of C{l} and C{u} are -Inf and Inf, 54 respectively. 55 - C{xmin, xmax} : optional lower and upper bounds on the 56 C{x} variables, defaults are -Inf and Inf, respectively. 57 - C{x0} : optional starting value of optimization vector C{x} 58 - C{opt} : optional options structure with the following fields, 59 all of which are also optional (default values shown in parentheses) 60 - C{verbose} (0) - controls level of progress output displayed 61 - 0 = no progress output 62 - 1 = some progress output 63 - 2 = verbose progress output 64 - C{max_it} (0) - maximum number of iterations allowed 65 - 0 = use algorithm default 66 - C{ipopt_opt} - options struct for IPOPT, values in 67 C{verbose} and C{max_it} override these options 68 - C{problem} : The inputs can alternatively be supplied in a single 69 C{problem} dict with fields corresponding to the input arguments 70 described above: C{H, c, A, l, u, xmin, xmax, x0, opt} 71 72 Outputs: 73 - C{x} : solution vector 74 - C{f} : final objective function value 75 - C{exitflag} : exit flag 76 - 1 = first order optimality conditions satisfied 77 - 0 = maximum number of iterations reached 78 - -1 = numerically failed 79 - C{output} : output struct with the following fields: 80 - C{iterations} - number of iterations performed 81 - C{hist} - dict list with trajectories of the following: 82 C{feascond}, C{gradcond}, C{compcond}, C{costcond}, C{gamma}, 83 C{stepsize}, C{obj}, C{alphap}, C{alphad} 84 - message - exit message 85 - C{lmbda} : dict containing the Langrange and Kuhn-Tucker 86 multipliers on the constraints, with fields: 87 - C{mu_l} - lower (left-hand) limit on linear constraints 88 - C{mu_u} - upper (right-hand) limit on linear constraints 89 - C{lower} - lower bound on optimization variables 90 - C{upper} - upper bound on optimization variables 91 92 Calling syntax options:: 93 x, f, exitflag, output, lmbda = \ 94 qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt) 95 96 x = qps_ipopt(H, c, A, l, u) 97 x = qps_ipopt(H, c, A, l, u, xmin, xmax) 98 x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0) 99 x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt) 100 x = qps_ipopt(problem), where problem is a struct with fields: 101 H, c, A, l, u, xmin, xmax, x0, opt 102 all fields except 'c', 'A' and 'l' or 'u' are optional 103 x = qps_ipopt(...) 104 x, f = qps_ipopt(...) 105 x, f, exitflag = qps_ipopt(...) 106 x, f, exitflag, output = qps_ipopt(...) 107 x, f, exitflag, output, lmbda = qps_ipopt(...) 108 109 Example:: 110 H = [ 1003.1 4.3 6.3 5.9; 111 4.3 2.2 2.1 3.9; 112 6.3 2.1 3.5 4.8; 113 5.9 3.9 4.8 10 ] 114 c = zeros((4, 1)) 115 A = [ 1 1 1 1 116 0.17 0.11 0.10 0.18 ] 117 l = [1, 0.10] 118 u = [1, Inf] 119 xmin = zeros((4, 1)) 120 x0 = [1, 0, 0, 1] 121 opt = {'verbose': 2) 122 x, f, s, out, lam = qps_ipopt(H, c, A, l, u, xmin, [], x0, opt) 123 124 Problem from U{http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm} 125 126 @see: C{pyipopt}, L{ipopt_options} 127 128 @author: Ray Zimmerman (PSERC Cornell) 129 @author: Richard Lincoln 130 """ 131 ##----- input argument handling ----- 132 ## gather inputs 133 if isinstance(H, dict): ## problem struct 134 p = H 135 if 'opt' in p: opt = p['opt'] 136 if 'x0' in p: x0 = p['x0'] 137 if 'xmax' in p: xmax = p['xmax'] 138 if 'xmin' in p: xmin = p['xmin'] 139 if 'u' in p: u = p['u'] 140 if 'l' in p: l = p['l'] 141 if 'A' in p: A = p['A'] 142 if 'c' in p: c = p['c'] 143 if 'H' in p: H = p['H'] 144 else: ## individual args 145 assert H is not None 146 assert c is not None 147 assert A is not None 148 assert l is not None 149 150 if opt is None: 151 opt = {} 152 # if x0 is None: 153 # x0 = array([]) 154 # if xmax is None: 155 # xmax = array([]) 156 # if xmin is None: 157 # xmin = array([]) 158 159 ## define nx, set default values for missing optional inputs 160 if len(H) == 0 or not any(any(H)): 161 if len(A) == 0 and len(xmin) == 0 and len(xmax) == 0: 162 stderr.write('qps_ipopt: LP problem must include constraints or variable bounds\n') 163 else: 164 if len(A) > 0: 165 nx = shape(A)[1] 166 elif len(xmin) > 0: 167 nx = len(xmin) 168 else: # if len(xmax) > 0 169 nx = len(xmax) 170 H = sparse((nx, nx)) 171 else: 172 nx = shape(H)[0] 173 174 if len(c) == 0: 175 c = zeros(nx) 176 177 if len(A) > 0 and (len(l) == 0 or all(l == -Inf)) and \ 178 (len(u) == 0 or all(u == Inf)): 179 A = None ## no limits => no linear constraints 180 181 nA = shape(A)[0] ## number of original linear constraints 182 if nA: 183 if len(u) == 0: ## By default, linear inequalities are ... 184 u = Inf * ones(nA) ## ... unbounded above and ... 185 186 if len(l) == 0: 187 l = -Inf * ones(nA) ## ... unbounded below. 188 189 if len(x0) == 0: 190 x0 = zeros(nx) 191 192 ## default options 193 if 'verbose' in opt: 194 verbose = opt['verbose'] 195 else: 196 verbose = 0 197 198 if 'max_it' in opt: 199 max_it = opt['max_it'] 200 else: 201 max_it = 0 202 203 ## make sure args are sparse/full as expected by IPOPT 204 if len(H) > 0: 205 if not issparse(H): 206 H = sparse(H) 207 208 if not issparse(A): 209 A = sparse(A) 210 211 ##----- run optimization ----- 212 ## set options dict for IPOPT 213 options = {} 214 if 'ipopt_opt' in opt: 215 options['ipopt'] = ipopt_options(opt['ipopt_opt']) 216 else: 217 options['ipopt'] = ipopt_options() 218 219 options['ipopt']['jac_c_constant'] = 'yes' 220 options['ipopt']['jac_d_constant'] = 'yes' 221 options['ipopt']['hessian_constant'] = 'yes' 222 options['ipopt']['least_square_init_primal'] = 'yes' 223 options['ipopt']['least_square_init_duals'] = 'yes' 224 # options['ipopt']['mehrotra_algorithm'] = 'yes' ## default 'no' 225 if verbose: 226 options['ipopt']['print_level'] = min(12, verbose * 2 + 1) 227 else: 228 options['ipopt']['print_level = 0'] 229 230 if max_it: 231 options['ipopt']['max_iter'] = max_it 232 233 ## define variable and constraint bounds, if given 234 if nA: 235 options['cu'] = u 236 options['cl'] = l 237 238 if len(xmin) > 0: 239 options['lb'] = xmin 240 241 if len(xmax) > 0: 242 options['ub'] = xmax 243 244 ## assign function handles 245 funcs = {} 246 funcs['objective'] = lambda x: 0.5 * x.T * H * x + c.T * x 247 funcs['gradient'] = lambda x: H * x + c 248 funcs['constraints'] = lambda x: A * x 249 funcs['jacobian'] = lambda x: A 250 funcs['jacobianstructure'] = lambda : A 251 funcs['hessian'] = lambda x, sigma, lmbda: tril(H) 252 funcs['hessianstructure'] = lambda : tril(H) 253 254 ## run the optimization 255 x, info = pyipopt(x0, funcs, options) 256 257 if info['status'] == 0 | info['status'] == 1: 258 eflag = 1 259 else: 260 eflag = 0 261 262 output = {} 263 if 'iter' in info: 264 output['iterations'] = info['iter'] 265 266 output['info'] = info['status'] 267 f = funcs['objective'](x) 268 269 ## repackage lmbdas 270 kl = find(info['lmbda'] < 0) ## lower bound binding 271 ku = find(info['lmbda'] > 0) ## upper bound binding 272 mu_l = zeros(nA) 273 mu_l[kl] = -info['lmbda'][kl] 274 mu_u = zeros(nA) 275 mu_u[ku] = info['lmbda'][ku] 276 277 lmbda = { 278 'mu_l': mu_l, 279 'mu_u': mu_u, 280 'lower': info['zl'], 281 'upper': info['zu'] 282 } 283 284 return x, f, eflag, output, lmbda

Source Code for Module pypower.qps_ipopt