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 % Copyright (c) 2010-2015 by Power System Engineering Research Center (PSERC) 0097 % by Ray Zimmerman, PSERC Cornell 0098 % 0099 % $Id: qps_mosek.m 2644 2015-03-11 19:34:22Z ray $ 0100 % 0101 % This file is part of MATPOWER. 0102 % Covered by the 3-clause BSD License (see LICENSE file for details). 0103 % See http://www.pserc.cornell.edu/matpower/ for more info. 0104 0105 %% check for Optimization Toolbox 0106 % if ~have_fcn('mosek') 0107 % error('qps_mosek: requires MOSEK'); 0108 % end 0109 0110 %%----- input argument handling ----- 0111 %% gather inputs 0112 if nargin == 1 && isstruct(H) %% problem struct 0113 p = H; 0114 else %% individual args 0115 p = struct('H', H, 'c', c, 'A', A, 'l', l, 'u', u); 0116 if nargin > 5 0117 p.xmin = xmin; 0118 if nargin > 6 0119 p.xmax = xmax; 0120 if nargin > 7 0121 p.x0 = x0; 0122 if nargin > 8 0123 p.opt = opt; 0124 end 0125 end 0126 end 0127 end 0128 end 0129 0130 %% define nx, set default values for H and c 0131 if ~isfield(p, 'H') || isempty(p.H) || ~any(any(p.H)) 0132 if (~isfield(p, 'A') || isempty(p.A)) && ... 0133 (~isfield(p, 'xmin') || isempty(p.xmin)) && ... 0134 (~isfield(p, 'xmax') || isempty(p.xmax)) 0135 error('qps_mosek: LP problem must include constraints or variable bounds'); 0136 else 0137 if isfield(p, 'A') && ~isempty(p.A) 0138 nx = size(p.A, 2); 0139 elseif isfield(p, 'xmin') && ~isempty(p.xmin) 0140 nx = length(p.xmin); 0141 else % if isfield(p, 'xmax') && ~isempty(p.xmax) 0142 nx = length(p.xmax); 0143 end 0144 end 0145 p.H = sparse(nx, nx); 0146 qp = 0; 0147 else 0148 nx = size(p.H, 1); 0149 qp = 1; 0150 end 0151 if ~isfield(p, 'c') || isempty(p.c) 0152 p.c = zeros(nx, 1); 0153 end 0154 if ~isfield(p, 'x0') || isempty(p.x0) 0155 p.x0 = zeros(nx, 1); 0156 end 0157 0158 %% default options 0159 if ~isfield(p, 'opt') 0160 p.opt = []; 0161 end 0162 if ~isempty(p.opt) && isfield(p.opt, 'verbose') && ~isempty(p.opt.verbose) 0163 verbose = p.opt.verbose; 0164 else 0165 verbose = 0; 0166 end 0167 if ~isempty(p.opt) && isfield(p.opt, 'mosek_opt') && ~isempty(p.opt.mosek_opt) 0168 mosek_opt = mosek_options(p.opt.mosek_opt); 0169 else 0170 mosek_opt = mosek_options; 0171 end 0172 if qp 0173 mosek_opt.MSK_IPAR_OPTIMIZER = 0; %% default solver only for QP 0174 end 0175 0176 %% set up problem struct for MOSEK 0177 prob.c = p.c; 0178 if qp 0179 [prob.qosubi, prob.qosubj, prob.qoval] = find(tril(sparse(p.H))); 0180 end 0181 if isfield(p, 'A') && ~isempty(p.A) 0182 prob.a = sparse(p.A); 0183 nA = size(p.A, 1); 0184 else 0185 nA = 0; 0186 end 0187 if isfield(p, 'l') && ~isempty(p.A) 0188 prob.blc = p.l; 0189 end 0190 if isfield(p, 'u') && ~isempty(p.A) 0191 prob.buc = p.u; 0192 end 0193 if isfield(p, 'xmin') && ~isempty(p.xmin) 0194 prob.blx = p.xmin; 0195 end 0196 if isfield(p, 'xmax') && ~isempty(p.xmax) 0197 prob.bux = p.xmax; 0198 end 0199 0200 %% A is not allowed to be empty 0201 if ~isfield(prob, 'a') || isempty(prob.a) 0202 unconstrained = 1; 0203 prob.a = sparse(1, 1, 1, 1, nx); 0204 prob.blc = -Inf; 0205 prob.buc = Inf; 0206 else 0207 unconstrained = 0; 0208 end 0209 0210 %%----- run optimization ----- 0211 if verbose 0212 s = have_fcn('mosek', 'all'); 0213 if s.vnum < 7 0214 methods = { 0215 'default', %% 0 : MSK_OPTIMIZER_FREE 0216 'interior point', %% 1 : MSK_OPTIMIZER_INTPNT 0217 '<conic>', %% 2 : MSK_OPTIMIZER_CONIC 0218 '<qcone>', %% 3 : MSK_OPTIMIZER_QCONE 0219 'primal simplex', %% 4 : MSK_OPTIMIZER_PRIMAL_SIMPLEX 0220 'dual simplex', %% 5 : MSK_OPTIMIZER_DUAL_SIMPLEX 0221 'primal dual simplex', %% 6 : MSK_OPTIMIZER_PRIMAL_DUAL_SIMPLEX 0222 'automatic simplex', %% 7 : MSK_OPTIMIZER_FREE_SIMPLEX 0223 '<mixed int>', %% 8 : MSK_OPTIMIZER_MIXED_INT 0224 '<nonconvex>', %% 9 : MSK_OPTIMIZER_NONCONVEX 0225 'concurrent' %% 10 : MSK_OPTIMIZER_CONCURRENT 0226 }; 0227 else 0228 methods = { 0229 'default', %% 0 : MSK_OPTIMIZER_FREE 0230 'interior point', %% 1 : MSK_OPTIMIZER_INTPNT 0231 '<conic>', %% 2 : MSK_OPTIMIZER_CONIC 0232 'primal simplex', %% 3 : MSK_OPTIMIZER_PRIMAL_SIMPLEX 0233 'dual simplex', %% 4 : MSK_OPTIMIZER_DUAL_SIMPLEX 0234 'primal dual simplex', %% 5 : MSK_OPTIMIZER_PRIMAL_DUAL_SIMPLEX 0235 'automatic simplex', %% 6 : MSK_OPTIMIZER_FREE_SIMPLEX 0236 'network simplex', %% 7 : MSK_OPTIMIZER_NETWORK_PRIMAL_SIMPLEX 0237 '<mixed int conic>', %% 8 : MSK_OPTIMIZER_MIXED_INT_CONIC 0238 '<mixed int>', %% 9 : MSK_OPTIMIZER_MIXED_INT 0239 'concurrent', %% 10 : MSK_OPTIMIZER_CONCURRENT 0240 '<nonconvex>' %% 11 : MSK_OPTIMIZER_NONCONVEX 0241 }; 0242 end 0243 if qp 0244 lpqp = 'QP'; 0245 else 0246 lpqp = 'LP'; 0247 end 0248 vn = have_fcn('mosek', 'vstr'); 0249 if isempty(vn) 0250 vn = '<unknown>'; 0251 end 0252 fprintf('MOSEK Version %s -- %s %s solver\n', ... 0253 vn, methods{mosek_opt.MSK_IPAR_OPTIMIZER+1}, lpqp); 0254 end 0255 cmd = sprintf('minimize echo(%d)', verbose); 0256 [r, res] = mosekopt(cmd, prob, mosek_opt); 0257 0258 %%----- repackage results ----- 0259 if isfield(res, 'sol') 0260 if isfield(res.sol, 'bas') 0261 sol = res.sol.bas; 0262 else 0263 sol = res.sol.itr; 0264 end 0265 x = sol.xx; 0266 else 0267 sol = []; 0268 x = NaN(nx, 1); 0269 end 0270 0271 %%----- process return codes ----- 0272 if isfield(res, 'symbcon') 0273 sc = res.symbcon; 0274 else 0275 sc = mosek_symbcon; 0276 end 0277 eflag = -r; 0278 msg = ''; 0279 switch (r) 0280 case sc.MSK_RES_OK 0281 if ~isempty(sol) 0282 % if sol.solsta == sc.MSK_SOL_STA_OPTIMAL 0283 if strcmp(sol.solsta, 'OPTIMAL') 0284 msg = 'The solution is optimal.'; 0285 eflag = 1; 0286 else 0287 eflag = -1; 0288 % if sol.prosta == sc.MSK_PRO_STA_PRIM_INFEAS 0289 if strcmp(sol.prosta, 'PRIMAL_INFEASIBLE') 0290 msg = 'The problem is primal infeasible.'; 0291 % elseif sol.prosta == sc.MSK_PRO_STA_DUAL_INFEAS 0292 elseif strcmp(sol.prosta, 'DUAL_INFEASIBLE') 0293 msg = 'The problem is dual infeasible.'; 0294 else 0295 msg = sol.solsta; 0296 end 0297 end 0298 end 0299 case sc.MSK_RES_TRM_MAX_ITERATIONS 0300 eflag = 0; 0301 msg = 'The optimizer terminated at the maximum number of iterations.'; 0302 otherwise 0303 if isfield(res, 'rmsg') && isfield(res, 'rcodestr') 0304 msg = sprintf('%s : %s', res.rcodestr, res.rmsg); 0305 else 0306 msg = sprintf('MOSEK return code = %d', r); 0307 end 0308 end 0309 0310 if (verbose || r == sc.MSK_RES_ERR_LICENSE || ... 0311 r == sc.MSK_RES_ERR_LICENSE_EXPIRED || ... 0312 r == sc.MSK_RES_ERR_LICENSE_VERSION || ... 0313 r == sc.MSK_RES_ERR_LICENSE_NO_SERVER_SUPPORT || ... 0314 r == sc.MSK_RES_ERR_LICENSE_FEATURE || ... 0315 r == sc.MSK_RES_ERR_LICENSE_INVALID_HOSTID || ... 0316 r == sc.MSK_RES_ERR_LICENSE_SERVER_VERSION || ... 0317 r == sc.MSK_RES_ERR_MISSING_LICENSE_FILE) ... 0318 && ~isempty(msg) %% always alert user of license problems 0319 fprintf('%s\n', msg); 0320 end 0321 0322 %%----- repackage results ----- 0323 if nargout > 1 0324 if r == 0 0325 f = p.c' * x; 0326 if ~isempty(p.H) 0327 f = 0.5 * x' * p.H * x + f; 0328 end 0329 else 0330 f = []; 0331 end 0332 if nargout > 3 0333 output.r = r; 0334 output.res = res; 0335 if nargout > 4 0336 if ~isempty(sol) 0337 if isfield(sol, 'slx') 0338 lambda.lower = sol.slx; 0339 else 0340 lambda.lower = []; 0341 end 0342 if isfield(sol, 'sux') 0343 lambda.upper = sol.sux; 0344 else 0345 lambda.upper = []; 0346 end 0347 if isfield(sol, 'slc') 0348 lambda.mu_l = sol.slc; 0349 else 0350 lambda.mu_l = []; 0351 end 0352 if isfield(sol, 'suc') 0353 lambda.mu_u = sol.suc; 0354 else 0355 lambda.mu_u = []; 0356 end 0357 else 0358 if isfield(p, 'xmin') && ~isempty(p.xmin) 0359 lambda.lower = NaN(nx, 1); 0360 else 0361 lambda.lower = []; 0362 end 0363 if isfield(p, 'xmax') && ~isempty(p.xmax) 0364 lambda.upper = NaN(nx, 1); 0365 else 0366 lambda.upper = []; 0367 end 0368 lambda.mu_l = NaN(nA, 1); 0369 lambda.mu_u = NaN(nA, 1); 0370 end 0371 if unconstrained 0372 lambda.mu_l = []; 0373 lambda.mu_u = []; 0374 end 0375 end 0376 end 0377 end