Home > matpower7.1 > mp-opt-model > lib > 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)
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_MOSEK(PROBLEM)
   A wrapper function providing a 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 https://v8doc.sas.com/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 QPS_MASTER, 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 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_MOSEK(PROBLEM)
0006 %   A wrapper function providing a standardized interface for using
0007 %   MOSEKOPT to solve the following QP (quadratic programming) problem:
0008 %
0009 %       min 1/2 X'*H*X + C'*X
0010 %        X
0011 %
0012 %   subject to
0013 %
0014 %       L <= A*X <= U       (linear constraints)
0015 %       XMIN <= X <= XMAX   (variable bounds)
0016 %
0017 %   Inputs (all optional except H, C, A and L):
0018 %       H : matrix (possibly sparse) of quadratic cost coefficients
0019 %       C : vector of linear cost coefficients
0020 %       A, L, U : define the optional linear constraints. Default
0021 %           values for the elements of L and U are -Inf and Inf,
0022 %           respectively.
0023 %       XMIN, XMAX : optional lower and upper bounds on the
0024 %           X variables, defaults are -Inf and Inf, respectively.
0025 %       X0 : optional starting value of optimization vector X
0026 %       OPT : optional options structure with the following fields,
0027 %           all of which are also optional (default values shown in
0028 %           parentheses)
0029 %           verbose (0) - controls level of progress output displayed
0030 %               0 = no progress output
0031 %               1 = some progress output
0032 %               2 = verbose progress output
0033 %           mosek_opt - options struct for MOSEK, value in verbose
0034 %                   overrides these options
0035 %       PROBLEM : The inputs can alternatively be supplied in a single
0036 %           PROBLEM struct with fields corresponding to the input arguments
0037 %           described above: H, c, A, l, u, xmin, xmax, x0, opt
0038 %
0039 %   Outputs:
0040 %       X : solution vector
0041 %       F : final objective function value
0042 %       EXITFLAG : exit flag
0043 %             1 = success
0044 %             0 = terminated at maximum number of iterations
0045 %            -1 = primal or dual infeasible
0046 %           < 0 = the negative of the MOSEK return code
0047 %       OUTPUT : output struct with the following fields:
0048 %           r - MOSEK return code
0049 %           res - MOSEK result struct
0050 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0051 %           multipliers on the constraints, with fields:
0052 %           mu_l - lower (left-hand) limit on linear constraints
0053 %           mu_u - upper (right-hand) limit on linear constraints
0054 %           lower - lower bound on optimization variables
0055 %           upper - upper bound on optimization variables
0056 %
0057 %   Note the calling syntax is almost identical to that of QUADPROG
0058 %   from MathWorks' Optimization Toolbox. The main difference is that
0059 %   the linear constraints are specified with A, L, U instead of
0060 %   A, B, Aeq, Beq.
0061 %
0062 %   Calling syntax options:
0063 %       [x, f, exitflag, output, lambda] = ...
0064 %           qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt)
0065 %
0066 %       x = qps_mosek(H, c, A, l, u)
0067 %       x = qps_mosek(H, c, A, l, u, xmin, xmax)
0068 %       x = qps_mosek(H, c, A, l, u, xmin, xmax, x0)
0069 %       x = qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt)
0070 %       x = qps_mosek(problem), where problem is a struct with fields:
0071 %                       H, c, A, l, u, xmin, xmax, x0, opt
0072 %                       all fields except 'c', 'A' and 'l' or 'u' are optional
0073 %       x = qps_mosek(...)
0074 %       [x, f] = qps_mosek(...)
0075 %       [x, f, exitflag] = qps_mosek(...)
0076 %       [x, f, exitflag, output] = qps_mosek(...)
0077 %       [x, f, exitflag, output, lambda] = qps_mosek(...)
0078 %
0079 %   Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm)
0080 %       H = [   1003.1  4.3     6.3     5.9;
0081 %               4.3     2.2     2.1     3.9;
0082 %               6.3     2.1     3.5     4.8;
0083 %               5.9     3.9     4.8     10  ];
0084 %       c = zeros(4,1);
0085 %       A = [   1       1       1       1;
0086 %               0.17    0.11    0.10    0.18    ];
0087 %       l = [1; 0.10];
0088 %       u = [1; Inf];
0089 %       xmin = zeros(4,1);
0090 %       x0 = [1; 0; 0; 1];
0091 %       opt = struct('verbose', 2);
0092 %       [x, f, s, out, lambda] = qps_mosek(H, c, A, l, u, xmin, [], x0, opt);
0093 %
0094 %   See also QPS_MASTER, MOSEKOPT.
0095 
0096 %   MP-Opt-Model
0097 %   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
0098 %   by Ray Zimmerman, PSERC Cornell
0099 %
0100 %   This file is part of MP-Opt-Model.
0101 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0102 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0103 
0104 %% check for MOSEK
0105 % if ~have_feature('mosek')
0106 %     error('qps_mosek: requires MOSEK');
0107 % end
0108 
0109 %%----- input argument handling  -----
0110 %% gather inputs
0111 if nargin == 1 && isstruct(H)       %% problem struct
0112     p = H;
0113 else                                %% individual args
0114     p = struct('H', H, 'c', c, 'A', A, 'l', l, 'u', u);
0115     if nargin > 5
0116         p.xmin = xmin;
0117         if nargin > 6
0118             p.xmax = xmax;
0119             if nargin > 7
0120                 p.x0 = x0;
0121                 if nargin > 8
0122                     p.opt = opt;
0123                 end
0124             end
0125         end
0126     end
0127 end
0128 
0129 %% define nx, set default values for H and c
0130 if ~isfield(p, 'H') || isempty(p.H) || ~any(any(p.H))
0131     if (~isfield(p, 'A') || isempty(p.A)) && ...
0132             (~isfield(p, 'xmin') || isempty(p.xmin)) && ...
0133             (~isfield(p, 'xmax') || isempty(p.xmax))
0134         error('qps_mosek: LP problem must include constraints or variable bounds');
0135     else
0136         if isfield(p, 'A') && ~isempty(p.A)
0137             nx = size(p.A, 2);
0138         elseif isfield(p, 'xmin') && ~isempty(p.xmin)
0139             nx = length(p.xmin);
0140         else    % if isfield(p, 'xmax') && ~isempty(p.xmax)
0141             nx = length(p.xmax);
0142         end
0143     end
0144     p.H = sparse(nx, nx);
0145     qp = 0;
0146 else
0147     nx = size(p.H, 1);
0148     qp = 1;
0149 end
0150 if ~isfield(p, 'c') || isempty(p.c)
0151     p.c = zeros(nx, 1);
0152 end
0153 if ~isfield(p, 'x0') || isempty(p.x0)
0154     p.x0 = zeros(nx, 1);
0155 end
0156 
0157 %% default options
0158 if ~isfield(p, 'opt')
0159     p.opt = [];
0160 end
0161 if ~isempty(p.opt) && isfield(p.opt, 'verbose') && ~isempty(p.opt.verbose)
0162     verbose = p.opt.verbose;
0163 else
0164     verbose = 0;
0165 end
0166 if ~isempty(p.opt) && isfield(p.opt, 'mosek_opt') && ~isempty(p.opt.mosek_opt)
0167     mosek_opt = mosek_options(p.opt.mosek_opt);
0168 else
0169     mosek_opt = mosek_options;
0170 end
0171 
0172 %% set up problem struct for MOSEK
0173 prob.c = p.c;
0174 if qp
0175    [prob.qosubi, prob.qosubj, prob.qoval] = find(tril(sparse(p.H)));
0176 end
0177 if isfield(p, 'A') && ~isempty(p.A)
0178     prob.a = sparse(p.A);
0179     nA = size(p.A, 1);
0180 else
0181     nA = 0;
0182 end
0183 if isfield(p, 'l') && ~isempty(p.A)
0184     prob.blc = p.l;
0185 end
0186 if isfield(p, 'u') && ~isempty(p.A)
0187     prob.buc = p.u;
0188 end
0189 if isfield(p, 'xmin') && ~isempty(p.xmin)
0190     prob.blx = p.xmin;
0191 end
0192 if isfield(p, 'xmax') && ~isempty(p.xmax)
0193     prob.bux = p.xmax;
0194 end
0195 
0196 %% A is not allowed to be empty
0197 if ~isfield(prob, 'a') || isempty(prob.a)
0198     unconstrained = 1;
0199     prob.a = sparse(1, 1, 1, 1, nx);
0200     prob.blc = -Inf;
0201     prob.buc =  Inf;
0202 else
0203     unconstrained = 0;
0204 end
0205 
0206 %%-----  run optimization  -----
0207 if verbose
0208     s = have_feature('mosek', 'all');
0209     if s.vnum < 7
0210         alg_names = {           %% version 6.x
0211             'default',              %%  0 : MSK_OPTIMIZER_FREE
0212             'interior point',       %%  1 : MSK_OPTIMIZER_INTPNT
0213             '<conic>',              %%  2 : MSK_OPTIMIZER_CONIC
0214             '<qcone>',              %%  3 : MSK_OPTIMIZER_QCONE
0215             'primal simplex',       %%  4 : MSK_OPTIMIZER_PRIMAL_SIMPLEX
0216             'dual simplex',         %%  5 : MSK_OPTIMIZER_DUAL_SIMPLEX
0217             'primal dual simplex',  %%  6 : MSK_OPTIMIZER_PRIMAL_DUAL_SIMPLEX
0218             'automatic simplex',    %%  7 : MSK_OPTIMIZER_FREE_SIMPLEX
0219             '<mixed int>',          %%  8 : MSK_OPTIMIZER_MIXED_INT
0220             '<nonconvex>',          %%  9 : MSK_OPTIMIZER_NONCONVEX
0221             'concurrent'            %% 10 : MSK_OPTIMIZER_CONCURRENT
0222         };
0223     elseif s.vnum < 8
0224         alg_names = {           %% version 7.x
0225             'default',              %%  0 : MSK_OPTIMIZER_FREE
0226             'interior point',       %%  1 : MSK_OPTIMIZER_INTPNT
0227             '<conic>',              %%  2 : MSK_OPTIMIZER_CONIC
0228             'primal simplex',       %%  3 : MSK_OPTIMIZER_PRIMAL_SIMPLEX
0229             'dual simplex',         %%  4 : MSK_OPTIMIZER_DUAL_SIMPLEX
0230             'primal dual simplex',  %%  5 : MSK_OPTIMIZER_PRIMAL_DUAL_SIMPLEX
0231             'automatic simplex',    %%  6 : MSK_OPTIMIZER_FREE_SIMPLEX
0232             'network simplex',      %%  7 : MSK_OPTIMIZER_NETWORK_PRIMAL_SIMPLEX
0233             '<mixed int conic>',    %%  8 : MSK_OPTIMIZER_MIXED_INT_CONIC
0234             '<mixed int>',          %%  9 : MSK_OPTIMIZER_MIXED_INT
0235             'concurrent',           %% 10 : MSK_OPTIMIZER_CONCURRENT
0236             '<nonconvex>'           %% 11 : MSK_OPTIMIZER_NONCONVEX
0237         };
0238     else
0239         alg_names = {           %% version 8.x
0240             '<conic>',              %%  0 : MSK_OPTIMIZER_CONIC
0241             'dual simplex',         %%  1 : MSK_OPTIMIZER_DUAL_SIMPLEX
0242             'default',              %%  2 : MSK_OPTIMIZER_FREE
0243             'automatic simplex',    %%  3 : MSK_OPTIMIZER_FREE_SIMPLEX
0244             'interior point',       %%  4 : MSK_OPTIMIZER_INTPNT
0245             '<mixed int>',          %%  5 : MSK_OPTIMIZER_MIXED_INT
0246             'primal simplex'        %%  6 : MSK_OPTIMIZER_PRIMAL_SIMPLEX
0247         };
0248     end
0249     if qp
0250         lpqp = 'QP';
0251     else
0252         lpqp = 'LP';
0253     end
0254     vn = have_feature('mosek', 'vstr');
0255     if isempty(vn)
0256         vn = '<unknown>';
0257     end
0258     fprintf('MOSEK Version %s -- %s %s solver\n', ...
0259             vn, alg_names{mosek_opt.MSK_IPAR_OPTIMIZER+1}, lpqp);
0260 end
0261 cmd = sprintf('minimize echo(%d)', verbose);
0262 [r, res] = mosekopt(cmd, prob, mosek_opt);
0263 
0264 %%-----  repackage results  -----
0265 if isfield(res, 'sol')
0266     if isfield(res.sol, 'bas')
0267         sol = res.sol.bas;
0268     else
0269         sol = res.sol.itr;
0270     end
0271     x = sol.xx;
0272 else
0273     sol = [];
0274     x = NaN(nx, 1);
0275 end
0276 
0277 %%-----  process return codes  -----
0278 if isfield(res, 'symbcon')
0279     sc = res.symbcon;
0280 else    
0281     sc = mosek_symbcon;
0282 end
0283 eflag = -r;
0284 msg = '';
0285 switch (r)
0286     case sc.MSK_RES_OK
0287         if ~isempty(sol)
0288 %            if sol.solsta == sc.MSK_SOL_STA_OPTIMAL
0289             if strcmp(sol.solsta, 'OPTIMAL')
0290                 msg = 'The solution is optimal.';
0291                 eflag = 1;
0292             else
0293                 eflag = -1;
0294 %                 if sol.prosta == sc.MSK_PRO_STA_PRIM_INFEAS
0295                 if strcmp(sol.prosta, 'PRIMAL_INFEASIBLE')
0296                     msg = 'The problem is primal infeasible.';
0297 %                 elseif sol.prosta == sc.MSK_PRO_STA_DUAL_INFEAS
0298                 elseif strcmp(sol.prosta, 'DUAL_INFEASIBLE')
0299                     msg = 'The problem is dual infeasible.';
0300                 else
0301                     msg = sol.solsta;
0302                 end
0303             end
0304         end
0305     case sc.MSK_RES_TRM_STALL
0306         if strcmp(sol.solsta, 'OPTIMAL')
0307             msg = 'Stalled at or near optimal solution.';
0308             eflag = 1;
0309         else
0310             msg = 'Stalled.';
0311         end
0312     case sc.MSK_RES_TRM_MAX_ITERATIONS
0313         eflag = 0;
0314         msg = 'The optimizer terminated at the maximum number of iterations.';
0315     otherwise
0316         if isfield(res, 'rmsg') && isfield(res, 'rcodestr')
0317             msg = sprintf('%s : %s', res.rcodestr, res.rmsg);
0318         else
0319             msg = sprintf('MOSEK return code = %d', r);
0320         end
0321 end
0322 
0323 if (verbose || r == sc.MSK_RES_ERR_LICENSE || ...
0324         r == sc.MSK_RES_ERR_LICENSE_EXPIRED || ...
0325         r == sc.MSK_RES_ERR_LICENSE_VERSION || ...
0326         r == sc.MSK_RES_ERR_LICENSE_NO_SERVER_SUPPORT || ...
0327         r == sc.MSK_RES_ERR_LICENSE_FEATURE || ...
0328         r == sc.MSK_RES_ERR_LICENSE_INVALID_HOSTID || ...
0329         r == sc.MSK_RES_ERR_LICENSE_SERVER_VERSION || ...
0330         r == sc.MSK_RES_ERR_MISSING_LICENSE_FILE) ...
0331         && ~isempty(msg)  %% always alert user of license problems
0332     fprintf('%s\n', msg);
0333 end
0334 
0335 %%-----  repackage results  -----
0336 if nargout > 1
0337     if eflag == 1
0338         f = p.c' * x;
0339         if ~isempty(p.H)
0340             f = 0.5 * x' * p.H * x + f;
0341         end
0342     else
0343         f = [];
0344     end
0345     if nargout > 3
0346         output.r = r;
0347         output.res = res;
0348         if nargout > 4
0349             if ~isempty(sol)
0350                 if isfield(sol, 'slx')
0351                     lambda.lower = sol.slx;
0352                 else
0353                     lambda.lower = [];
0354                 end
0355                 if isfield(sol, 'sux')
0356                     lambda.upper = sol.sux;
0357                 else
0358                     lambda.upper = [];
0359                 end
0360                 if isfield(sol, 'slc')
0361                     lambda.mu_l  = sol.slc;
0362                 else
0363                     lambda.mu_l  = [];
0364                 end
0365                 if isfield(sol, 'suc')
0366                     lambda.mu_u  = sol.suc;
0367                 else
0368                     lambda.mu_u  = [];
0369                 end
0370             else
0371                 if isfield(p, 'xmin') && ~isempty(p.xmin)
0372                     lambda.lower = NaN(nx, 1);
0373                 else
0374                     lambda.lower = [];
0375                 end
0376                 if isfield(p, 'xmax') && ~isempty(p.xmax)
0377                     lambda.upper = NaN(nx, 1);
0378                 else
0379                     lambda.upper = [];
0380                 end
0381                 lambda.mu_l = NaN(nA, 1);
0382                 lambda.mu_u = NaN(nA, 1);
0383             end
0384             if unconstrained
0385                 lambda.mu_l  = [];
0386                 lambda.mu_u  = [];
0387             end
0388         end
0389     end
0390 end

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005