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

Generated on Mon 24-Jun-2019 15:58:45 by m2html © 2005