Home > matpower5.1 > qps_ot.m

qps_ot

PURPOSE ^

QPS_OT Quadratic Program Solver based on QUADPROG/LINPROG.

SYNOPSIS ^

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

DESCRIPTION ^

QPS_OT  Quadratic Program Solver based on QUADPROG/LINPROG.
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
       QPS_OT(H, C, A, L, U, XMIN, XMAX, X0, OPT)
   A wrapper function providing a MATPOWER standardized interface for using
   QUADPROG or LINPROG from the Optimization Toolbox 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
           linprog_opt - options struct for LINPROG, value in
               verbose overrides these options
           quadprog_opt - options struct for QUADPROG, 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 : QUADPROG/LINPROG exit flag
           (see QUADPROG and LINPROG documentation for details)
       OUTPUT : QUADPROG/LINPROG output struct
           (see QUADPROG and LINPROG documentation for details)
       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_ot(H, c, A, l, u, xmin, xmax, x0, opt)

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


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

   See also QUADPROG, LINPROG.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x, f, eflag, output, lambda] = qps_ot(H, c, A, l, u, xmin, xmax, x0, opt)
0002 %QPS_OT  Quadratic Program Solver based on QUADPROG/LINPROG.
0003 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
0004 %       QPS_OT(H, C, A, L, U, XMIN, XMAX, X0, OPT)
0005 %   A wrapper function providing a MATPOWER standardized interface for using
0006 %   QUADPROG or LINPROG from the Optimization Toolbox to solve the
0007 %   following QP (quadratic programming) problem:
0008 %
0009 %       min 1/2 X'*H*X + C'*X
0010 %        X
0011 %
0012 %   subject to
0013 %
0014 %       L <= A*X <= U       (linear constraints)
0015 %       XMIN <= X <= XMAX   (variable bounds)
0016 %
0017 %   Inputs (all optional except H, C, A and L):
0018 %       H : matrix (possibly sparse) of quadratic cost coefficients
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
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 %           linprog_opt - options struct for LINPROG, value in
0034 %               verbose overrides these options
0035 %           quadprog_opt - options struct for QUADPROG, value in
0036 %               verbose overrides these options
0037 %       PROBLEM : The inputs can alternatively be supplied in a single
0038 %           PROBLEM struct with fields corresponding to the input arguments
0039 %           described above: H, c, A, l, u, xmin, xmax, x0, opt
0040 %
0041 %   Outputs:
0042 %       X : solution vector
0043 %       F : final objective function value
0044 %       EXITFLAG : QUADPROG/LINPROG exit flag
0045 %           (see QUADPROG and LINPROG documentation for details)
0046 %       OUTPUT : QUADPROG/LINPROG output struct
0047 %           (see QUADPROG and LINPROG documentation for details)
0048 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0049 %           multipliers on the constraints, with fields:
0050 %           mu_l - lower (left-hand) limit on linear constraints
0051 %           mu_u - upper (right-hand) limit on linear constraints
0052 %           lower - lower bound on optimization variables
0053 %           upper - upper bound on optimization variables
0054 %
0055 %   Note the calling syntax is almost identical to that of QUADPROG
0056 %   from MathWorks' Optimization Toolbox. The main difference is that
0057 %   the linear constraints are specified with A, L, U instead of
0058 %   A, B, Aeq, Beq.
0059 %
0060 %   Calling syntax options:
0061 %       [x, f, exitflag, output, lambda] = ...
0062 %           qps_ot(H, c, A, l, u, xmin, xmax, x0, opt)
0063 %
0064 %       x = qps_ot(H, c, A, l, u)
0065 %       x = qps_ot(H, c, A, l, u, xmin, xmax)
0066 %       x = qps_ot(H, c, A, l, u, xmin, xmax, x0)
0067 %       x = qps_ot(H, c, A, l, u, xmin, xmax, x0, opt)
0068 %       x = qps_ot(problem), where problem is a struct with fields:
0069 %                       H, c, A, l, u, xmin, xmax, x0, opt
0070 %                       all fields except 'c', 'A' and 'l' or 'u' are optional
0071 %       x = qps_ot(...)
0072 %       [x, f] = qps_ot(...)
0073 %       [x, f, exitflag] = qps_ot(...)
0074 %       [x, f, exitflag, output] = qps_ot(...)
0075 %       [x, f, exitflag, output, lambda] = qps_ot(...)
0076 %
0077 %
0078 %   Example: (problem from from http://www.jmu.edu/docs/sasdoc/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_ot(H, c, A, l, u, xmin, [], x0, opt);
0092 %
0093 %   See also QUADPROG, LINPROG.
0094 
0095 %   MATPOWER
0096 %   Copyright (c) 2010-2015 by Power System Engineering Research Center (PSERC)
0097 %   by Ray Zimmerman, PSERC Cornell
0098 %
0099 %   $Id: qps_ot.m 2661 2015-03-20 17:02:46Z ray $
0100 %
0101 %   This file is part of MATPOWER.
0102 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0103 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0104 
0105 %% check for Optimization Toolbox
0106 % if ~have_fcn('quadprog')
0107 %     error('qps_ot: requires the Optimization Toolbox');
0108 % end
0109 
0110 %%----- input argument handling  -----
0111 %% gather inputs
0112 if nargin == 1 && isstruct(H)       %% problem struct
0113     p = H;
0114     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0115     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0116     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0117     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0118     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0119     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0120     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0121     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0122     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0123 else                                %% individual args
0124     if nargin < 9
0125         opt = [];
0126         if nargin < 8
0127             x0 = [];
0128             if nargin < 7
0129                 xmax = [];
0130                 if nargin < 6
0131                     xmin = [];
0132                 end
0133             end
0134         end
0135     end
0136 end
0137 
0138 %% define nx, set default values for missing optional inputs
0139 if isempty(H) || ~any(any(H))
0140     if isempty(A) && isempty(xmin) && isempty(xmax)
0141         error('qps_ot: LP problem must include constraints or variable bounds');
0142     else
0143         if ~isempty(A)
0144             nx = size(A, 2);
0145         elseif ~isempty(xmin)
0146             nx = length(xmin);
0147         else    % if ~isempty(xmax)
0148             nx = length(xmax);
0149         end
0150     end
0151 else
0152     nx = size(H, 1);
0153 end
0154 if isempty(c)
0155     c = zeros(nx, 1);
0156 end
0157 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0158                                  (isempty(u) || all(u == Inf)))
0159     A = sparse(0,nx);           %% no limits => no linear constraints
0160 end
0161 nA = size(A, 1);                %% number of original linear constraints
0162 if isempty(u)                   %% By default, linear inequalities are ...
0163     u = Inf(nA, 1);             %% ... unbounded above and ...
0164 end
0165 if isempty(l)
0166     l = -Inf(nA, 1);            %% ... unbounded below.
0167 end
0168 if isempty(xmin)                %% By default, optimization variables are ...
0169     xmin = -Inf(nx, 1);         %% ... unbounded below and ...
0170 end
0171 if isempty(xmax)
0172     xmax = Inf(nx, 1);          %% ... unbounded above.
0173 end
0174 if isempty(x0)
0175     x0 = zeros(nx, 1);
0176 end
0177 if isempty(H) || ~any(any(H))
0178     isLP = 1;   %% it's an LP
0179 else
0180     isLP = 0;   %% nope, it's a QP
0181 end
0182 
0183 %% default options
0184 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0185     verbose = opt.verbose;
0186 else
0187     verbose = 0;
0188 end
0189 
0190 %% split up linear constraints
0191 ieq = find( abs(u-l) <= eps );          %% equality
0192 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0193 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0194 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0195 Ae = A(ieq, :);
0196 be = u(ieq);
0197 Ai  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0198 bi  = [ u(ilt);    -l(igt);    u(ibx);    -l(ibx)];
0199 
0200 %% grab some dimensions
0201 nlt = length(ilt);      %% number of upper bounded linear inequalities
0202 ngt = length(igt);      %% number of lower bounded linear inequalities
0203 nbx = length(ibx);      %% number of doubly bounded linear inequalities
0204 
0205 %% set up options
0206 if verbose > 1
0207     vrb = 'iter';       %% seems to be same as 'final'
0208 elseif verbose == 1
0209     vrb = 'final';
0210 else
0211     vrb = 'off';
0212 end
0213 if have_fcn('optimoptions')     %% Optimization Tbx 6.3 + (R2013a +)
0214     %% could use optimset for everything, except some options are not
0215     %% recognized by optimset, only optimoptions, such as
0216     %% ot_opt.Algorithm = 'dual-simplex'
0217     if isLP
0218         ot_opt = optimoptions('linprog');
0219         if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt)
0220             ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt);
0221         end
0222     else
0223         ot_opt = optimoptions('quadprog');
0224         if have_fcn('quadprog_ls')
0225             ot_opt.Algorithm = 'interior-point-convex';
0226         else
0227             ot_opt.LargeScale = 'off';
0228         end
0229         if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt)
0230             ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt);
0231         end
0232     end
0233     ot_opt = optimoptions(ot_opt, 'Display', vrb);
0234 else                            %% need to use optimset()
0235     if isLP
0236         ot_opt = optimset('linprog');
0237         if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt)
0238             ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt);
0239         end
0240     else
0241         ot_opt = optimset('quadprog');
0242         if have_fcn('quadprog_ls')
0243             ot_opt = optimset(ot_opt, 'Algorithm', 'interior-point-convex');
0244         else
0245             ot_opt = optimset(ot_opt, 'LargeScale', 'off');
0246         end
0247         if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt)
0248             ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt);
0249         end
0250     end
0251     ot_opt = optimset(ot_opt, 'Display', vrb);
0252 end
0253 
0254 %% call the solver
0255 if isLP
0256     [x, f, eflag, output, lam] = ...
0257         linprog(c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt);
0258 else
0259     [x, f, eflag, output, lam] = ...
0260         quadprog(H, c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt);
0261 end
0262 
0263 %% repackage lambdas
0264 if isempty(x)
0265     x = NaN(nx, 1);
0266 end
0267 if isempty(lam) || (isempty(lam.eqlin) && isempty(lam.ineqlin) && ...
0268                     isempty(lam.lower) && isempty(lam.upper))
0269     lambda = struct( ...
0270         'mu_l', NaN(nA, 1), ...
0271         'mu_u', NaN(nA, 1), ...
0272         'lower', NaN(nx, 1), ...
0273         'upper', NaN(nx, 1) ...
0274     );
0275 else
0276     kl = find(lam.eqlin < 0);   %% lower bound binding
0277     ku = find(lam.eqlin > 0);   %% upper bound binding
0278 
0279     mu_l = zeros(nA, 1);
0280     mu_l(ieq(kl)) = -lam.eqlin(kl);
0281     mu_l(igt) = lam.ineqlin(nlt+(1:ngt));
0282     mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0283 
0284     mu_u = zeros(nA, 1);
0285     mu_u(ieq(ku)) = lam.eqlin(ku);
0286     mu_u(ilt) = lam.ineqlin(1:nlt);
0287     mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx));
0288 
0289     lambda = struct( ...
0290         'mu_l', mu_l, ...
0291         'mu_u', mu_u, ...
0292         'lower', lam.lower(1:nx), ...
0293         'upper', lam.upper(1:nx) ...
0294     );
0295 end

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