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