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
max_it (0) - maximum number of iterations allowed
0 = use algorithm default
mosek_opt - options struct for MOSEK, values in
verbose and max_it override 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, lam] = 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 % max_it (0) - maximum number of iterations allowed 0033 % 0 = use algorithm default 0034 % mosek_opt - options struct for MOSEK, values in 0035 % verbose and max_it override these options 0036 % PROBLEM : The inputs can alternatively be supplied in a single 0037 % PROBLEM struct with fields corresponding to the input arguments 0038 % described above: H, c, A, l, u, xmin, xmax, x0, opt 0039 % 0040 % Outputs: 0041 % X : solution vector 0042 % F : final objective function value 0043 % EXITFLAG : exit flag 0044 % 1 = success 0045 % 0 = terminated at maximum number of iterations 0046 % -1 = primal or dual infeasible 0047 % < 0 = the negative of the MOSEK return code 0048 % OUTPUT : output struct with the following fields: 0049 % r - MOSEK return code 0050 % res - MOSEK result struct 0051 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0052 % multipliers on the constraints, with fields: 0053 % mu_l - lower (left-hand) limit on linear constraints 0054 % mu_u - upper (right-hand) limit on linear constraints 0055 % lower - lower bound on optimization variables 0056 % upper - upper bound on optimization variables 0057 % 0058 % Note the calling syntax is almost identical to that of QUADPROG 0059 % from MathWorks' Optimization Toolbox. The main difference is that 0060 % the linear constraints are specified with A, L, U instead of 0061 % A, B, Aeq, Beq. 0062 % 0063 % Calling syntax options: 0064 % [x, f, exitflag, output, lambda] = ... 0065 % qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt) 0066 % 0067 % x = qps_mosek(H, c, A, l, u) 0068 % x = qps_mosek(H, c, A, l, u, xmin, xmax) 0069 % x = qps_mosek(H, c, A, l, u, xmin, xmax, x0) 0070 % x = qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt) 0071 % x = qps_mosek(problem), where problem is a struct with fields: 0072 % H, c, A, l, u, xmin, xmax, x0, opt 0073 % all fields except 'c', 'A' and 'l' or 'u' are optional 0074 % x = qps_mosek(...) 0075 % [x, f] = qps_mosek(...) 0076 % [x, f, exitflag] = qps_mosek(...) 0077 % [x, f, exitflag, output] = qps_mosek(...) 0078 % [x, f, exitflag, output, lambda] = qps_mosek(...) 0079 % 0080 % Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm) 0081 % H = [ 1003.1 4.3 6.3 5.9; 0082 % 4.3 2.2 2.1 3.9; 0083 % 6.3 2.1 3.5 4.8; 0084 % 5.9 3.9 4.8 10 ]; 0085 % c = zeros(4,1); 0086 % A = [ 1 1 1 1; 0087 % 0.17 0.11 0.10 0.18 ]; 0088 % l = [1; 0.10]; 0089 % u = [1; Inf]; 0090 % xmin = zeros(4,1); 0091 % x0 = [1; 0; 0; 1]; 0092 % opt = struct('verbose', 2); 0093 % [x, f, s, out, lam] = qps_mosek(H, c, A, l, u, xmin, [], x0, opt); 0094 % 0095 % See also MOSEKOPT. 0096 0097 % MATPOWER 0098 % $Id: qps_mosek.m,v 1.2 2010/12/15 18:42:47 cvs Exp $ 0099 % by Ray Zimmerman, PSERC Cornell 0100 % Copyright (c) 2010 by Power System Engineering Research Center (PSERC) 0101 % 0102 % This file is part of MATPOWER. 0103 % See http://www.pserc.cornell.edu/matpower/ for more info. 0104 % 0105 % MATPOWER is free software: you can redistribute it and/or modify 0106 % it under the terms of the GNU General Public License as published 0107 % by the Free Software Foundation, either version 3 of the License, 0108 % or (at your option) any later version. 0109 % 0110 % MATPOWER is distributed in the hope that it will be useful, 0111 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0112 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0113 % GNU General Public License for more details. 0114 % 0115 % You should have received a copy of the GNU General Public License 0116 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0117 % 0118 % Additional permission under GNU GPL version 3 section 7 0119 % 0120 % If you modify MATPOWER, or any covered work, to interface with 0121 % other modules (such as MATLAB code and MEX-files) available in a 0122 % MATLAB(R) or comparable environment containing parts covered 0123 % under other licensing terms, the licensors of MATPOWER grant 0124 % you additional permission to convey the resulting work. 0125 0126 %% check for Optimization Toolbox 0127 % if ~have_fcn('mosek') 0128 % error('qps_mosek: requires MOSEK'); 0129 % end 0130 0131 %%----- input argument handling ----- 0132 %% gather inputs 0133 if nargin == 1 && isstruct(H) %% problem struct 0134 p = H; 0135 else %% individual args 0136 p = struct('H', H, 'c', c, 'A', A, 'l', l, 'u', u); 0137 if nargin > 5 0138 p.xmin = xmin; 0139 if nargin > 6 0140 p.xmax = xmax; 0141 if nargin > 7 0142 p.x0 = x0; 0143 if nargin > 8 0144 p.opt = opt; 0145 end 0146 end 0147 end 0148 end 0149 end 0150 0151 %% define nx, set default values for H and c 0152 if ~isfield(p, 'H') || isempty(p.H) || ~any(any(p.H)) 0153 if (~isfield(p, 'A') || isempty(p.A)) && ... 0154 (~isfield(p, 'xmin') || isempty(p.xmin)) && ... 0155 (~isfield(p, 'xmax') || isempty(p.xmax)) 0156 error('qps_mosek: LP problem must include constraints or variable bounds'); 0157 else 0158 if isfield(p, 'A') && ~isempty(p.A) 0159 nx = size(p.A, 2); 0160 elseif isfield(p, 'xmin') && ~isempty(p.xmin) 0161 nx = length(p.xmin); 0162 else % if isfield(p, 'xmax') && ~isempty(p.xmax) 0163 nx = length(p.xmax); 0164 end 0165 end 0166 p.H = sparse(nx, nx); 0167 qp = 0; 0168 else 0169 nx = size(p.H, 1); 0170 qp = 1; 0171 end 0172 if ~isfield(p, 'c') || isempty(p.c) 0173 p.c = zeros(nx, 1); 0174 end 0175 if ~isfield(p, 'x0') || isempty(p.x0) 0176 p.x0 = zeros(nx, 1); 0177 end 0178 0179 %% default options 0180 if ~isfield(p, 'opt') 0181 p.opt = []; 0182 end 0183 if ~isempty(p.opt) && isfield(p.opt, 'verbose') && ~isempty(p.opt.verbose) 0184 verbose = p.opt.verbose; 0185 else 0186 verbose = 0; 0187 end 0188 if ~isempty(p.opt) && isfield(p.opt, 'max_it') && ~isempty(p.opt.max_it) 0189 max_it = p.opt.max_it; 0190 else 0191 max_it = 0; 0192 end 0193 if ~isempty(p.opt) && isfield(p.opt, 'mosek_opt') && ~isempty(p.opt.mosek_opt) 0194 mosek_opt = mosek_options(p.opt.mosek_opt); 0195 else 0196 mosek_opt = mosek_options; 0197 end 0198 if max_it 0199 mosek_opt.MSK_IPAR_INTPNT_MAX_ITERATIONS = max_it; 0200 end 0201 if qp 0202 mosek_opt.MSK_IPAR_OPTIMIZER = 0; %% default solver only for QP 0203 end 0204 0205 %% set up problem struct for MOSEK 0206 prob.c = p.c; 0207 if qp 0208 [prob.qosubi, prob.qosubj, prob.qoval] = find(tril(sparse(p.H))); 0209 end 0210 if isfield(p, 'A') && ~isempty(p.A) 0211 prob.a = sparse(p.A); 0212 end 0213 if isfield(p, 'l') && ~isempty(p.A) 0214 prob.blc = p.l; 0215 end 0216 if isfield(p, 'u') && ~isempty(p.A) 0217 prob.buc = p.u; 0218 end 0219 if isfield(p, 'xmin') && ~isempty(p.xmin) 0220 prob.blx = p.xmin; 0221 end 0222 if isfield(p, 'xmax') && ~isempty(p.xmax) 0223 prob.bux = p.xmax; 0224 end 0225 0226 %% A is not allowed to be empty 0227 if ~isfield(prob, 'a') || isempty(prob.a) 0228 unconstrained = 1; 0229 prob.a = sparse(1, 1, 1, 1, nx); 0230 prob.blc = -Inf; 0231 prob.buc = Inf; 0232 else 0233 unconstrained = 0; 0234 end 0235 0236 %%----- run optimization ----- 0237 cmd = sprintf('minimize echo(%d)', verbose); 0238 [r, res] = mosekopt(cmd, prob, mosek_opt); 0239 0240 %%----- repackage results ----- 0241 if isfield(res, 'sol') 0242 if isfield(res.sol, 'bas') 0243 sol = res.sol.bas; 0244 else 0245 sol = res.sol.itr; 0246 end 0247 x = sol.xx; 0248 else 0249 sol = []; 0250 x = []; 0251 end 0252 0253 %%----- process return codes ----- 0254 if isfield(res, 'symbcon') 0255 sc = res.symbcon; 0256 else 0257 [r2, res2] = mosekopt('symbcon echo(0)'); 0258 sc = res2.symbcon; 0259 end 0260 eflag = -r; 0261 msg = ''; 0262 switch (r) 0263 case sc.MSK_RES_OK 0264 if ~isempty(sol) 0265 % if sol.solsta == sc.MSK_SOL_STA_OPTIMAL 0266 if strcmp(sol.solsta, 'OPTIMAL') 0267 msg = 'The solution is optimal.'; 0268 eflag = 1; 0269 else 0270 eflag = -1; 0271 % if sol.prosta == sc.MSK_PRO_STA_PRIM_INFEAS 0272 if strcmp(sol.prosta, 'PRIMAL_INFEASIBLE') 0273 msg = 'The problem is primal infeasible.'; 0274 % elseif sol.prosta == sc.MSK_PRO_STA_DUAL_INFEAS 0275 elseif strcmp(sol.prosta, 'DUAL_INFEASIBLE') 0276 msg = 'The problem is dual infeasible.'; 0277 else 0278 msg = sol.solsta; 0279 end 0280 end 0281 end 0282 case sc.MSK_RES_TRM_MAX_ITERATIONS 0283 eflag = 0; 0284 msg = 'The optimizer terminated at the maximum number of iterations.'; 0285 otherwise 0286 if isfield(res, 'rmsg') && isfield(res, 'rcodestr') 0287 msg = sprintf('%s : %s', res.rcodestr, res.rmsg); 0288 else 0289 msg = sprintf('MOSEK return code = %d', r); 0290 end 0291 end 0292 0293 if verbose && ~isempty(msg) 0294 fprintf('%s\n', msg); 0295 end 0296 0297 %%----- repackage results ----- 0298 if nargout > 1 0299 if r == 0 0300 f = p.c' * x; 0301 if ~isempty(p.H) 0302 f = 0.5 * x' * p.H * x + f; 0303 end 0304 else 0305 f = []; 0306 end 0307 if nargout > 3 0308 output.r = r; 0309 output.res = res; 0310 if nargout > 4 0311 if isfield(res, 'sol') 0312 lambda.lower = sol.slx; 0313 lambda.upper = sol.sux; 0314 lambda.mu_l = sol.slc; 0315 lambda.mu_u = sol.suc; 0316 if unconstrained 0317 lambda.mu_l = []; 0318 lambda.mu_u = []; 0319 end 0320 else 0321 lambda = []; 0322 end 0323 end 0324 end 0325 end