Home > matpower6.0 > qps_matpower.m

qps_matpower

PURPOSE ^

QPS_MATPOWER Quadratic Program Solver for MATPOWER.

SYNOPSIS ^

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

DESCRIPTION ^

QPS_MATPOWER  Quadratic Program Solver for MATPOWER.
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
       QPS_MATPOWER(H, C, A, L, U, XMIN, XMAX, X0, OPT)
   A common wrapper function for various QP solvers. 
   Solves 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)
           alg ('DEFAULT') : determines which solver to use, can be either
                   a (new-style) string or an (old-style) numerical alg code
               'DEFAULT' : (or 0) automatic, first available of CPLEX,
                       Gurobi, MOSEK, Opt Tbx (if Matlab), GLPK (LPs only),
                       BPMPD, MIPS
               'MIPS'    : (or 200) MIPS, MATPOWER Interior Point Solver
                        pure MATLAB implementation of a primal-dual
                        interior point method, if mips_opt.step_control = 1
                        (or alg=250) it uses MIPS-sc, a step controlled
                        variant of MIPS
               'BPMPD'   : (or 100) BPMPD_MEX
               'CLP'     : CLP
               'CPLEX'   : (or 500) CPLEX
               'GLPK'    : GLPK, (LP problems only, i.e. empty H matrix)
               'GUROBI'  : (or 700) Gurobi
               'IPOPT'   : (or 400) IPOPT
               'MOSEK'   : (or 600) MOSEK
               'OT'      : (or 300) Optimization Toolbox, QUADPROG or LINPROG
           verbose (0) - controls level of progress output displayed
               0 = no progress output
               1 = some progress output
               2 = verbose progress output
           bp_opt      - options vector for BP
           cplex_opt   - options struct for CPLEX
           glpk_opt    - options struct for GLPK
           grb_opt     - options struct for GBUROBI_MEX
           ipopt_opt   - options struct for IPOPT
           linprog_opt - options struct for LINPROG
           mips_opt    - options struct for QPS_MIPS
           mosek_opt   - options struct for MOSEK
           quadprog_opt - options struct for QUADPROG
       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 = converged
           0 or negative values = algorithm specific failure codes
       OUTPUT : output struct with the following fields:
           alg - algorithm code of solver used
           (others) - algorithm specific fields
       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_matpower(H, c, A, l, u, xmin, xmax, x0, opt)

       x = qps_matpower(H, c, A, l, u)
       x = qps_matpower(H, c, A, l, u, xmin, xmax)
       x = qps_matpower(H, c, A, l, u, xmin, xmax, x0)
       x = qps_matpower(H, c, A, l, u, xmin, xmax, x0, opt)
       x = qps_matpower(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_matpower(...)
       [x, f] = qps_matpower(...)
       [x, f, exitflag] = qps_matpower(...)
       [x, f, exitflag, output] = qps_matpower(...)
       [x, f, exitflag, output, lambda] = qps_matpower(...)

   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_matpower(H, c, A, l, u, xmin, [], x0, opt);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x, f, eflag, output, lambda] = qps_matpower(H, c, A, l, u, xmin, xmax, x0, opt)
