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-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 http://www.pserc.cornell.edu/matpower/ 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_MAX_ITERATIONS 0305 eflag = 0; 0306 msg = 'The optimizer terminated at the maximum number of iterations.'; 0307 otherwise 0308 if isfield(res, 'rmsg') && isfield(res, 'rcodestr') 0309 msg = sprintf('%s : %s', res.rcodestr, res.rmsg); 0310 else 0311 msg = sprintf('MOSEK return code = %d', r); 0312 end 0313 end 0314 0315 if (verbose || r == sc.MSK_RES_ERR_LICENSE || ... 0316 r == sc.MSK_RES_ERR_LICENSE_EXPIRED || ... 0317 r == sc.MSK_RES_ERR_LICENSE_VERSION || ... 0318 r == sc.MSK_RES_ERR_LICENSE_NO_SERVER_SUPPORT || ... 0319 r == sc.MSK_RES_ERR_LICENSE_FEATURE || ... 0320 r == sc.MSK_RES_ERR_LICENSE_INVALID_HOSTID || ... 0321 r == sc.MSK_RES_ERR_LICENSE_SERVER_VERSION || ... 0322 r == sc.MSK_RES_ERR_MISSING_LICENSE_FILE) ... 0323 && ~isempty(msg) %% always alert user of license problems 0324 fprintf('%s\n', msg); 0325 end 0326 0327 %%----- repackage results ----- 0328 if nargout > 1 0329 if r == 0 0330 f = p.c' * x; 0331 if ~isempty(p.H) 0332 f = 0.5 * x' * p.H * x + f; 0333 end 0334 else 0335 f = []; 0336 end 0337 if nargout > 3 0338 output.r = r; 0339 output.res = res; 0340 if nargout > 4 0341 if ~isempty(sol) 0342 if isfield(sol, 'slx') 0343 lambda.lower = sol.slx; 0344 else 0345 lambda.lower = []; 0346 end 0347 if isfield(sol, 'sux') 0348 lambda.upper = sol.sux; 0349 else 0350 lambda.upper = []; 0351 end 0352 if isfield(sol, 'slc') 0353 lambda.mu_l = sol.slc; 0354 else 0355 lambda.mu_l = []; 0356 end 0357 if isfield(sol, 'suc') 0358 lambda.mu_u = sol.suc; 0359 else 0360 lambda.mu_u = []; 0361 end 0362 else 0363 if isfield(p, 'xmin') && ~isempty(p.xmin) 0364 lambda.lower = NaN(nx, 1); 0365 else 0366 lambda.lower = []; 0367 end 0368 if isfield(p, 'xmax') && ~isempty(p.xmax) 0369 lambda.upper = NaN(nx, 1); 0370 else 0371 lambda.upper = []; 0372 end 0373 lambda.mu_l = NaN(nA, 1); 0374 lambda.mu_u = NaN(nA, 1); 0375 end 0376 if unconstrained 0377 lambda.mu_l = []; 0378 lambda.mu_u = []; 0379 end 0380 end 0381 end 0382 end