Home > matpower5.0 > qps_glpk.m

qps_glpk

PURPOSE ^

QPS_GLPK Quadratic Program Solver based on QUADPROG/LINPROG.

SYNOPSIS ^

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

DESCRIPTION ^

QPS_GLPK  Quadratic Program Solver based on QUADPROG/LINPROG.
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
       QPS_GLPK(H, C, A, L, U, XMIN, XMAX, X0, OPT)
   A wrapper function providing a MATPOWER standardized interface for using
   GLKP to solve the following LP (linear programming) problem:

       min 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 : dummy matrix (possibly sparse) of quadratic cost coefficients
           for QP problems, which GLPK does not handle
       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 (NOT USED)
       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
           glpk_opt - options struct for GLPK, 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 - optimal, <= 0 - infeasible, unbounded or other
       OUTPUT : struct with errnum and status fields from GLPK output args
       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 GLPK. 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_glpk(H, c, A, l, u, xmin, xmax, x0, opt)

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


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

   See also GLPK.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x, f, eflag, output, lambda] = qps_glpk(H, c, A, l, u, xmin, xmax, x0, opt)
0002 %QPS_GLPK  Quadratic Program Solver based on QUADPROG/LINPROG.
0003 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
0004 %       QPS_GLPK(H, C, A, L, U, XMIN, XMAX, X0, OPT)
0005 %   A wrapper function providing a MATPOWER standardized interface for using
0006 %   GLKP to solve the following LP (linear programming) problem:
0007 %
0008 %       min 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 : dummy matrix (possibly sparse) of quadratic cost coefficients
0018 %           for QP problems, which GLPK does not handle
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 (NOT USED)
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 %           glpk_opt - options struct for GLPK, value in
0034 %               verbose 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, 1 - optimal, <= 0 - infeasible, unbounded or other
0043 %       OUTPUT : struct with errnum and status fields from GLPK output args
0044 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0045 %           multipliers on the constraints, with fields:
0046 %           mu_l - lower (left-hand) limit on linear constraints
0047 %           mu_u - upper (right-hand) limit on linear constraints
0048 %           lower - lower bound on optimization variables
0049 %           upper - upper bound on optimization variables
0050 %
0051 %   Note the calling syntax is almost identical to that of GLPK. The main
0052 %   difference is that the linear constraints are specified with A, L, U
0053 %   instead of A, B, Aeq, Beq.
0054 %
0055 %   Calling syntax options:
0056 %       [x, f, exitflag, output, lambda] = ...
0057 %           qps_glpk(H, c, A, l, u, xmin, xmax, x0, opt)
0058 %
0059 %       x = qps_glpk(H, c, A, l, u)
0060 %       x = qps_glpk(H, c, A, l, u, xmin, xmax)
0061 %       x = qps_glpk(H, c, A, l, u, xmin, xmax, x0)
0062 %       x = qps_glpk(H, c, A, l, u, xmin, xmax, x0, opt)
0063 %       x = qps_glpk(problem), where problem is a struct with fields:
0064 %                       H, c, A, l, u, xmin, xmax, x0, opt
0065 %                       all fields except 'c', 'A' and 'l' or 'u' are optional
0066 %       x = qps_glpk(...)
0067 %       [x, f] = qps_glpk(...)
0068 %       [x, f, exitflag] = qps_glpk(...)
0069 %       [x, f, exitflag, output] = qps_glpk(...)
0070 %       [x, f, exitflag, output, lambda] = qps_glpk(...)
0071 %
0072 %
0073 %   Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm)
0074 %       H = [   1003.1  4.3     6.3     5.9;
0075 %               4.3     2.2     2.1     3.9;
0076 %               6.3     2.1     3.5     4.8;
0077 %               5.9     3.9     4.8     10  ];
0078 %       c = zeros(4,1);
0079 %       A = [   1       1       1       1;
0080 %               0.17    0.11    0.10    0.18    ];
0081 %       l = [1; 0.10];
0082 %       u = [1; Inf];
0083 %       xmin = zeros(4,1);
0084 %       x0 = [1; 0; 0; 1];
0085 %       opt = struct('verbose', 2);
0086 %       [x, f, s, out, lambda] = qps_glpk(H, c, A, l, u, xmin, [], x0, opt);
0087 %
0088 %   See also GLPK.
0089 
0090 %   MATPOWER
0091 %   $Id: qps_glpk.m 2497 2014-12-17 19:56:20Z ray $
0092 %   by Ray Zimmerman, PSERC Cornell
0093 %   Copyright (c) 2010-2014 by Power System Engineering Research Center (PSERC)
0094 %
0095 %   This file is part of MATPOWER.
0096 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0097 %
0098 %   MATPOWER is free software: you can redistribute it and/or modify
0099 %   it under the terms of the GNU General Public License as published
0100 %   by the Free Software Foundation, either version 3 of the License,
0101 %   or (at your option) any later version.
0102 %
0103 %   MATPOWER is distributed in the hope that it will be useful,
0104 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0105 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0106 %   GNU General Public License for more details.
0107 %
0108 %   You should have received a copy of the GNU General Public License
0109 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0110 %
0111 %   Additional permission under GNU GPL version 3 section 7
0112 %
0113 %   If you modify MATPOWER, or any covered work, to interface with
0114 %   other modules (such as MATLAB code and MEX-files) available in a
0115 %   MATLAB(R) or comparable environment containing parts covered
0116 %   under other licensing terms, the licensors of MATPOWER grant
0117 %   you additional permission to convey the resulting work.
0118 
0119 %% check for Optimization Toolbox
0120 % if ~have_fcn('quadprog')
0121 %     error('qps_glpk: requires the MEX interface to GLPK');
0122 % end
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 %% define nx, set default values for missing optional inputs
0153 if isempty(H) || ~any(any(H))
0154     if isempty(A) && isempty(xmin) && isempty(xmax)
0155         error('qps_glpk: LP problem must include constraints or variable bounds');
0156     else
0157         if ~isempty(A)
0158             nx = size(A, 2);
0159         elseif ~isempty(xmin)
0160             nx = length(xmin);
0161         else    % if ~isempty(xmax)
0162             nx = length(xmax);
0163         end
0164     end
0165 else
0166     error('qps_glpk: GLPK handles only LP problems, not QP problems');
0167     nx = size(H, 1);
0168 end
0169 if isempty(c)
0170     c = zeros(nx, 1);
0171 end
0172 if  ~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0173                    (isempty(u) || all(u == Inf))
0174     A = sparse(0,nx);           %% no limits => no linear constraints
0175 end
0176 nA = size(A, 1);                %% number of original linear constraints
0177 if isempty(u)                   %% By default, linear inequalities are ...
0178     u = Inf * ones(nA, 1);      %% ... unbounded above and ...
0179 end
0180 if isempty(l)
0181     l = -Inf * ones(nA, 1);     %% ... unbounded below.
0182 end
0183 if isempty(xmin)                %% By default, optimization variables are ...
0184     xmin = -Inf * ones(nx, 1);  %% ... unbounded below and ...
0185 end
0186 if isempty(xmax)
0187     xmax = Inf * ones(nx, 1);   %% ... unbounded above.
0188 end
0189 if isempty(x0)
0190     x0 = zeros(nx, 1);
0191 end
0192 
0193 %% default options
0194 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0195     verbose = opt.verbose;
0196 else
0197     verbose = 0;
0198 end
0199 
0200 %% split up linear constraints
0201 ieq = find( abs(u-l) <= eps );          %% equality
0202 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0203 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0204 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0205 AA = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0206 bb = [ u(ieq);    u(ilt);    -l(igt);    u(ibx);    -l(ibx)];
0207 
0208 %% grab some dimensions
0209 nlt = length(ilt);      %% number of upper bounded linear inequalities
0210 ngt = length(igt);      %% number of lower bounded linear inequalities
0211 nbx = length(ibx);      %% number of doubly bounded linear inequalities
0212 neq = length(ieq);      %% number of equalities
0213 nie = nlt+ngt+2*nbx;    %% number of inequalities
0214 
0215 ctype = [repmat('S', neq, 1); repmat('U', nlt+ngt+2*nbx, 1)];
0216 vtype = repmat('C', nx, 1);
0217 
0218 %% set options struct for GLPK
0219 if ~isempty(opt) && isfield(opt, 'glpk_opt') && ~isempty(opt.glpk_opt)
0220     glpk_opt = glpk_options(opt.glpk_opt);
0221 else
0222     glpk_opt = glpk_options;
0223 end
0224 glpk_opt.msglev = verbose;
0225 
0226 %% call the solver
0227 [x, f, errnum, extra] = ...
0228     glpk(c, AA, bb, xmin, xmax, ctype, vtype, 1, glpk_opt);
0229 
0230 %% set exit flag
0231 if isfield(extra, 'status')
0232     output.errnum = errnum;
0233     output.status = extra.status;
0234     eflag = -errnum;
0235     if eflag == 0 && extra.status == 5
0236         eflag = 1;
0237     end
0238 else
0239     output.errnum = [];
0240     output.status = errnum;
0241     if have_fcn('octave')
0242         if errnum == 180 || errnum == 151
0243             eflag = 1;
0244         else
0245             eflag = 0;
0246         end
0247     else
0248         if errnum == 5
0249             eflag = 1;
0250         else
0251             eflag = 0;
0252         end
0253     end
0254 end
0255 
0256 %% repackage lambdas
0257 if isempty(extra) || ~isfield(extra, 'lambda') || isempty(extra.lambda)
0258     lambda = struct( ...
0259         'mu_l', zeros(nA, 1), ...
0260         'mu_u', zeros(nA, 1), ...
0261         'lower', zeros(nx, 1), ...
0262         'upper', zeros(nx, 1) ...
0263     );
0264 else
0265     lam.eqlin = extra.lambda(1:neq);
0266     lam.ineqlin = extra.lambda(neq+(1:nie));
0267     lam.lower = extra.redcosts;
0268     lam.upper = -extra.redcosts;
0269     lam.lower(lam.lower < 0) = 0;
0270     lam.upper(lam.upper < 0) = 0;
0271     kl = find(lam.eqlin > 0);   %% lower bound binding
0272     ku = find(lam.eqlin < 0);   %% upper bound binding
0273 
0274     mu_l = zeros(nA, 1);
0275     mu_l(ieq(kl)) = lam.eqlin(kl);
0276     mu_l(igt) = lam.ineqlin(nlt+(1:ngt));
0277     mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0278 
0279     mu_u = zeros(nA, 1);
0280     mu_u(ieq(ku)) = -lam.eqlin(ku);
0281     mu_u(ilt) = -lam.ineqlin(1:nlt);
0282     mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx));
0283 
0284     lambda = struct( ...
0285         'mu_l', mu_l, ...
0286         'mu_u', mu_u, ...
0287         'lower', lam.lower(1:nx), ...
0288         'upper', lam.upper(1:nx) ...
0289     );
0290 end

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