QPS_MOSEK Quadratic Program Solver based on MOSEK.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
QPS_MOSEK(H, C, A, L, U, XMIN, XMAX, X0, OPT)
A wrapper function providing a MATPOWER standardized interface for using
MOSEKOPT to solve the following QP (quadratic programming) problem:
min 1/2 X'*H*X + C'*X
X
subject to
L <= A*X <= U (linear constraints)
XMIN <= X <= XMAX (variable bounds)
Inputs (all optional except H, C, A and L):
H : matrix (possibly sparse) of quadratic cost coefficients
C : vector of linear cost coefficients
A, L, U : define the optional linear constraints. Default
values for the elements of L and U are -Inf and Inf,
respectively.
XMIN, XMAX : optional lower and upper bounds on the
X variables, defaults are -Inf and Inf, respectively.
X0 : optional starting value of optimization vector X
OPT : optional options structure with the following fields,
all of which are also optional (default values shown in
parentheses)
verbose (0) - controls level of progress output displayed
0 = no progress output
1 = some progress output
2 = verbose progress output
max_it (0) - maximum number of iterations allowed
0 = use algorithm default
mosek_opt - options struct for MOSEK, values in
verbose and max_it override these options
PROBLEM : The inputs can alternatively be supplied in a single
PROBLEM struct with fields corresponding to the input arguments
described above: H, c, A, l, u, xmin, xmax, x0, opt
Outputs:
X : solution vector
F : final objective function value
EXITFLAG : exit flag
1 = success
0 = terminated at maximum number of iterations
-1 = primal or dual infeasible
< 0 = the negative of the MOSEK return code
OUTPUT : output struct with the following fields:
r - MOSEK return code
res - MOSEK result struct
LAMBDA : struct containing the Langrange and Kuhn-Tucker
multipliers on the constraints, with fields:
mu_l - lower (left-hand) limit on linear constraints
mu_u - upper (right-hand) limit on linear constraints
lower - lower bound on optimization variables
upper - upper bound on optimization variables
Note the calling syntax is almost identical to that of QUADPROG
from MathWorks' Optimization Toolbox. The main difference is that
the linear constraints are specified with A, L, U instead of
A, B, Aeq, Beq.
Calling syntax options:
[x, f, exitflag, output, lambda] = ...
qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_mosek(H, c, A, l, u)
x = qps_mosek(H, c, A, l, u, xmin, xmax)
x = qps_mosek(H, c, A, l, u, xmin, xmax, x0)
x = qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_mosek(problem), where problem is a struct with fields:
H, c, A, l, u, xmin, xmax, x0, opt
all fields except 'c', 'A' and 'l' or 'u' are optional
x = qps_mosek(...)
[x, f] = qps_mosek(...)
[x, f, exitflag] = qps_mosek(...)
[x, f, exitflag, output] = qps_mosek(...)
[x, f, exitflag, output, lambda] = qps_mosek(...)
Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm)
H = [ 1003.1 4.3 6.3 5.9;
4.3 2.2 2.1 3.9;
6.3 2.1 3.5 4.8;
5.9 3.9 4.8 10 ];
c = zeros(4,1);
A = [ 1 1 1 1;
0.17 0.11 0.10 0.18 ];
l = [1; 0.10];
u = [1; Inf];
xmin = zeros(4,1);
x0 = [1; 0; 0; 1];
opt = struct('verbose', 2);
[x, f, s, out, lambda] = qps_mosek(H, c, A, l, u, xmin, [], x0, opt);
See also MOSEKOPT.0001 function [x, f, eflag, output, lambda] = qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt) 0002 %QPS_MOSEK Quadratic Program Solver based on MOSEK. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % QPS_MOSEK(H, C, A, L, U, XMIN, XMAX, X0, OPT) 0005 % A wrapper function providing a MATPOWER standardized interface for using 0006 % MOSEKOPT to solve the following QP (quadratic programming) problem: 0007 % 0008 % min 1/2 X'*H*X + C'*X 0009 % X 0010 % 0011 % subject to 0012 % 0013 % L <= A*X <= U (linear constraints) 0014 % XMIN <= X <= XMAX (variable bounds) 0015 % 0016 % Inputs (all optional except H, C, A and L): 0017 % H : matrix (possibly sparse) of quadratic cost coefficients 0018 % C : vector of linear cost coefficients 0019 % A, L, U : define the optional linear constraints. Default 0020 % values for the elements of L and U are -Inf and Inf, 0021 % respectively. 0022 % XMIN, XMAX : optional lower and upper bounds on the 0023 % X variables, defaults are -Inf and Inf, respectively. 0024 % X0 : optional starting value of optimization vector X 0025 % OPT : optional options structure with the following fields, 0026 % all of which are also optional (default values shown in 0027 % parentheses) 0028 % verbose (0) - controls level of progress output displayed 0029 % 0 = no progress output 0030 % 1 = some progress output 0031 % 2 = verbose progress output 0032 % max_it (0) - maximum number of iterations allowed 0033 % 0 = use algorithm default 0034 % mosek_opt - options struct for MOSEK, values in 0035 % verbose and max_it override these options 0036 % PROBLEM : The inputs can alternatively be supplied in a single 0037 % PROBLEM struct with fields corresponding to the input arguments 0038 % described above: H, c, A, l, u, xmin, xmax, x0, opt 0039 % 0040 % Outputs: 0041 % X : solution vector 0042 % F : final objective function value 0043 % EXITFLAG : exit flag 0044 % 1 = success 0045 % 0 = terminated at maximum number of iterations 0046 % -1 = primal or dual infeasible 0047 % < 0 = the negative of the MOSEK return code 0048 % OUTPUT : output struct with the following fields: 0049 % r - MOSEK return code 0050 % res - MOSEK result struct 0051 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0052 % multipliers on the constraints, with fields: 0053 % mu_l - lower (left-hand) limit on linear constraints 0054 % mu_u - upper (right-hand) limit on linear constraints 0055 % lower - lower bound on optimization variables 0056 % upper - upper bound on optimization variables 0057 % 0058 % Note the calling syntax is almost identical to that of QUADPROG 0059 % from MathWorks' Optimization Toolbox. The main difference is that 0060 % the linear constraints are specified with A, L, U instead of 0061 % A, B, Aeq, Beq. 0062 % 0063 % Calling syntax options: 0064 % [x, f, exitflag, output, lambda] = ... 0065 % qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt) 0066 % 0067 % x = qps_mosek(H, c, A, l, u) 0068 % x = qps_mosek(H, c, A, l, u, xmin, xmax) 0069 % x = qps_mosek(H, c, A, l, u, xmin, xmax, x0) 0070 % x = qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt) 0071 % x = qps_mosek(problem), where problem is a struct with fields: 0072 % H, c, A, l, u, xmin, xmax, x0, opt 0073 % all fields except 'c', 'A' and 'l' or 'u' are optional 0074 % x = qps_mosek(...) 0075 % [x, f] = qps_mosek(...) 0076 % [x, f, exitflag] = qps_mosek(...) 0077 % [x, f, exitflag, output] = qps_mosek(...) 0078 % [x, f, exitflag, output, lambda] = qps_mosek(...) 0079 % 0080 % Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm) 0081 % H = [ 1003.1 4.3 6.3 5.9; 0082 % 4.3 2.2 2.1 3.9; 0083 % 6.3 2.1 3.5 4.8; 0084 % 5.9 3.9 4.8 10 ]; 0085 % c = zeros(4,1); 0086 % A = [ 1 1 1 1; 0087 % 0.17 0.11 0.10 0.18 ]; 0088 % l = [1; 0.10]; 0089 % u = [1; Inf]; 0090 % xmin = zeros(4,1); 0091 % x0 = [1; 0; 0; 1]; 0092 % opt = struct('verbose', 2); 0093 % [x, f, s, out, lambda] = qps_mosek(H, c, A, l, u, xmin, [], x0, opt); 0094 % 0095 % See also MOSEKOPT. 0096 0097 % MATPOWER 0098 % $Id: qps_mosek.m,v 1.5 2011/09/09 15:26:08 cvs Exp $ 0099 % by Ray Zimmerman, PSERC Cornell 0100 % Copyright (c) 2010 by Power System Engineering Research Center (PSERC) 0101 % 0102 % This file is part of MATPOWER. 0103 % See http://www.pserc.cornell.edu/matpower/ for more info. 0104 % 0105 % MATPOWER is free software: you can redistribute it and/or modify 0106 % it under the terms of the GNU General Public License as published 0107 % by the Free Software Foundation, either version 3 of the License, 0108 % or (at your option) any later version. 0109 % 0110 % MATPOWER is distributed in the hope that it will be useful, 0111 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0112 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0113 % GNU General Public License for more details. 0114 % 0115 % You should have received a copy of the GNU General Public License 0116 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0117 % 0118 % Additional permission under GNU GPL version 3 section 7 0119 % 0120 % If you modify MATPOWER, or any covered work, to interface with 0121 % other modules (such as MATLAB code and MEX-files) available in a 0122 % MATLAB(R) or comparable environment containing parts covered 0123 % under other licensing terms, the licensors of MATPOWER grant 0124 % you additional permission to convey the resulting work. 0125 0126 %% check for Optimization Toolbox 0127 % if ~have_fcn('mosek') 0128 % error('qps_mosek: requires MOSEK'); 0129 % end 0130 0131 %%----- input argument handling ----- 0132 %% gather inputs 0133 if nargin == 1 && isstruct(H) %% problem struct 0134 p = H; 0135 else %% individual args 0136 p = struct('H', H, 'c', c, 'A', A, 'l', l, 'u', u); 0137 if nargin > 5 0138 p.xmin = xmin; 0139 if nargin > 6 0140 p.xmax = xmax; 0141 if nargin > 7 0142 p.x0 = x0; 0143 if nargin > 8 0144 p.opt = opt; 0145 end 0146 end 0147 end 0148 end 0149 end 0150 0151 %% define nx, set default values for H and c 0152 if ~isfield(p, 'H') || isempty(p.H) || ~any(any(p.H)) 0153 if (~isfield(p, 'A') || isempty(p.A)) && ... 0154 (~isfield(p, 'xmin') || isempty(p.xmin)) && ... 0155 (~isfield(p, 'xmax') || isempty(p.xmax)) 0156 error('qps_mosek: LP problem must include constraints or variable bounds'); 0157 else 0158 if isfield(p, 'A') && ~isempty(p.A) 0159 nx = size(p.A, 2); 0160 elseif isfield(p, 'xmin') && ~isempty(p.xmin) 0161 nx = length(p.xmin); 0162 else % if isfield(p, 'xmax') && ~isempty(p.xmax) 0163 nx = length(p.xmax); 0164 end 0165 end 0166 p.H = sparse(nx, nx); 0167 qp = 0; 0168 else 0169 nx = size(p.H, 1); 0170 qp = 1; 0171 end 0172 if ~isfield(p, 'c') || isempty(p.c) 0173 p.c = zeros(nx, 1); 0174 end 0175 if ~isfield(p, 'x0') || isempty(p.x0) 0176 p.x0 = zeros(nx, 1); 0177 end 0178 0179 %% default options 0180 if ~isfield(p, 'opt') 0181 p.opt = []; 0182 end 0183 if ~isempty(p.opt) && isfield(p.opt, 'verbose') && ~isempty(p.opt.verbose) 0184 verbose = p.opt.verbose; 0185 else 0186 verbose = 0; 0187 end 0188 if ~isempty(p.opt) && isfield(p.opt, 'max_it') && ~isempty(p.opt.max_it) 0189 max_it = p.opt.max_it; 0190 else 0191 max_it = 0; 0192 end 0193 if ~isempty(p.opt) && isfield(p.opt, 'mosek_opt') && ~isempty(p.opt.mosek_opt) 0194 mosek_opt = mosek_options(p.opt.mosek_opt); 0195 else 0196 mosek_opt = mosek_options; 0197 end 0198 if max_it 0199 mosek_opt.MSK_IPAR_INTPNT_MAX_ITERATIONS = max_it; 0200 end 0201 if qp 0202 mosek_opt.MSK_IPAR_OPTIMIZER = 0; %% default solver only for QP 0203 end 0204 0205 %% set up problem struct for MOSEK 0206 prob.c = p.c; 0207 if qp 0208 [prob.qosubi, prob.qosubj, prob.qoval] = find(tril(sparse(p.H))); 0209 end 0210 if isfield(p, 'A') && ~isempty(p.A) 0211 prob.a = sparse(p.A); 0212 end 0213 if isfield(p, 'l') && ~isempty(p.A) 0214 prob.blc = p.l; 0215 end 0216 if isfield(p, 'u') && ~isempty(p.A) 0217 prob.buc = p.u; 0218 end 0219 if isfield(p, 'xmin') && ~isempty(p.xmin) 0220 prob.blx = p.xmin; 0221 end 0222 if isfield(p, 'xmax') && ~isempty(p.xmax) 0223 prob.bux = p.xmax; 0224 end 0225 0226 %% A is not allowed to be empty 0227 if ~isfield(prob, 'a') || isempty(prob.a) 0228 unconstrained = 1; 0229 prob.a = sparse(1, 1, 1, 1, nx); 0230 prob.blc = -Inf; 0231 prob.buc = Inf; 0232 else 0233 unconstrained = 0; 0234 end 0235 0236 %%----- run optimization ----- 0237 if verbose 0238 methods = { 0239 'default', 0240 'interior point', 0241 '<default>', 0242 '<default>', 0243 'primal simplex', 0244 'dual simplex', 0245 'primal dual simplex', 0246 'automatic simplex', 0247 '<default>', 0248 '<default>', 0249 'concurrent' 0250 }; 0251 if isempty(H) || ~any(any(H)) 0252 lpqp = 'LP'; 0253 else 0254 lpqp = 'QP'; 0255 end 0256 % (this code is also in mpver.m) 0257 % MOSEK Version 6.0.0.93 (Build date: 2010-10-26 13:03:27) 0258 % MOSEK Version 6.0.0.106 (Build date: 2011-3-17 10:46:54) 0259 % pat = 'Version (\.*\d)+.*Build date: (\d\d\d\d-\d\d-\d\d)'; 0260 pat = 'Version (\.*\d)+.*Build date: (\d+-\d+-\d+)'; 0261 [s,e,tE,m,t] = regexp(evalc('mosekopt'), pat); 0262 if isempty(t) 0263 vn = '<unknown>'; 0264 else 0265 vn = t{1}{1}; 0266 end 0267 fprintf('MOSEK Version %s -- %s %s solver\n', ... 0268 vn, methods{mosek_opt.MSK_IPAR_OPTIMIZER+1}, lpqp); 0269 end 0270 cmd = sprintf('minimize echo(%d)', verbose); 0271 [r, res] = mosekopt(cmd, prob, mosek_opt); 0272 0273 %%----- repackage results ----- 0274 if isfield(res, 'sol') 0275 if isfield(res.sol, 'bas') 0276 sol = res.sol.bas; 0277 else 0278 sol = res.sol.itr; 0279 end 0280 x = sol.xx; 0281 else 0282 sol = []; 0283 x = []; 0284 end 0285 0286 %%----- process return codes ----- 0287 if isfield(res, 'symbcon') 0288 sc = res.symbcon; 0289 else 0290 [r2, res2] = mosekopt('symbcon echo(0)'); 0291 sc = res2.symbcon; 0292 end 0293 eflag = -r; 0294 msg = ''; 0295 switch (r) 0296 case sc.MSK_RES_OK 0297 if ~isempty(sol) 0298 % if sol.solsta == sc.MSK_SOL_STA_OPTIMAL 0299 if strcmp(sol.solsta, 'OPTIMAL') 0300 msg = 'The solution is optimal.'; 0301 eflag = 1; 0302 else 0303 eflag = -1; 0304 % if sol.prosta == sc.MSK_PRO_STA_PRIM_INFEAS 0305 if strcmp(sol.prosta, 'PRIMAL_INFEASIBLE') 0306 msg = 'The problem is primal infeasible.'; 0307 % elseif sol.prosta == sc.MSK_PRO_STA_DUAL_INFEAS 0308 elseif strcmp(sol.prosta, 'DUAL_INFEASIBLE') 0309 msg = 'The problem is dual infeasible.'; 0310 else 0311 msg = sol.solsta; 0312 end 0313 end 0314 end 0315 case sc.MSK_RES_TRM_MAX_ITERATIONS 0316 eflag = 0; 0317 msg = 'The optimizer terminated at the maximum number of iterations.'; 0318 otherwise 0319 if isfield(res, 'rmsg') && isfield(res, 'rcodestr') 0320 msg = sprintf('%s : %s', res.rcodestr, res.rmsg); 0321 else 0322 msg = sprintf('MOSEK return code = %d', r); 0323 end 0324 end 0325 0326 if (verbose || r == 1001) && ~isempty(msg) %% always alert user if license is expired 0327 fprintf('%s\n', msg); 0328 end 0329 0330 %%----- repackage results ----- 0331 if nargout > 1 0332 if r == 0 0333 f = p.c' * x; 0334 if ~isempty(p.H) 0335 f = 0.5 * x' * p.H * x + f; 0336 end 0337 else 0338 f = []; 0339 end 0340 if nargout > 3 0341 output.r = r; 0342 output.res = res; 0343 if nargout > 4 0344 if isfield(res, 'sol') 0345 lambda.lower = sol.slx; 0346 lambda.upper = sol.sux; 0347 lambda.mu_l = sol.slc; 0348 lambda.mu_u = sol.suc; 0349 if unconstrained 0350 lambda.mu_l = []; 0351 lambda.mu_u = []; 0352 end 0353 else 0354 lambda = []; 0355 end 0356 end 0357 end 0358 end