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

Generated on Fri 20-Mar-2015 18:23:34 by m2html © 2005