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

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005