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