pypower.qps

38 """Quadratic Program Solver based on MOSEK. 39 40 A wrapper function providing a PYPOWER standardized interface for using 41 MOSEKOPT to solve the following QP (quadratic programming) problem:: 42 43 min 1/2 x'*H*x + c'*x 44 x 45 46 subject to:: 47 48 l <= A*x <= u (linear constraints) 49 xmin <= x <= xmax (variable bounds) 50 51 Inputs (all optional except C{H}, C{C}, C{A} and C{L}): 52 - C{H} : matrix (possibly sparse) of quadratic cost coefficients 53 - C{C} : vector of linear cost coefficients 54 - C{A, l, u} : define the optional linear constraints. Default 55 values for the elements of L and U are -Inf and Inf, respectively. 56 - xmin, xmax : optional lower and upper bounds on the 57 C{x} variables, defaults are -Inf and Inf, respectively. 58 - C{x0} : optional starting value of optimization vector C{x} 59 - C{opt} : optional options structure with the following fields, 60 all of which are also optional (default values shown in parentheses) 61 - C{verbose} (0) - controls level of progress output displayed 62 - 0 = no progress output 63 - 1 = some progress output 64 - 2 = verbose progress output 65 - C{max_it} (0) - maximum number of iterations allowed 66 - 0 = use algorithm default 67 - C{mosek_opt} - options struct for MOSEK, values in 68 C{verbose} and C{max_it} override these options 69 - C{problem} : The inputs can alternatively be supplied in a single 70 C{problem} struct with fields corresponding to the input arguments 71 described above: C{H, c, A, l, u, xmin, xmax, x0, opt} 72 73 Outputs: 74 - C{x} : solution vector 75 - C{f} : final objective function value 76 - C{exitflag} : exit flag 77 - 1 = success 78 - 0 = terminated at maximum number of iterations 79 - -1 = primal or dual infeasible 80 < 0 = the negative of the MOSEK return code 81 - C{output} : output dict with the following fields: 82 - C{r} - MOSEK return code 83 - C{res} - MOSEK result dict 84 - C{lmbda} : dict containing the Langrange and Kuhn-Tucker 85 multipliers on the constraints, with fields: 86 - C{mu_l} - lower (left-hand) limit on linear constraints 87 - C{mu_u} - upper (right-hand) limit on linear constraints 88 - C{lower} - lower bound on optimization variables 89 - C{upper} - upper bound on optimization variables 90 91 @author: Ray Zimmerman (PSERC Cornell) 92 @author: Richard Lincoln 93 """ 94 ##----- input argument handling ----- 95 ## gather inputs 96 if isinstance(H, dict): ## problem struct 97 p = H 98 else: ## individual args 99 p = {'H': H, 'c': c, 'A': A, 'l': l, 'u': u} 100 if xmin is not None: 101 p['xmin'] = xmin 102 if xmax is not None: 103 p['xmax'] = xmax 104 if x0 is not None: 105 p['x0'] = x0 106 if opt is not None: 107 p['opt'] = opt 108 109 ## define nx, set default values for H and c 110 if 'H' not in p or len(p['H']) or not any(any(p['H'])): 111 if ('A' not in p) | len(p['A']) == 0 & \ 112 ('xmin' not in p) | len(p['xmin']) == 0 & \ 113 ('xmax' not in p) | len(p['xmax']) == 0: 114 stderr.write('qps_mosek: LP problem must include constraints or variable bounds\n') 115 else: 116 if 'A' in p & len(p['A']) > 0: 117 nx = shape(p['A'])[1] 118 elif 'xmin' in p & len(p['xmin']) > 0: 119 nx = len(p['xmin']) 120 else: # if isfield(p, 'xmax') && ~isempty(p.xmax) 121 nx = len(p['xmax']) 122 p['H'] = sparse((nx, nx)) 123 qp = 0 124 else: 125 nx = shape(p['H'])[0] 126 qp = 1 127 128 if 'c' not in p | len(p['c']) == 0: 129 p['c'] = zeros(nx) 130 131 if 'x0' not in p | len(p['x0']) == 0: 132 p['x0'] = zeros(nx) 133 134 ## default options 135 if 'opt' not in p: 136 p['opt'] = [] 137 138 if 'verbose' in p['opt']: 139 verbose = p['opt']['verbose'] 140 else: 141 verbose = 0 142 143 if 'max_it' in p['opt']: 144 max_it = p['opt']['max_it'] 145 else: 146 max_it = 0 147 148 if 'mosek_opt' in p['opt']: 149 mosek_opt = mosek_options(p['opt']['mosek_opt']) 150 else: 151 mosek_opt = mosek_options() 152 153 if max_it: 154 mosek_opt['MSK_IPAR_INTPNT_MAX_ITERATIONS'] = max_it 155 156 if qp: 157 mosek_opt['MSK_IPAR_OPTIMIZER'] = 0 ## default solver only for QP 158 159 ## set up problem struct for MOSEK 160 prob = {} 161 prob['c'] = p['c'] 162 if qp: 163 prob['qosubi'], prob['qosubj'], prob['qoval'] = find(tril(sparse(p['H']))) 164 165 if 'A' in p & len(p['A']) > 0: 166 prob['a'] = sparse(p['A']) 167 168 if 'l' in p & len(p['A']) > 0: 169 prob['blc'] = p['l'] 170 171 if 'u' in p & len(p['A']) > 0: 172 prob['buc'] = p['u'] 173 174 if 'xmin' in p & len(p['xmin']) > 0: 175 prob['blx'] = p['xmin'] 176 177 if 'xmax' in p & len(p['xmax']) > 0: 178 prob['bux'] = p['xmax'] 179 180 ## A is not allowed to be empty 181 if 'a' not in prob | len(prob['a']) == 0: 182 unconstrained = True 183 prob['a'] = sparse((1, (1, 1)), (1, nx)) 184 prob.blc = -Inf 185 prob.buc = Inf 186 else: 187 unconstrained = False 188 189 ##----- run optimization ----- 190 cmd = 'minimize echo(%d)' % verbose 191 r, res = mosekopt(cmd, prob, mosek_opt) 192 193 ##----- repackage results ----- 194 if 'sol' in res: 195 if 'bas' in res['sol']: 196 sol = res['sol.bas'] 197 else: 198 sol = res['sol.itr'] 199 x = sol['xx'] 200 else: 201 sol = array([]) 202 x = array([]) 203 204 ##----- process return codes ----- 205 if 'symbcon' in res: 206 sc = res['symbcon'] 207 else: 208 r2, res2 = mosekopt('symbcon echo(0)') 209 sc = res2['symbcon'] 210 211 eflag = -r 212 msg = '' 213 if r == sc.MSK_RES_OK: 214 if len(sol) > 0: 215 # if sol['solsta'] == sc.MSK_SOL_STA_OPTIMAL: 216 if sol['solsta'] == 'OPTIMAL': 217 msg = 'The solution is optimal.' 218 eflag = 1 219 else: 220 eflag = -1 221 # if sol['prosta'] == sc['MSK_PRO_STA_PRIM_INFEAS']: 222 if sol['prosta'] == 'PRIMAL_INFEASIBLE': 223 msg = 'The problem is primal infeasible.' 224 # elif sol['prosta'] == sc['MSK_PRO_STA_DUAL_INFEAS']: 225 elif sol['prosta'] == 'DUAL_INFEASIBLE': 226 msg = 'The problem is dual infeasible.' 227 else: 228 msg = sol['solsta'] 229 230 elif r == sc['MSK_RES_TRM_MAX_ITERATIONS']: 231 eflag = 0 232 msg = 'The optimizer terminated at the maximum number of iterations.' 233 else: 234 if 'rmsg' in res and 'rcodestr' in res: 235 msg = '%s : %s' % (res['rcodestr'], res['rmsg']) 236 else: 237 msg = 'MOSEK return code = %d' % r 238 239 if verbose and len(msg) > 0: 240 stdout.write('%s\n' % msg) 241 242 ##----- repackage results ----- 243 if r == 0: 244 f = p['c'].T * x 245 if len(p['H']) > 0: 246 f = 0.5 * x.T * p['H'] * x + f 247 else: 248 f = array([]) 249 250 output = {} 251 output['r'] = r 252 output['res'] = res 253 254 if 'sol' in res: 255 lmbda = {} 256 lmbda['lower'] = sol['slx'] 257 lmbda['upper'] = sol['sux'] 258 lmbda['mu_l'] = sol['slc'] 259 lmbda['mu_u'] = sol['suc'] 260 if unconstrained: 261 lmbda['mu_l'] = array([]) 262 lmbda['mu_u'] = array([]) 263 else: 264 lmbda = array([]) 265 266 return x, f, eflag, output, lmbda

Source Code for Module pypower.qps_mosek