Home > matpower5.0 > qps_mosek.m

qps_mosek

PURPOSE ^

QPS_MOSEK Quadratic Program Solver based on MOSEK.

SYNOPSIS ^

function [x, f, eflag, output, lambda] = qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt)

DESCRIPTION ^

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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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 %   $Id: qps_mosek.m 2340 2014-06-27 19:15:10Z ray $
0097 %   by Ray Zimmerman, PSERC Cornell
0098 %   Copyright (c) 2010 by Power System Engineering Research Center (PSERC)
0099 %
0100 %   This file is part of MATPOWER.
0101 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0102 %
0103 %   MATPOWER is free software: you can redistribute it and/or modify
0104 %   it under the terms of the GNU General Public License as published
0105 %   by the Free Software Foundation, either version 3 of the License,
0106 %   or (at your option) any later version.
0107 %
0108 %   MATPOWER is distributed in the hope that it will be useful,
0109 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0110 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0111 %   GNU General Public License for more details.
0112 %
0113 %   You should have received a copy of the GNU General Public License
0114 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0115 %
0116 %   Additional permission under GNU GPL version 3 section 7
0117 %
0118 %   If you modify MATPOWER, or any covered work, to interface with
0119 %   other modules (such as MATLAB code and MEX-files) available in a
0120 %   MATLAB(R) or comparable environment containing parts covered
0121 %   under other licensing terms, the licensors of MATPOWER grant
0122 %   you additional permission to convey the resulting work.
0123 
0124 %% check for Optimization Toolbox
0125 % if ~have_fcn('mosek')
0126 %     error('qps_mosek: requires MOSEK');
0127 % end
0128 
0129 %%----- input argument handling  -----
0130 %% gather inputs
0131 if nargin == 1 && isstruct(H)       %% problem struct
0132     p = H;
0133 else                                %% individual args
0134     p = struct('H', H, 'c', c, 'A', A, 'l', l, 'u', u);
0135     if nargin > 5
0136         p.xmin = xmin;
0137         if nargin > 6
0138             p.xmax = xmax;
0139             if nargin > 7
0140                 p.x0 = x0;
0141                 if nargin > 8
0142                     p.opt = opt;
0143                 end
0144             end
0145         end
0146     end
0147 end
0148 
0149 %% define nx, set default values for H and c
0150 if ~isfield(p, 'H') || isempty(p.H) || ~any(any(p.H))
0151     if (~isfield(p, 'A') || isempty(p.A)) && ...
0152             (~isfield(p, 'xmin') || isempty(p.xmin)) && ...
0153             (~isfield(p, 'xmax') || isempty(p.xmax))
0154         error('qps_mosek: LP problem must include constraints or variable bounds');
0155     else
0156         if isfield(p, 'A') && ~isempty(p.A)
0157             nx = size(p.A, 2);
0158         elseif isfield(p, 'xmin') && ~isempty(p.xmin)
0159             nx = length(p.xmin);
0160         else    % if isfield(p, 'xmax') && ~isempty(p.xmax)
0161             nx = length(p.xmax);
0162         end
0163     end
0164     p.H = sparse(nx, nx);
0165     qp = 0;
0166 else
0167     nx = size(p.H, 1);
0168     qp = 1;
0169 end
0170 if ~isfield(p, 'c') || isempty(p.c)
0171     p.c = zeros(nx, 1);
0172 end
0173 if ~isfield(p, 'x0') || isempty(p.x0)
0174     p.x0 = zeros(nx, 1);
0175 end
0176 
0177 %% default options
0178 if ~isfield(p, 'opt')
0179     p.opt = [];
0180 end
0181 if ~isempty(p.opt) && isfield(p.opt, 'verbose') && ~isempty(p.opt.verbose)
0182     verbose = p.opt.verbose;
0183 else
0184     verbose = 0;
0185 end
0186 if ~isempty(p.opt) && isfield(p.opt, 'mosek_opt') && ~isempty(p.opt.mosek_opt)
0187     mosek_opt = mosek_options(p.opt.mosek_opt);
0188 else
0189     mosek_opt = mosek_options;
0190 end
0191 if qp
0192     mosek_opt.MSK_IPAR_OPTIMIZER = 0;   %% default solver only for QP
0193 end
0194 
0195 %% set up problem struct for MOSEK
0196 prob.c = p.c;
0197 if qp
0198    [prob.qosubi, prob.qosubj, prob.qoval] = find(tril(sparse(p.H)));
0199 end
0200 if isfield(p, 'A') && ~isempty(p.A)
0201     prob.a = sparse(p.A);
0202     nA = size(p.A, 1);
0203 else
0204     nA = 0;
0205 end
0206 if isfield(p, 'l') && ~isempty(p.A)
0207     prob.blc = p.l;
0208 end
0209 if isfield(p, 'u') && ~isempty(p.A)
0210     prob.buc = p.u;
0211 end
0212 if isfield(p, 'xmin') && ~isempty(p.xmin)
0213     prob.blx = p.xmin;
0214 end
0215 if isfield(p, 'xmax') && ~isempty(p.xmax)
0216     prob.bux = p.xmax;
0217 end
0218 
0219 %% A is not allowed to be empty
0220 if ~isfield(prob, 'a') || isempty(prob.a)
0221     unconstrained = 1;
0222     prob.a = sparse(1, 1, 1, 1, nx);
0223     prob.blc = -Inf;
0224     prob.buc =  Inf;
0225 else
0226     unconstrained = 0;
0227 end
0228 
0229 %%-----  run optimization  -----
0230 if verbose
0231     methods = {
0232         'default',
0233         'interior point',
0234         '<default>',
0235         '<default>',
0236         'primal simplex',
0237         'dual simplex',
0238         'primal dual simplex',
0239         'automatic simplex',
0240         '<default>',
0241         '<default>',
0242         'concurrent'
0243     };
0244     if qp
0245         lpqp = 'QP';
0246     else
0247         lpqp = 'LP';
0248     end
0249     % (this code is also in mpver.m)
0250     % MOSEK Version 6.0.0.93 (Build date: 2010-10-26 13:03:27)
0251     % MOSEK Version 6.0.0.106 (Build date: 2011-3-17 10:46:54)
0252 %    pat = 'Version (\.*\d)+.*Build date: (\d\d\d\d-\d\d-\d\d)';
0253     pat = 'Version (\.*\d)+.*Build date: (\d+-\d+-\d+)';
0254     [s,e,tE,m,t] = regexp(evalc('mosekopt'), pat);
0255     if isempty(t)
0256         vn = '<unknown>';
0257     else
0258         vn = t{1}{1};
0259     end
0260     fprintf('MOSEK Version %s -- %s %s solver\n', ...
0261             vn, methods{mosek_opt.MSK_IPAR_OPTIMIZER+1}, lpqp);
0262 end
0263 cmd = sprintf('minimize echo(%d)', verbose);
0264 [r, res] = mosekopt(cmd, prob, mosek_opt);
0265 
0266 %%-----  repackage results  -----
0267 if isfield(res, 'sol')
0268     if isfield(res.sol, 'bas')
0269         sol = res.sol.bas;
0270     else
0271         sol = res.sol.itr;
0272     end
0273     x = sol.xx;
0274 else
0275     sol = [];
0276     x = NaN(nx, 1);
0277 end
0278 
0279 %%-----  process return codes  -----
0280 if isfield(res, 'symbcon')
0281     sc = res.symbcon;
0282 else    
0283     [r2, res2] = mosekopt('symbcon echo(0)');
0284     sc = res2.symbcon;
0285 end
0286 eflag = -r;
0287 msg = '';
0288 switch (r)
0289     case sc.MSK_RES_OK
0290         if ~isempty(sol)
0291 %            if sol.solsta == sc.MSK_SOL_STA_OPTIMAL
0292             if strcmp(sol.solsta, 'OPTIMAL')
0293                 msg = 'The solution is optimal.';
0294                 eflag = 1;
0295             else
0296                 eflag = -1;
0297 %                 if sol.prosta == sc.MSK_PRO_STA_PRIM_INFEAS
0298                 if strcmp(sol.prosta, 'PRIMAL_INFEASIBLE')
0299                     msg = 'The problem is primal infeasible.';
0300 %                 elseif sol.prosta == sc.MSK_PRO_STA_DUAL_INFEAS
0301                 elseif strcmp(sol.prosta, 'DUAL_INFEASIBLE')
0302                     msg = 'The problem is dual infeasible.';
0303                 else
0304                     msg = sol.solsta;
0305                 end
0306             end
0307         end
0308     case sc.MSK_RES_TRM_MAX_ITERATIONS
0309         eflag = 0;
0310         msg = 'The optimizer terminated at the maximum number of iterations.';
0311     otherwise
0312         if isfield(res, 'rmsg') && isfield(res, 'rcodestr')
0313             msg = sprintf('%s : %s', res.rcodestr, res.rmsg);
0314         else
0315             msg = sprintf('MOSEK return code = %d', r);
0316         end
0317 end
0318 
0319 if (verbose || r == 1001) && ~isempty(msg)  %% always alert user if license is expired
0320     fprintf('%s\n', msg);
0321 end
0322 
0323 %%-----  repackage results  -----
0324 if nargout > 1
0325     if r == 0
0326         f = p.c' * x;
0327         if ~isempty(p.H)
0328             f = 0.5 * x' * p.H * x + f;
0329         end
0330     else
0331         f = [];
0332     end
0333     if nargout > 3
0334         output.r = r;
0335         output.res = res;
0336         if nargout > 4
0337             if ~isempty(sol)
0338                 lambda.lower = sol.slx;
0339                 lambda.upper = sol.sux;
0340                 lambda.mu_l  = sol.slc;
0341                 lambda.mu_u  = sol.suc;
0342             else
0343                 if isfield(p, 'xmin') && ~isempty(p.xmin)
0344                     lambda.lower = NaN(nx, 1);
0345                 else
0346                     lambda.lower = [];
0347                 end
0348                 if isfield(p, 'xmax') && ~isempty(p.xmax)
0349                     lambda.upper = NaN(nx, 1);
0350                 else
0351                     lambda.upper = [];
0352                 end
0353                 lambda.mu_l = NaN(nA, 1);
0354                 lambda.mu_u = NaN(nA, 1);
0355             end
0356             if unconstrained
0357                 lambda.mu_l  = [];
0358                 lambda.mu_u  = [];
0359             end
0360         end
0361     end
0362 end

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005