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

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