MIQPS_MOSEK Mixed Integer Quadratic Program Solver based on MOSEK.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
MIQPS_MOSEK(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, 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
VTYPE : character string of length NX (number of elements in X),
or 1 (value applies to all variables in x),
allowed values are 'C' (continuous), 'B' (binary),
'I' (integer).
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
skip_prices (0) - flag that specifies whether or not to
skip the price computation stage, in which the problem
is re-solved for only the continuous variables, with all
others being constrained to their solved values
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, vtype, 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] = ...
miqps_mosek(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
x = miqps_mosek(H, c, A, l, u)
x = miqps_mosek(H, c, A, l, u, xmin, xmax)
x = miqps_mosek(H, c, A, l, u, xmin, xmax, x0)
x = miqps_mosek(H, c, A, l, u, xmin, xmax, x0, vtype)
x = miqps_mosek(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
x = miqps_mosek(problem), where problem is a struct with fields:
H, c, A, l, u, xmin, xmax, x0, vtype, opt
all fields except 'c', 'A' and 'l' or 'u' are optional
x = miqps_mosek(...)
[x, f] = miqps_mosek(...)
[x, f, exitflag] = miqps_mosek(...)
[x, f, exitflag, output] = miqps_mosek(...)
[x, f, exitflag, output, lambda] = miqps_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] = miqps_mosek(H, c, A, l, u, xmin, [], x0, vtype, opt);
See also MOSEKOPT.0001 function [x, f, eflag, output, lambda] = miqps_mosek(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0002 %MIQPS_MOSEK Mixed Integer Quadratic Program Solver based on MOSEK. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % MIQPS_MOSEK(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, 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 % VTYPE : character string of length NX (number of elements in X), 0026 % or 1 (value applies to all variables in x), 0027 % allowed values are 'C' (continuous), 'B' (binary), 0028 % 'I' (integer). 0029 % OPT : optional options structure with the following fields, 0030 % all of which are also optional (default values shown in 0031 % parentheses) 0032 % verbose (0) - controls level of progress output displayed 0033 % 0 = no progress output 0034 % 1 = some progress output 0035 % 2 = verbose progress output 0036 % skip_prices (0) - flag that specifies whether or not to 0037 % skip the price computation stage, in which the problem 0038 % is re-solved for only the continuous variables, with all 0039 % others being constrained to their solved values 0040 % mosek_opt - options struct for MOSEK, value in verbose 0041 % overrides these options 0042 % PROBLEM : The inputs can alternatively be supplied in a single 0043 % PROBLEM struct with fields corresponding to the input arguments 0044 % described above: H, c, A, l, u, xmin, xmax, x0, vtype, opt 0045 % 0046 % Outputs: 0047 % X : solution vector 0048 % F : final objective function value 0049 % EXITFLAG : exit flag 0050 % 1 = success 0051 % 0 = terminated at maximum number of iterations 0052 % -1 = primal or dual infeasible 0053 % < 0 = the negative of the MOSEK return code 0054 % OUTPUT : output struct with the following fields: 0055 % r - MOSEK return code 0056 % res - MOSEK result struct 0057 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0058 % multipliers on the constraints, with fields: 0059 % mu_l - lower (left-hand) limit on linear constraints 0060 % mu_u - upper (right-hand) limit on linear constraints 0061 % lower - lower bound on optimization variables 0062 % upper - upper bound on optimization variables 0063 % 0064 % Note the calling syntax is almost identical to that of QUADPROG 0065 % from MathWorks' Optimization Toolbox. The main difference is that 0066 % the linear constraints are specified with A, L, U instead of 0067 % A, B, Aeq, Beq. 0068 % 0069 % Calling syntax options: 0070 % [x, f, exitflag, output, lambda] = ... 0071 % miqps_mosek(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0072 % 0073 % x = miqps_mosek(H, c, A, l, u) 0074 % x = miqps_mosek(H, c, A, l, u, xmin, xmax) 0075 % x = miqps_mosek(H, c, A, l, u, xmin, xmax, x0) 0076 % x = miqps_mosek(H, c, A, l, u, xmin, xmax, x0, vtype) 0077 % x = miqps_mosek(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0078 % x = miqps_mosek(problem), where problem is a struct with fields: 0079 % H, c, A, l, u, xmin, xmax, x0, vtype, opt 0080 % all fields except 'c', 'A' and 'l' or 'u' are optional 0081 % x = miqps_mosek(...) 0082 % [x, f] = miqps_mosek(...) 0083 % [x, f, exitflag] = miqps_mosek(...) 0084 % [x, f, exitflag, output] = miqps_mosek(...) 0085 % [x, f, exitflag, output, lambda] = miqps_mosek(...) 0086 % 0087 % Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm) 0088 % H = [ 1003.1 4.3 6.3 5.9; 0089 % 4.3 2.2 2.1 3.9; 0090 % 6.3 2.1 3.5 4.8; 0091 % 5.9 3.9 4.8 10 ]; 0092 % c = zeros(4,1); 0093 % A = [ 1 1 1 1; 0094 % 0.17 0.11 0.10 0.18 ]; 0095 % l = [1; 0.10]; 0096 % u = [1; Inf]; 0097 % xmin = zeros(4,1); 0098 % x0 = [1; 0; 0; 1]; 0099 % opt = struct('verbose', 2); 0100 % [x, f, s, out, lambda] = miqps_mosek(H, c, A, l, u, xmin, [], x0, vtype, opt); 0101 % 0102 % See also MOSEKOPT. 0103 0104 % MATPOWER 0105 % Copyright (c) 2010-2015 by Power System Engineering Research Center (PSERC) 0106 % by Ray Zimmerman, PSERC Cornell 0107 % 0108 % $Id: miqps_mosek.m 2644 2015-03-11 19:34:22Z ray $ 0109 % 0110 % This file is part of MATPOWER. 0111 % Covered by the 3-clause BSD License (see LICENSE file for details). 0112 % See http://www.pserc.cornell.edu/matpower/ for more info. 0113 0114 %% check for Optimization Toolbox 0115 % if ~have_fcn('mosek') 0116 % error('miqps_mosek: requires MOSEK'); 0117 % end 0118 0119 %%----- input argument handling ----- 0120 %% gather inputs 0121 if nargin == 1 && isstruct(H) %% problem struct 0122 p = H; 0123 else %% individual args 0124 p = struct('H', H, 'c', c, 'A', A, 'l', l, 'u', u); 0125 if nargin > 5 0126 p.xmin = xmin; 0127 if nargin > 6 0128 p.xmax = xmax; 0129 if nargin > 7 0130 p.x0 = x0; 0131 if nargin > 8 0132 p.vtype = vtype; 0133 if nargin > 9 0134 p.opt = opt; 0135 end 0136 end 0137 end 0138 end 0139 end 0140 end 0141 0142 %% define nx, set default values for H and c 0143 if ~isfield(p, 'H') || isempty(p.H) || ~any(any(p.H)) 0144 if (~isfield(p, 'A') || isempty(p.A)) && ... 0145 (~isfield(p, 'xmin') || isempty(p.xmin)) && ... 0146 (~isfield(p, 'xmax') || isempty(p.xmax)) 0147 error('miqps_mosek: LP problem must include constraints or variable bounds'); 0148 else 0149 if isfield(p, 'A') && ~isempty(p.A) 0150 nx = size(p.A, 2); 0151 elseif isfield(p, 'xmin') && ~isempty(p.xmin) 0152 nx = length(p.xmin); 0153 else % if isfield(p, 'xmax') && ~isempty(p.xmax) 0154 nx = length(p.xmax); 0155 end 0156 end 0157 p.H = sparse(nx, nx); 0158 qp = 0; 0159 else 0160 nx = size(p.H, 1); 0161 qp = 1; 0162 end 0163 if ~isfield(p, 'c') || isempty(p.c) 0164 p.c = zeros(nx, 1); 0165 end 0166 if ~isfield(p, 'x0') || isempty(p.x0) 0167 p.x0 = zeros(nx, 1); 0168 end 0169 if ~isfield(p, 'vtype') || isempty(p.vtype) 0170 p.vtype = ''; 0171 end 0172 0173 %% default options 0174 if ~isfield(p, 'opt') 0175 p.opt = []; 0176 end 0177 if ~isempty(p.opt) && isfield(p.opt, 'verbose') && ~isempty(p.opt.verbose) 0178 verbose = p.opt.verbose; 0179 else 0180 verbose = 0; 0181 end 0182 if ~isempty(p.opt) && isfield(p.opt, 'mosek_opt') && ~isempty(p.opt.mosek_opt) 0183 mosek_opt = mosek_options(p.opt.mosek_opt); 0184 else 0185 mosek_opt = mosek_options; 0186 end 0187 if qp 0188 mosek_opt.MSK_IPAR_OPTIMIZER = 0; %% default solver only for QP 0189 end 0190 0191 %% set up problem struct for MOSEK 0192 prob.c = p.c; 0193 if qp 0194 [prob.qosubi, prob.qosubj, prob.qoval] = find(tril(sparse(p.H))); 0195 end 0196 if isfield(p, 'A') && ~isempty(p.A) 0197 prob.a = sparse(p.A); 0198 nA = size(p.A, 1); 0199 else 0200 nA = 0; 0201 end 0202 if isfield(p, 'l') && ~isempty(p.A) 0203 prob.blc = p.l; 0204 end 0205 if isfield(p, 'u') && ~isempty(p.A) 0206 prob.buc = p.u; 0207 end 0208 if ~isempty(p.vtype) 0209 if length(p.vtype) == 1 0210 if p.vtype == 'I' 0211 prob.ints.sub = (1:nx); 0212 elseif p.vtype == 'B' 0213 prob.ints.sub = (1:nx); 0214 p.xmin = zeros(nx, 1); 0215 p.xmax = ones(nx, 1); 0216 end 0217 else 0218 k = find(p.vtype == 'B' | p.vtype == 'I'); 0219 prob.ints.sub = k; 0220 k = find(p.vtype == 'B'); 0221 if ~isempty(k) 0222 if isempty(p.xmin) 0223 p.xmin = -Inf(nx, 1); 0224 end 0225 if isempty(p.xmax) 0226 p.xmax = Inf(nx, 1); 0227 end 0228 p.xmin(k) = 0; 0229 p.xmax(k) = 1; 0230 end 0231 end 0232 end 0233 if isfield(p, 'xmin') && ~isempty(p.xmin) 0234 prob.blx = p.xmin; 0235 end 0236 if isfield(p, 'xmax') && ~isempty(p.xmax) 0237 prob.bux = p.xmax; 0238 end 0239 0240 %% A is not allowed to be empty 0241 if ~isfield(prob, 'a') || isempty(prob.a) 0242 unconstrained = 1; 0243 prob.a = sparse(1, 1, 1, 1, nx); 0244 prob.blc = -Inf; 0245 prob.buc = Inf; 0246 else 0247 unconstrained = 0; 0248 end 0249 if isfield(prob, 'ints') && isfield(prob.ints, 'sub') && ~isempty(prob.ints.sub) 0250 mi = 1; 0251 else 0252 mi = 0; 0253 end 0254 0255 %%----- run optimization ----- 0256 if verbose 0257 s = have_fcn('mosek', 'all'); 0258 if s.vnum < 7 0259 methods = { 0260 'default', %% 0 : MSK_OPTIMIZER_FREE 0261 'interior point', %% 1 : MSK_OPTIMIZER_INTPNT 0262 '<conic>', %% 2 : MSK_OPTIMIZER_CONIC 0263 '<qcone>', %% 3 : MSK_OPTIMIZER_QCONE 0264 'primal simplex', %% 4 : MSK_OPTIMIZER_PRIMAL_SIMPLEX 0265 'dual simplex', %% 5 : MSK_OPTIMIZER_DUAL_SIMPLEX 0266 'primal dual simplex', %% 6 : MSK_OPTIMIZER_PRIMAL_DUAL_SIMPLEX 0267 'automatic simplex', %% 7 : MSK_OPTIMIZER_FREE_SIMPLEX 0268 '<mixed int>', %% 8 : MSK_OPTIMIZER_MIXED_INT 0269 '<nonconvex>', %% 9 : MSK_OPTIMIZER_NONCONVEX 0270 'concurrent' %% 10 : MSK_OPTIMIZER_CONCURRENT 0271 }; 0272 else 0273 methods = { 0274 'default', %% 0 : MSK_OPTIMIZER_FREE 0275 'interior point', %% 1 : MSK_OPTIMIZER_INTPNT 0276 '<conic>', %% 2 : MSK_OPTIMIZER_CONIC 0277 'primal simplex', %% 3 : MSK_OPTIMIZER_PRIMAL_SIMPLEX 0278 'dual simplex', %% 4 : MSK_OPTIMIZER_DUAL_SIMPLEX 0279 'primal dual simplex', %% 5 : MSK_OPTIMIZER_PRIMAL_DUAL_SIMPLEX 0280 'automatic simplex', %% 6 : MSK_OPTIMIZER_FREE_SIMPLEX 0281 'network simplex', %% 7 : MSK_OPTIMIZER_NETWORK_PRIMAL_SIMPLEX 0282 '<mixed int conic>', %% 8 : MSK_OPTIMIZER_MIXED_INT_CONIC 0283 '<mixed int>', %% 9 : MSK_OPTIMIZER_MIXED_INT 0284 'concurrent', %% 10 : MSK_OPTIMIZER_CONCURRENT 0285 '<nonconvex>' %% 11 : MSK_OPTIMIZER_NONCONVEX 0286 }; 0287 end 0288 if qp 0289 lpqp = 'QP'; 0290 else 0291 lpqp = 'LP'; 0292 end 0293 if mi 0294 lpqp = ['MI' lpqp]; 0295 end 0296 vn = have_fcn('mosek', 'vstr'); 0297 if isempty(vn) 0298 vn = '<unknown>'; 0299 end 0300 fprintf('MOSEK Version %s -- %s %s solver\n', ... 0301 vn, methods{mosek_opt.MSK_IPAR_OPTIMIZER+1}, lpqp); 0302 end 0303 cmd = sprintf('minimize echo(%d)', verbose); 0304 [r, res] = mosekopt(cmd, prob, mosek_opt); 0305 0306 %%----- repackage results ----- 0307 if isfield(res, 'sol') 0308 if isfield(res.sol, 'int') 0309 sol = res.sol.int; 0310 elseif isfield(res.sol, 'bas') 0311 sol = res.sol.bas; 0312 else 0313 sol = res.sol.itr; 0314 end 0315 x = sol.xx; 0316 else 0317 sol = []; 0318 x = NaN(nx, 1); 0319 end 0320 0321 %%----- process return codes ----- 0322 if isfield(res, 'symbcon') 0323 sc = res.symbcon; 0324 else 0325 sc = mosek_symbcon; 0326 end 0327 eflag = -r; 0328 msg = ''; 0329 switch (r) 0330 case sc.MSK_RES_OK 0331 if ~isempty(sol) 0332 % if sol.solsta == sc.MSK_SOL_STA_OPTIMAL 0333 if strcmp(sol.solsta, 'OPTIMAL') || strcmp(sol.solsta, 'INTEGER_OPTIMAL') 0334 msg = 'The solution is optimal.'; 0335 eflag = 1; 0336 else 0337 eflag = -1; 0338 % if sol.prosta == sc.MSK_PRO_STA_PRIM_INFEAS 0339 if strcmp(sol.prosta, 'PRIMAL_INFEASIBLE') 0340 msg = 'The problem is primal infeasible.'; 0341 % elseif sol.prosta == sc.MSK_PRO_STA_DUAL_INFEAS 0342 elseif strcmp(sol.prosta, 'DUAL_INFEASIBLE') 0343 msg = 'The problem is dual infeasible.'; 0344 else 0345 msg = sol.solsta; 0346 end 0347 end 0348 end 0349 case sc.MSK_RES_TRM_MAX_ITERATIONS 0350 eflag = 0; 0351 msg = 'The optimizer terminated at the maximum number of iterations.'; 0352 otherwise 0353 if isfield(res, 'rmsg') && isfield(res, 'rcodestr') 0354 msg = sprintf('%s : %s', res.rcodestr, res.rmsg); 0355 else 0356 msg = sprintf('MOSEK return code = %d', r); 0357 end 0358 end 0359 0360 if (verbose || r == sc.MSK_RES_ERR_LICENSE || ... 0361 r == sc.MSK_RES_ERR_LICENSE_EXPIRED || ... 0362 r == sc.MSK_RES_ERR_LICENSE_VERSION || ... 0363 r == sc.MSK_RES_ERR_LICENSE_NO_SERVER_SUPPORT || ... 0364 r == sc.MSK_RES_ERR_LICENSE_FEATURE || ... 0365 r == sc.MSK_RES_ERR_LICENSE_INVALID_HOSTID || ... 0366 r == sc.MSK_RES_ERR_LICENSE_SERVER_VERSION || ... 0367 r == sc.MSK_RES_ERR_MISSING_LICENSE_FILE) ... 0368 && ~isempty(msg) %% always alert user of license problems 0369 fprintf('%s\n', msg); 0370 end 0371 0372 %%----- repackage results ----- 0373 if nargout > 1 0374 if r == 0 0375 f = p.c' * x; 0376 if ~isempty(p.H) 0377 f = 0.5 * x' * p.H * x + f; 0378 end 0379 else 0380 f = []; 0381 end 0382 if nargout > 3 0383 output.r = r; 0384 output.res = res; 0385 if nargout > 4 0386 if ~isempty(sol) 0387 if isfield(sol, 'slx') 0388 lambda.lower = sol.slx; 0389 else 0390 lambda.lower = []; 0391 end 0392 if isfield(sol, 'sux') 0393 lambda.upper = sol.sux; 0394 else 0395 lambda.upper = []; 0396 end 0397 if isfield(sol, 'slc') 0398 lambda.mu_l = sol.slc; 0399 else 0400 lambda.mu_l = []; 0401 end 0402 if isfield(sol, 'suc') 0403 lambda.mu_u = sol.suc; 0404 else 0405 lambda.mu_u = []; 0406 end 0407 else 0408 if isfield(p, 'xmin') && ~isempty(p.xmin) 0409 lambda.lower = NaN(nx, 1); 0410 else 0411 lambda.lower = []; 0412 end 0413 if isfield(p, 'xmax') && ~isempty(p.xmax) 0414 lambda.upper = NaN(nx, 1); 0415 else 0416 lambda.upper = []; 0417 end 0418 lambda.mu_l = NaN(nA, 1); 0419 lambda.mu_u = NaN(nA, 1); 0420 end 0421 if unconstrained 0422 lambda.mu_l = []; 0423 lambda.mu_u = []; 0424 end 0425 end 0426 end 0427 end 0428 0429 if mi && eflag == 1 && (~isfield(p.opt, 'skip_prices') || ~p.opt.skip_prices) 0430 if verbose 0431 fprintf('--- Integer stage complete, starting price computation stage ---\n'); 0432 end 0433 tol = 1e-7; 0434 qp = p; 0435 x(prob.ints.sub) = round(x(prob.ints.sub)); 0436 qp.xmin(prob.ints.sub) = x(prob.ints.sub); 0437 qp.xmax(prob.ints.sub) = x(prob.ints.sub); 0438 qp.x0 = x; 0439 qp.opt.mosek_opt.MSK_IPAR_OPTIMIZER = sc.MSK_OPTIMIZER_PRIMAL_SIMPLEX; 0440 [x_, f_, eflag_, output_, lambda] = qps_mosek(qp); 0441 if eflag ~= eflag_ 0442 error('miqps_mosek: EXITFLAG from price computation stage = %d', eflag_); 0443 end 0444 if abs(f - f_)/max(abs(f), 1) > tol 0445 warning('miqps_mosek: relative mismatch in objective function value from price computation stage = %g', abs(f - f_)/max(abs(f), 1)); 0446 end 0447 xn = x; 0448 xn(abs(xn)<1) = 1; 0449 [mx, k] = max(abs(x - x_) ./ xn); 0450 if mx > tol 0451 warning('miqps_mosek: max relative mismatch in x from price computation stage = %g (%g)', mx, x(k)); 0452 end 0453 output.price_stage = output_; 0454 end