QPS_MOSEK Quadratic Program Solver based on MOSEK.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
QPS_MOSEK(H, C, A, L, U, XMIN, XMAX, X0, OPT)
[X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_MOSEK(PROBLEM)
A wrapper function providing a 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 https://v8doc.sas.com/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 QPS_MASTER, 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 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_MOSEK(PROBLEM) 0006 % A wrapper function providing a standardized interface for using 0007 % MOSEKOPT to solve the following QP (quadratic programming) problem: 0008 % 0009 % min 1/2 X'*H*X + C'*X 0010 % X 0011 % 0012 % subject to 0013 % 0014 % L <= A*X <= U (linear constraints) 0015 % XMIN <= X <= XMAX (variable bounds) 0016 % 0017 % Inputs (all optional except H, C, A and L): 0018 % H : matrix (possibly sparse) of quadratic cost coefficients 0019 % C : vector of linear cost coefficients 0020 % A, L, U : define the optional linear constraints. Default 0021 % values for the elements of L and U are -Inf and Inf, 0022 % respectively. 0023 % XMIN, XMAX : optional lower and upper bounds on the 0024 % X variables, defaults are -Inf and Inf, respectively. 0025 % X0 : optional starting value of optimization vector X 0026 % OPT : optional options structure with the following fields, 0027 % all of which are also optional (default values shown in 0028 % parentheses) 0029 % verbose (0) - controls level of progress output displayed 0030 % 0 = no progress output 0031 % 1 = some progress output 0032 % 2 = verbose progress output 0033 % mosek_opt - options struct for MOSEK, value in verbose 0034 % overrides these options 0035 % PROBLEM : The inputs can alternatively be supplied in a single 0036 % PROBLEM struct with fields corresponding to the input arguments 0037 % described above: H, c, A, l, u, xmin, xmax, x0, opt 0038 % 0039 % Outputs: 0040 % X : solution vector 0041 % F : final objective function value 0042 % EXITFLAG : exit flag 0043 % 1 = success 0044 % 0 = terminated at maximum number of iterations 0045 % -1 = primal or dual infeasible 0046 % < 0 = the negative of the MOSEK return code 0047 % OUTPUT : output struct with the following fields: 0048 % r - MOSEK return code 0049 % res - MOSEK result struct 0050 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0051 % multipliers on the constraints, with fields: 0052 % mu_l - lower (left-hand) limit on linear constraints 0053 % mu_u - upper (right-hand) limit on linear constraints 0054 % lower - lower bound on optimization variables 0055 % upper - upper bound on optimization variables 0056 % 0057 % Note the calling syntax is almost identical to that of QUADPROG 0058 % from MathWorks' Optimization Toolbox. The main difference is that 0059 % the linear constraints are specified with A, L, U instead of 0060 % A, B, Aeq, Beq. 0061 % 0062 % Calling syntax options: 0063 % [x, f, exitflag, output, lambda] = ... 0064 % qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt) 0065 % 0066 % x = qps_mosek(H, c, A, l, u) 0067 % x = qps_mosek(H, c, A, l, u, xmin, xmax) 0068 % x = qps_mosek(H, c, A, l, u, xmin, xmax, x0) 0069 % x = qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt) 0070 % x = qps_mosek(problem), where problem is a struct with fields: 0071 % H, c, A, l, u, xmin, xmax, x0, opt 0072 % all fields except 'c', 'A' and 'l' or 'u' are optional 0073 % x = qps_mosek(...) 0074 % [x, f] = qps_mosek(...) 0075 % [x, f, exitflag] = qps_mosek(...) 0076 % [x, f, exitflag, output] = qps_mosek(...) 0077 % [x, f, exitflag, output, lambda] = qps_mosek(...) 0078 % 0079 % Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm) 0080 % H = [ 1003.1 4.3 6.3 5.9; 0081 % 4.3 2.2 2.1 3.9; 0082 % 6.3 2.1 3.5 4.8; 0083 % 5.9 3.9 4.8 10 ]; 0084 % c = zeros(4,1); 0085 % A = [ 1 1 1 1; 0086 % 0.17 0.11 0.10 0.18 ]; 0087 % l = [1; 0.10]; 0088 % u = [1; Inf]; 0089 % xmin = zeros(4,1); 0090 % x0 = [1; 0; 0; 1]; 0091 % opt = struct('verbose', 2); 0092 % [x, f, s, out, lambda] = qps_mosek(H, c, A, l, u, xmin, [], x0, opt); 0093 % 0094 % See also QPS_MASTER, MOSEKOPT. 0095 0096 % MP-Opt-Model 0097 % Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC) 0098 % by Ray Zimmerman, PSERC Cornell 0099 % 0100 % This file is part of MP-Opt-Model. 0101 % Covered by the 3-clause BSD License (see LICENSE file for details). 0102 % See https://github.com/MATPOWER/mp-opt-model for more info. 0103 0104 %% check for MOSEK 0105 % if ~have_feature('mosek') 0106 % error('qps_mosek: requires MOSEK'); 0107 % end 0108 0109 %%----- input argument handling ----- 0110 %% gather inputs 0111 if nargin == 1 && isstruct(H) %% problem struct 0112 p = H; 0113 else %% individual args 0114 p = struct('H', H, 'c', c, 'A', A, 'l', l, 'u', u); 0115 if nargin > 5 0116 p.xmin = xmin; 0117 if nargin > 6 0118 p.xmax = xmax; 0119 if nargin > 7 0120 p.x0 = x0; 0121 if nargin > 8 0122 p.opt = opt; 0123 end 0124 end 0125 end 0126 end 0127 end 0128 0129 %% define nx, set default values for H and c 0130 if ~isfield(p, 'H') || isempty(p.H) || ~any(any(p.H)) 0131 if (~isfield(p, 'A') || isempty(p.A)) && ... 0132 (~isfield(p, 'xmin') || isempty(p.xmin)) && ... 0133 (~isfield(p, 'xmax') || isempty(p.xmax)) 0134 error('qps_mosek: LP problem must include constraints or variable bounds'); 0135 else 0136 if isfield(p, 'A') && ~isempty(p.A) 0137 nx = size(p.A, 2); 0138 elseif isfield(p, 'xmin') && ~isempty(p.xmin) 0139 nx = length(p.xmin); 0140 else % if isfield(p, 'xmax') && ~isempty(p.xmax) 0141 nx = length(p.xmax); 0142 end 0143 end 0144 p.H = sparse(nx, nx); 0145 qp = 0; 0146 else 0147 nx = size(p.H, 1); 0148 qp = 1; 0149 end 0150 if ~isfield(p, 'c') || isempty(p.c) 0151 p.c = zeros(nx, 1); 0152 end 0153 if ~isfield(p, 'x0') || isempty(p.x0) 0154 p.x0 = zeros(nx, 1); 0155 end 0156 0157 %% default options 0158 if ~isfield(p, 'opt') 0159 p.opt = []; 0160 end 0161 if ~isempty(p.opt) && isfield(p.opt, 'verbose') && ~isempty(p.opt.verbose) 0162 verbose = p.opt.verbose; 0163 else 0164 verbose = 0; 0165 end 0166 if ~isempty(p.opt) && isfield(p.opt, 'mosek_opt') && ~isempty(p.opt.mosek_opt) 0167 mosek_opt = mosek_options(p.opt.mosek_opt); 0168 else 0169 mosek_opt = mosek_options; 0170 end 0171 0172 %% set up problem struct for MOSEK 0173 prob.c = p.c; 0174 if qp 0175 [prob.qosubi, prob.qosubj, prob.qoval] = find(tril(sparse(p.H))); 0176 end 0177 if isfield(p, 'A') && ~isempty(p.A) 0178 prob.a = sparse(p.A); 0179 nA = size(p.A, 1); 0180 else 0181 nA = 0; 0182 end 0183 if isfield(p, 'l') && ~isempty(p.A) 0184 prob.blc = p.l; 0185 end 0186 if isfield(p, 'u') && ~isempty(p.A) 0187 prob.buc = p.u; 0188 end 0189 if isfield(p, 'xmin') && ~isempty(p.xmin) 0190 prob.blx = p.xmin; 0191 end 0192 if isfield(p, 'xmax') && ~isempty(p.xmax) 0193 prob.bux = p.xmax; 0194 end 0195 0196 %% A is not allowed to be empty 0197 if ~isfield(prob, 'a') || isempty(prob.a) 0198 unconstrained = 1; 0199 prob.a = sparse(1, 1, 1, 1, nx); 0200 prob.blc = -Inf; 0201 prob.buc = Inf; 0202 else 0203 unconstrained = 0; 0204 end 0205 0206 %%----- run optimization ----- 0207 if verbose 0208 s = have_feature('mosek', 'all'); 0209 if s.vnum < 7 0210 alg_names = { %% version 6.x 0211 'default', %% 0 : MSK_OPTIMIZER_FREE 0212 'interior point', %% 1 : MSK_OPTIMIZER_INTPNT 0213 '<conic>', %% 2 : MSK_OPTIMIZER_CONIC 0214 '<qcone>', %% 3 : MSK_OPTIMIZER_QCONE 0215 'primal simplex', %% 4 : MSK_OPTIMIZER_PRIMAL_SIMPLEX 0216 'dual simplex', %% 5 : MSK_OPTIMIZER_DUAL_SIMPLEX 0217 'primal dual simplex', %% 6 : MSK_OPTIMIZER_PRIMAL_DUAL_SIMPLEX 0218 'automatic simplex', %% 7 : MSK_OPTIMIZER_FREE_SIMPLEX 0219 '<mixed int>', %% 8 : MSK_OPTIMIZER_MIXED_INT 0220 '<nonconvex>', %% 9 : MSK_OPTIMIZER_NONCONVEX 0221 'concurrent' %% 10 : MSK_OPTIMIZER_CONCURRENT 0222 }; 0223 elseif s.vnum < 8 0224 alg_names = { %% version 7.x 0225 'default', %% 0 : MSK_OPTIMIZER_FREE 0226 'interior point', %% 1 : MSK_OPTIMIZER_INTPNT 0227 '<conic>', %% 2 : MSK_OPTIMIZER_CONIC 0228 'primal simplex', %% 3 : MSK_OPTIMIZER_PRIMAL_SIMPLEX 0229 'dual simplex', %% 4 : MSK_OPTIMIZER_DUAL_SIMPLEX 0230 'primal dual simplex', %% 5 : MSK_OPTIMIZER_PRIMAL_DUAL_SIMPLEX 0231 'automatic simplex', %% 6 : MSK_OPTIMIZER_FREE_SIMPLEX 0232 'network simplex', %% 7 : MSK_OPTIMIZER_NETWORK_PRIMAL_SIMPLEX 0233 '<mixed int conic>', %% 8 : MSK_OPTIMIZER_MIXED_INT_CONIC 0234 '<mixed int>', %% 9 : MSK_OPTIMIZER_MIXED_INT 0235 'concurrent', %% 10 : MSK_OPTIMIZER_CONCURRENT 0236 '<nonconvex>' %% 11 : MSK_OPTIMIZER_NONCONVEX 0237 }; 0238 else 0239 alg_names = { %% version 8.x 0240 '<conic>', %% 0 : MSK_OPTIMIZER_CONIC 0241 'dual simplex', %% 1 : MSK_OPTIMIZER_DUAL_SIMPLEX 0242 'default', %% 2 : MSK_OPTIMIZER_FREE 0243 'automatic simplex', %% 3 : MSK_OPTIMIZER_FREE_SIMPLEX 0244 'interior point', %% 4 : MSK_OPTIMIZER_INTPNT 0245 '<mixed int>', %% 5 : MSK_OPTIMIZER_MIXED_INT 0246 'primal simplex' %% 6 : MSK_OPTIMIZER_PRIMAL_SIMPLEX 0247 }; 0248 end 0249 if qp 0250 lpqp = 'QP'; 0251 else 0252 lpqp = 'LP'; 0253 end 0254 vn = have_feature('mosek', 'vstr'); 0255 if isempty(vn) 0256 vn = '<unknown>'; 0257 end 0258 fprintf('MOSEK Version %s -- %s %s solver\n', ... 0259 vn, alg_names{mosek_opt.MSK_IPAR_OPTIMIZER+1}, lpqp); 0260 end 0261 cmd = sprintf('minimize echo(%d)', verbose); 0262 [r, res] = mosekopt(cmd, prob, mosek_opt); 0263 0264 %%----- repackage results ----- 0265 if isfield(res, 'sol') 0266 if isfield(res.sol, 'bas') 0267 sol = res.sol.bas; 0268 else 0269 sol = res.sol.itr; 0270 end 0271 x = sol.xx; 0272 else 0273 sol = []; 0274 x = NaN(nx, 1); 0275 end 0276 0277 %%----- process return codes ----- 0278 if isfield(res, 'symbcon') 0279 sc = res.symbcon; 0280 else 0281 sc = mosek_symbcon; 0282 end 0283 eflag = -r; 0284 msg = ''; 0285 switch (r) 0286 case sc.MSK_RES_OK 0287 if ~isempty(sol) 0288 % if sol.solsta == sc.MSK_SOL_STA_OPTIMAL 0289 if strcmp(sol.solsta, 'OPTIMAL') 0290 msg = 'The solution is optimal.'; 0291 eflag = 1; 0292 else 0293 eflag = -1; 0294 % if sol.prosta == sc.MSK_PRO_STA_PRIM_INFEAS 0295 if strcmp(sol.prosta, 'PRIMAL_INFEASIBLE') 0296 msg = 'The problem is primal infeasible.'; 0297 % elseif sol.prosta == sc.MSK_PRO_STA_DUAL_INFEAS 0298 elseif strcmp(sol.prosta, 'DUAL_INFEASIBLE') 0299 msg = 'The problem is dual infeasible.'; 0300 else 0301 msg = sol.solsta; 0302 end 0303 end 0304 end 0305 case sc.MSK_RES_TRM_STALL 0306 if strcmp(sol.solsta, 'OPTIMAL') 0307 msg = 'Stalled at or near optimal solution.'; 0308 eflag = 1; 0309 else 0310 msg = 'Stalled.'; 0311 end 0312 case sc.MSK_RES_TRM_MAX_ITERATIONS 0313 eflag = 0; 0314 msg = 'The optimizer terminated at the maximum number of iterations.'; 0315 otherwise 0316 if isfield(res, 'rmsg') && isfield(res, 'rcodestr') 0317 msg = sprintf('%s : %s', res.rcodestr, res.rmsg); 0318 else 0319 msg = sprintf('MOSEK return code = %d', r); 0320 end 0321 end 0322 0323 if (verbose || r == sc.MSK_RES_ERR_LICENSE || ... 0324 r == sc.MSK_RES_ERR_LICENSE_EXPIRED || ... 0325 r == sc.MSK_RES_ERR_LICENSE_VERSION || ... 0326 r == sc.MSK_RES_ERR_LICENSE_NO_SERVER_SUPPORT || ... 0327 r == sc.MSK_RES_ERR_LICENSE_FEATURE || ... 0328 r == sc.MSK_RES_ERR_LICENSE_INVALID_HOSTID || ... 0329 r == sc.MSK_RES_ERR_LICENSE_SERVER_VERSION || ... 0330 r == sc.MSK_RES_ERR_MISSING_LICENSE_FILE) ... 0331 && ~isempty(msg) %% always alert user of license problems 0332 fprintf('%s\n', msg); 0333 end 0334 0335 %%----- repackage results ----- 0336 if nargout > 1 0337 if eflag == 1 0338 f = p.c' * x; 0339 if ~isempty(p.H) 0340 f = 0.5 * x' * p.H * x + f; 0341 end 0342 else 0343 f = []; 0344 end 0345 if nargout > 3 0346 output.r = r; 0347 output.res = res; 0348 if nargout > 4 0349 if ~isempty(sol) 0350 if isfield(sol, 'slx') 0351 lambda.lower = sol.slx; 0352 else 0353 lambda.lower = []; 0354 end 0355 if isfield(sol, 'sux') 0356 lambda.upper = sol.sux; 0357 else 0358 lambda.upper = []; 0359 end 0360 if isfield(sol, 'slc') 0361 lambda.mu_l = sol.slc; 0362 else 0363 lambda.mu_l = []; 0364 end 0365 if isfield(sol, 'suc') 0366 lambda.mu_u = sol.suc; 0367 else 0368 lambda.mu_u = []; 0369 end 0370 else 0371 if isfield(p, 'xmin') && ~isempty(p.xmin) 0372 lambda.lower = NaN(nx, 1); 0373 else 0374 lambda.lower = []; 0375 end 0376 if isfield(p, 'xmax') && ~isempty(p.xmax) 0377 lambda.upper = NaN(nx, 1); 0378 else 0379 lambda.upper = []; 0380 end 0381 lambda.mu_l = NaN(nA, 1); 0382 lambda.mu_u = NaN(nA, 1); 0383 end 0384 if unconstrained 0385 lambda.mu_l = []; 0386 lambda.mu_u = []; 0387 end 0388 end 0389 end 0390 end