Home > matpower6.0 > 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-2016, Power Systems Engineering Research Center (PSERC)
0097 %   by Ray Zimmerman, PSERC Cornell
0098 %
0099 %   This file is part of MATPOWER.
0100 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0101 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0102 
0103 %% check for Optimization Toolbox
0104 % if ~have_fcn('quadprog')
0105 %     error('qps_ot: requires the Optimization Toolbox');
0106 % end
0107 
0108 %%----- input argument handling  -----
0109 %% gather inputs
0110 if nargin == 1 && isstruct(H)       %% problem struct
0111     p = H;
0112     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0113     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0114     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0115     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0116     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0117     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0118     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0119     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0120     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0121 else                                %% individual args
0122     if nargin < 9
0123         opt = [];
0124         if nargin < 8
0125             x0 = [];
0126             if nargin < 7
0127                 xmax = [];
0128                 if nargin < 6
0129                     xmin = [];
0130                 end
0131             end
0132         end
0133     end
0134 end
0135 
0136 %% define nx, set default values for missing optional inputs
0137 if isempty(H) || ~any(any(H))
0138     if isempty(A) && isempty(xmin) && isempty(xmax)
0139         error('qps_ot: LP problem must include constraints or variable bounds');
0140     else
0141         if ~isempty(A)
0142             nx = size(A, 2);
0143         elseif ~isempty(xmin)
0144             nx = length(xmin);
0145         else    % if ~isempty(xmax)
0146             nx = length(xmax);
0147         end
0148     end
0149 else
0150     nx = size(H, 1);
0151 end
0152 if isempty(c)
0153     c = zeros(nx, 1);
0154 end
0155 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0156                                  (isempty(u) || all(u == Inf)))
0157     A = sparse(0,nx);           %% no limits => no linear constraints
0158 end
0159 nA = size(A, 1);                %% number of original linear constraints
0160 if isempty(u)                   %% By default, linear inequalities are ...
0161     u = Inf(nA, 1);             %% ... unbounded above and ...
0162 end
0163 if isempty(l)
0164     l = -Inf(nA, 1);            %% ... unbounded below.
0165 end
0166 if isempty(xmin)                %% By default, optimization variables are ...
0167     xmin = -Inf(nx, 1);         %% ... unbounded below and ...
0168 end
0169 if isempty(xmax)
0170     xmax = Inf(nx, 1);          %% ... unbounded above.
0171 end
0172 if isempty(x0)
0173     x0 = zeros(nx, 1);
0174 end
0175 if isempty(H) || ~any(any(H))
0176     isLP = 1;   %% it's an LP
0177 else
0178     isLP = 0;   %% nope, it's a QP
0179 end
0180 
0181 %% default options
0182 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0183     verbose = opt.verbose;
0184 else
0185     verbose = 0;
0186 end
0187 %% Matlab or Octave
0188 matlab = have_fcn('matlab');
0189 otver = have_fcn('quadprog', 'vnum');
0190 
0191 %% split up linear constraints
0192 ieq = find( abs(u-l) <= eps );          %% equality
0193 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0194 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0195 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0196 Ae = A(ieq, :);
0197 be = u(ieq);
0198 Ai  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0199 bi  = [ u(ilt);    -l(igt);    u(ibx);    -l(ibx)];
0200 
0201 %% grab some dimensions
0202 nlt = length(ilt);      %% number of upper bounded linear inequalities
0203 ngt = length(igt);      %% number of lower bounded linear inequalities
0204 nbx = length(ibx);      %% number of doubly bounded linear inequalities
0205 
0206 %% set up options
0207 if verbose > 1
0208     vrb = 'iter';       %% seems to be same as 'final'
0209 elseif verbose == 1
0210     vrb = 'final';
0211 else
0212     vrb = 'off';
0213 end
0214 if have_fcn('optimoptions')     %% Optimization Tbx 6.3 + (R2013a +)
0215     %% could use optimset for everything, except some options are not
0216     %% recognized by optimset, only optimoptions, such as
0217     %% ot_opt.Algorithm = 'dual-simplex'
0218     if isLP
0219         ot_opt = optimoptions('linprog');
0220         if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt)
0221             ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt);
0222         end
0223     else
0224         ot_opt = optimoptions('quadprog');
0225         if have_fcn('quadprog_ls')
0226             ot_opt.Algorithm = 'interior-point-convex';
0227         else
0228             ot_opt.LargeScale = 'off';
0229         end
0230         if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt)
0231             ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt);
0232         end
0233     end
0234     ot_opt = optimoptions(ot_opt, 'Display', vrb);
0235 else                            %% need to use optimset()
0236     if isLP
0237         if matlab
0238             ot_opt = optimset('linprog');
0239         else
0240             ot_opt = optimset();
0241         end
0242         if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt)
0243             ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt);
0244         end
0245     else
0246         if matlab
0247             ot_opt = optimset('quadprog');
0248             if have_fcn('quadprog_ls')
0249                 ot_opt = optimset(ot_opt, 'Algorithm', 'interior-point-convex');
0250             else
0251                 ot_opt = optimset(ot_opt, 'LargeScale', 'off');
0252             end
0253         else
0254             ot_opt = optimset();
0255         end
0256         if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt)
0257             ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt);
0258         end
0259     end
0260     ot_opt = optimset(ot_opt, 'Display', vrb);
0261 end
0262 
0263 %% call the solver
0264 if isLP
0265     if matlab
0266         [x, f, eflag, output, lam] = ...
0267             linprog(c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt);
0268     else
0269 % don't use linprog under Octave (using GLPK directly is recommended)
0270 %         [x, f] = linprog(c, Ai, bi, Ae, be, xmin, xmax);
0271 %         eflag = [];
0272 %         output = [];
0273 %         lam = [];
0274         [x, f, eflag, output, lam] = ...
0275             quadprog(sparse(nx,nx), c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt);
0276     end
0277 else
0278     [x, f, eflag, output, lam] = ...
0279         quadprog(H, c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt);
0280 end
0281 
0282 %% repackage lambdas
0283 if isempty(x)
0284     x = NaN(nx, 1);
0285 end
0286 if isempty(lam) || (isempty(lam.eqlin) && isempty(lam.ineqlin) && ...
0287                     isempty(lam.lower) && isempty(lam.upper))
0288     lambda = struct( ...
0289         'mu_l', NaN(nA, 1), ...
0290         'mu_u', NaN(nA, 1), ...
0291         'lower', NaN(nx, 1), ...
0292         'upper', NaN(nx, 1) ...
0293     );
0294 else
0295     kl = find(lam.eqlin < 0);   %% lower bound binding
0296     ku = find(lam.eqlin > 0);   %% upper bound binding
0297 
0298     mu_l = zeros(nA, 1);
0299 %     %% workaround for Octave optim 1.5.0 and earlier, which
0300 %     %% has opposite sign convention for equality multipliers
0301 %     if ~matlab && otver <= 1.005
0302 %         mu_l(ieq(ku)) = lam.eqlin(ku);
0303 %     else
0304         mu_l(ieq(kl)) = -lam.eqlin(kl);
0305 %     end
0306     mu_l(igt) = lam.ineqlin(nlt+(1:ngt));
0307     mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0308 
0309     mu_u = zeros(nA, 1);
0310 %     %% workaround for Octave optim 1.5.0 and earlier, which
0311 %     %% has opposite sign convention for equality multipliers
0312 %     if ~matlab && otver <= 1.005
0313 %         mu_u(ieq(kl)) = -lam.eqlin(kl);
0314 %     else
0315         mu_u(ieq(ku)) = lam.eqlin(ku);
0316 %     end
0317     mu_u(ilt) = lam.ineqlin(1:nlt);
0318     mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx));
0319 
0320     %% workaround for Octave optim 1.5.0 and earlier, which
0321     %% has opposite sign convention for equality multipliers
0322     % if ~matlab && otver <= 1.005
0323         %% there are also issues with variable bounds that are
0324         %% converted to equalities, and maybe other issues
0325     % end
0326 
0327     lambda = struct( ...
0328         'mu_l', mu_l, ...
0329         'mu_u', mu_u, ...
0330         'lower', lam.lower(1:nx), ...
0331         'upper', lam.upper(1:nx) ...
0332     );
0333 end

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