0002 %QPS_MATPOWER  Quadratic Program Solver for MATPOWER.
0003 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
0004 %       QPS_MATPOWER(H, C, A, L, U, XMIN, XMAX, X0, OPT)
0005 %   A common wrapper function for various QP solvers.
0006 %   Solves 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 %           alg ('DEFAULT') : determines which solver to use, can be either
0029 %                   a (new-style) string or an (old-style) numerical alg code
0030 %               'DEFAULT' : (or 0) automatic, first available of CPLEX,
0031 %                       Gurobi, MOSEK, Opt Tbx (if Matlab), GLPK (LPs only),
0032 %                       BPMPD, MIPS
0033 %               'MIPS'    : (or 200) MIPS, MATPOWER Interior Point Solver
0034 %                        pure MATLAB implementation of a primal-dual
0035 %                        interior point method, if mips_opt.step_control = 1
0036 %                        (or alg=250) it uses MIPS-sc, a step controlled
0037 %                        variant of MIPS
0038 %               'BPMPD'   : (or 100) BPMPD_MEX
0039 %               'CLP'     : CLP
0040 %               'CPLEX'   : (or 500) CPLEX
0041 %               'GLPK'    : GLPK, (LP problems only, i.e. empty H matrix)
0042 %               'GUROBI'  : (or 700) Gurobi
0043 %               'IPOPT'   : (or 400) IPOPT
0044 %               'MOSEK'   : (or 600) MOSEK
0045 %               'OT'      : (or 300) Optimization Toolbox, QUADPROG or LINPROG
0046 %           verbose (0) - controls level of progress output displayed
0047 %               0 = no progress output
0048 %               1 = some progress output
0049 %               2 = verbose progress output
0050 %           bp_opt      - options vector for BP
0051 %           cplex_opt   - options struct for CPLEX
0052 %           glpk_opt    - options struct for GLPK
0053 %           grb_opt     - options struct for GBUROBI_MEX
0054 %           ipopt_opt   - options struct for IPOPT
0055 %           linprog_opt - options struct for LINPROG
0056 %           mips_opt    - options struct for QPS_MIPS
0057 %           mosek_opt   - options struct for MOSEK
0058 %           quadprog_opt - options struct for QUADPROG
0059 %       PROBLEM : The inputs can alternatively be supplied in a single
0060 %           PROBLEM struct with fields corresponding to the input arguments
0061 %           described above: H, c, A, l, u, xmin, xmax, x0, opt
0062 %
0063 %   Outputs:
0064 %       X : solution vector
0065 %       F : final objective function value
0066 %       EXITFLAG : exit flag
0067 %           1 = converged
0068 %           0 or negative values = algorithm specific failure codes
0069 %       OUTPUT : output struct with the following fields:
0070 %           alg - algorithm code of solver used
0071 %           (others) - algorithm specific fields
0072 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0073 %           multipliers on the constraints, with fields:
0074 %           mu_l - lower (left-hand) limit on linear constraints
0075 %           mu_u - upper (right-hand) limit on linear constraints
0076 %           lower - lower bound on optimization variables
0077 %           upper - upper bound on optimization variables
0078 %
0079 %   Note the calling syntax is almost identical to that of QUADPROG
0080 %   from MathWorks' Optimization Toolbox. The main difference is that
0081 %   the linear constraints are specified with A, L, U instead of
0082 %   A, B, Aeq, Beq.
0083 %
0084 %   Calling syntax options:
0085 %       [x, f, exitflag, output, lambda] = ...
0086 %           qps_matpower(H, c, A, l, u, xmin, xmax, x0, opt)
0087 %
0088 %       x = qps_matpower(H, c, A, l, u)
0089 %       x = qps_matpower(H, c, A, l, u, xmin, xmax)
0090 %       x = qps_matpower(H, c, A, l, u, xmin, xmax, x0)
0091 %       x = qps_matpower(H, c, A, l, u, xmin, xmax, x0, opt)
0092 %       x = qps_matpower(problem), where problem is a struct with fields:
0093 %                       H, c, A, l, u, xmin, xmax, x0, opt
0094 %                       all fields except 'c', 'A' and 'l' or 'u' are optional
0095 %       x = qps_matpower(...)
0096 %       [x, f] = qps_matpower(...)
0097 %       [x, f, exitflag] = qps_matpower(...)
0098 %       [x, f, exitflag, output] = qps_matpower(...)
0099 %       [x, f, exitflag, output, lambda] = qps_matpower(...)
0100 %
0101 %   Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm)
0102 %       H = [   1003.1  4.3     6.3     5.9;
0103 %               4.3     2.2     2.1     3.9;
0104 %               6.3     2.1     3.5     4.8;
0105 %               5.9     3.9     4.8     10  ];
0106 %       c = zeros(4,1);
0107 %       A = [   1       1       1       1;
0108 %               0.17    0.11    0.10    0.18    ];
0109 %       l = [1; 0.10];
0110 %       u = [1; Inf];
0111 %       xmin = zeros(4,1);
0112 %       x0 = [1; 0; 0; 1];
0113 %       opt = struct('verbose', 2);
0114 %       [x, f, s, out, lambda] = qps_matpower(H, c, A, l, u, xmin, [], x0, opt);
0115 
0116 %   MATPOWER
0117 %   Copyright (c) 2010-2016, Power Systems Engineering Research Center (PSERC)
0118 %   by Ray Zimmerman, PSERC Cornell
0119 %
0120 %   This file is part of MATPOWER.
0121 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0122 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0123 
0124 %%----- input argument handling  -----
0125 %% gather inputs
0126 if nargin == 1 && isstruct(H)       %% problem struct
0127     p = H;
0128     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0129     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0130     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0131     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0132     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0133     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0134     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0135     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0136     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0137 else                                %% individual args
0138     if nargin < 9
0139         opt = [];
0140         if nargin < 8
0141             x0 = [];
0142             if nargin < 7
0143                 xmax = [];
0144                 if nargin < 6
0145                     xmin = [];
0146                 end
0147             end
0148         end
0149     end
0150 end
0151 
0152 %% default options
0153 if ~isempty(opt) && isfield(opt, 'alg') && ~isempty(opt.alg)
0154     alg = opt.alg;
0155     %% convert integer codes to string values
0156     if ~ischar(alg)
0157         switch alg
0158             case 0
0159                 alg = 'DEFAULT';
0160             case 100
0161                 alg = 'BPMPD';
0162             case 200
0163                 alg = 'MIPS';
0164                 opt.mips_opt.step_control = 0;
0165             case 250
0166                 alg = 'MIPS';
0167                 opt.mips_opt.step_control = 1;
0168             case 300
0169                 alg = 'OT';
0170             case 400
0171                 alg = 'IPOPT';
0172             case 500
0173                 alg = 'CPLEX';
0174             case 600
0175                 alg = 'MOSEK';
0176             case 700
0177                 alg = 'GUROBI';
0178             otherwise
0179                 error('qps_matpower: %d is not a valid algorithm code', alg);
0180         end
0181     end
0182 else
0183     alg = 'DEFAULT';
0184 end
0185 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0186     verbose = opt.verbose;
0187 else
0188     verbose = 0;
0189 end
0190 if strcmp(alg, 'DEFAULT')
0191     if have_fcn('gurobi')       %% use Gurobi by default, if available
0192         alg = 'GUROBI';
0193     elseif have_fcn('cplex')    %% if not, then CPLEX, if available
0194         alg = 'CPLEX';
0195     elseif have_fcn('mosek')    %% if not, then MOSEK, if available
0196         alg = 'MOSEK';
0197     elseif have_fcn('quadprog') && have_fcn('matlab')   %% if not, then Opt Tbx, if available in Matlab
0198         alg = 'OT';
0199     elseif (isempty(H) || ~any(any(H))) && have_fcn('glpk') %% if not, and
0200         alg = 'GLPK';           %% prob is LP (not QP), then GLPK, if available
0201     elseif have_fcn('bpmpd')    %% if not, then BPMPD_MEX, if available
0202         alg = 'BPMPD';
0203     else                        %% otherwise MIPS
0204         alg = 'MIPS';
0205     end
0206 end
0207 
0208 %%----- call the appropriate solver  -----
0209 switch alg
0210     case 'BPMPD'                    %% use BPMPD_MEX
0211         [x, f, eflag, output, lambda] = ...
0212             qps_bpmpd(H, c, A, l, u, xmin, xmax, x0, opt);
0213 
0214         if eflag == -99
0215             if verbose
0216                 fprintf('         Retrying with QPS_MIPS solver ...\n\n');
0217             end
0218             %% save (incorrect) solution from BPMPD
0219             bpmpd = struct('x', x, 'f', f, 'eflag', eflag, ...
0220                             'output', output, 'lambda', lambda);
0221             opt.alg = 'MIPS';
0222             [x, f, eflag, output, lambda] = ...
0223                 qps_matpower(H, c, A, l, u, xmin, xmax, x0, opt);
0224             output.bpmpd = bpmpd;
0225         end
0226     case 'CLP'
0227         [x, f, eflag, output, lambda] = ...
0228             qps_clp(H, c, A, l, u, xmin, xmax, x0, opt);
0229     case 'CPLEX'
0230         [x, f, eflag, output, lambda] = ...
0231             qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt);
0232     case 'GLPK'
0233         [x, f, eflag, output, lambda] = ...
0234             qps_glpk(H, c, A, l, u, xmin, xmax, x0, opt);
0235     case 'GUROBI'
0236         [x, f, eflag, output, lambda] = ...
0237             qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt);
0238     case 'IPOPT'
0239         [x, f, eflag, output, lambda] = ...
0240             qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt);
0241     case 'MIPS'
0242         %% set up options
0243         if ~isempty(opt) && isfield(opt, 'mips_opt') && ~isempty(opt.mips_opt)
0244             mips_opt = opt.mips_opt;
0245         else
0246             mips_opt = [];
0247         end
0248         mips_opt.verbose = verbose;
0249         
0250         %% call solver
0251         [x, f, eflag, output, lambda] = ...
0252             qps_mips(H, c, A, l, u, xmin, xmax, x0, mips_opt);
0253     case 'MOSEK'
0254         [x, f, eflag, output, lambda] = ...
0255             qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt);
0256     case 'OT'                   %% use QUADPROG or LINPROG from Opt Tbx ver 2.x+
0257         [x, f, eflag, output, lambda] = ...
0258             qps_ot(H, c, A, l, u, xmin, xmax, x0, opt);
0259     otherwise
0260         error('qps_matpower: ''%s'' is not a valid algorithm code', alg);
0261 end
0262 if ~isfield(output, 'alg') || isempty(output.alg)
0263     output.alg = alg;
0264 end

Generated on Fri 16-Dec-2016 12:45:37 by m2html © 2005