Home > matpower7.1 > mp-opt-model > lib > 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)
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_OT(PROBLEM)
   A wrapper function providing a 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 https://v8doc.sas.com/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 QPS_MASTER, 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 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_OT(PROBLEM)
0006 %   A wrapper function providing a standardized interface for using
0007 %   QUADPROG or LINPROG from the Optimization Toolbox to solve the
0008 %   following QP (quadratic programming) problem:
0009 %
0010 %       min 1/2 X'*H*X + C'*X
0011 %        X
0012 %
0013 %   subject to
0014 %
0015 %       L <= A*X <= U       (linear constraints)
0016 %       XMIN <= X <= XMAX   (variable bounds)
0017 %
0018 %   Inputs (all optional except H, C, A and L):
0019 %       H : matrix (possibly sparse) of quadratic cost coefficients
0020 %       C : vector of linear cost coefficients
0021 %       A, L, U : define the optional linear constraints. Default
0022 %           values for the elements of L and U are -Inf and Inf,
0023 %           respectively.
0024 %       XMIN, XMAX : optional lower and upper bounds on the
0025 %           X variables, defaults are -Inf and Inf, respectively.
0026 %       X0 : optional starting value of optimization vector X
0027 %       OPT : optional options structure with the following fields,
0028 %           all of which are also optional (default values shown in
0029 %           parentheses)
0030 %           verbose (0) - controls level of progress output displayed
0031 %               0 = no progress output
0032 %               1 = some progress output
0033 %               2 = verbose progress output
0034 %           linprog_opt - options struct for LINPROG, value in verbose
0035 %                   overrides these options
0036 %           quadprog_opt - options struct for QUADPROG, value in verbose
0037 %                   overrides these options
0038 %       PROBLEM : The inputs can alternatively be supplied in a single
0039 %           PROBLEM struct with fields corresponding to the input arguments
0040 %           described above: H, c, A, l, u, xmin, xmax, x0, opt
0041 %
0042 %   Outputs:
0043 %       X : solution vector
0044 %       F : final objective function value
0045 %       EXITFLAG : QUADPROG/LINPROG exit flag
0046 %           (see QUADPROG and LINPROG documentation for details)
0047 %       OUTPUT : QUADPROG/LINPROG output struct
0048 %           (see QUADPROG and LINPROG documentation for details)
0049 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0050 %           multipliers on the constraints, with fields:
0051 %           mu_l - lower (left-hand) limit on linear constraints
0052 %           mu_u - upper (right-hand) limit on linear constraints
0053 %           lower - lower bound on optimization variables
0054 %           upper - upper bound on optimization variables
0055 %
0056 %   Note the calling syntax is almost identical to that of QUADPROG
0057 %   from MathWorks' Optimization Toolbox. The main difference is that
0058 %   the linear constraints are specified with A, L, U instead of
0059 %   A, B, Aeq, Beq.
0060 %
0061 %   Calling syntax options:
0062 %       [x, f, exitflag, output, lambda] = ...
0063 %           qps_ot(H, c, A, l, u, xmin, xmax, x0, opt)
0064 %
0065 %       x = qps_ot(H, c, A, l, u)
0066 %       x = qps_ot(H, c, A, l, u, xmin, xmax)
0067 %       x = qps_ot(H, c, A, l, u, xmin, xmax, x0)
0068 %       x = qps_ot(H, c, A, l, u, xmin, xmax, x0, opt)
0069 %       x = qps_ot(problem), where problem is a struct with fields:
0070 %                       H, c, A, l, u, xmin, xmax, x0, opt
0071 %                       all fields except 'c', 'A' and 'l' or 'u' are optional
0072 %       x = qps_ot(...)
0073 %       [x, f] = qps_ot(...)
0074 %       [x, f, exitflag] = qps_ot(...)
0075 %       [x, f, exitflag, output] = qps_ot(...)
0076 %       [x, f, exitflag, output, lambda] = qps_ot(...)
0077 %
0078 %
0079 %   Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm)
0080 %       H = [   1003.1  4.3     6.3     5.9;
0081 %               4.3     2.2     2.1     3.9;
0082 %               6.3     2.1     3.5     4.8;
0083 %               5.9     3.9     4.8     10  ];
0084 %       c = zeros(4,1);
0085 %       A = [   1       1       1       1;
0086 %               0.17    0.11    0.10    0.18    ];
0087 %       l = [1; 0.10];
0088 %       u = [1; Inf];
0089 %       xmin = zeros(4,1);
0090 %       x0 = [1; 0; 0; 1];
0091 %       opt = struct('verbose', 2);
0092 %       [x, f, s, out, lambda] = qps_ot(H, c, A, l, u, xmin, [], x0, opt);
0093 %
0094 %   See also QPS_MASTER, QUADPROG, LINPROG.
0095 
0096 %   MP-Opt-Model
0097 %   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
0098 %   by Ray Zimmerman, PSERC Cornell
0099 %
0100 %   This file is part of MP-Opt-Model.
0101 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0102 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0103 
0104 %% check for Optimization Toolbox
0105 % if ~have_feature('quadprog')
0106 %     error('qps_ot: requires the Optimization Toolbox');
0107 % end
0108 
0109 %%----- input argument handling  -----
0110 %% gather inputs
0111 if nargin == 1 && isstruct(H)       %% problem struct
0112     p = H;
0113     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0114     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0115     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0116     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0117     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0118     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0119     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0120     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0121     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0122 else                                %% individual args
0123     if nargin < 9
0124         opt = [];
0125         if nargin < 8
0126             x0 = [];
0127             if nargin < 7
0128                 xmax = [];
0129                 if nargin < 6
0130                     xmin = [];
0131                 end
0132             end
0133         end
0134     end
0135 end
0136 
0137 %% define nx, set default values for missing optional inputs
0138 if isempty(H) || ~any(any(H))
0139     if isempty(A) && isempty(xmin) && isempty(xmax)
0140         error('qps_ot: LP problem must include constraints or variable bounds');
0141     else
0142         if ~isempty(A)
0143             nx = size(A, 2);
0144         elseif ~isempty(xmin)
0145             nx = length(xmin);
0146         else    % if ~isempty(xmax)
0147             nx = length(xmax);
0148         end
0149     end
0150 else
0151     nx = size(H, 1);
0152 end
0153 if isempty(c)
0154     c = zeros(nx, 1);
0155 end
0156 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0157                                  (isempty(u) || all(u == Inf)))
0158     A = sparse(0,nx);           %% no limits => no linear constraints
0159 end
0160 nA = size(A, 1);                %% number of original linear constraints
0161 if isempty(u)                   %% By default, linear inequalities are ...
0162     u = Inf(nA, 1);             %% ... unbounded above and ...
0163 end
0164 if isempty(l)
0165     l = -Inf(nA, 1);            %% ... unbounded below.
0166 end
0167 if isempty(xmin)                %% By default, optimization variables are ...
0168     xmin = -Inf(nx, 1);         %% ... unbounded below and ...
0169 end
0170 if isempty(xmax)
0171     xmax = Inf(nx, 1);          %% ... unbounded above.
0172 end
0173 if isempty(x0)
0174     x0 = zeros(nx, 1);
0175 end
0176 if isempty(H) || ~any(any(H))
0177     isLP = 1;   %% it's an LP
0178 else
0179     isLP = 0;   %% nope, it's a QP
0180 end
0181 
0182 %% default options
0183 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0184     verbose = opt.verbose;
0185 else
0186     verbose = 0;
0187 end
0188 %% MATLAB or Octave
0189 matlab = have_feature('matlab');
0190 otver = have_feature('quadprog', 'vnum');
0191 
0192 %% split up linear constraints
0193 ieq = find( abs(u-l) <= eps );          %% equality
0194 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0195 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0196 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0197 Ae = A(ieq, :);
0198 be = u(ieq);
0199 Ai  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0200 bi  = [ u(ilt);    -l(igt);    u(ibx);    -l(ibx)];
0201 
0202 %% grab some dimensions
0203 nlt = length(ilt);      %% number of upper bounded linear inequalities
0204 ngt = length(igt);      %% number of lower bounded linear inequalities
0205 nbx = length(ibx);      %% number of doubly bounded linear inequalities
0206 
0207 %% set up options
0208 if verbose > 1
0209     vrb = 'iter';       %% seems to be same as 'final'
0210 elseif verbose == 1
0211     vrb = 'final';
0212 else
0213     vrb = 'off';
0214 end
0215 if have_feature('optimoptions')     %% Optimization Tbx 6.3 + (R2013a +)
0216     %% could use optimset for everything, except some options are not
0217     %% recognized by optimset, only optimoptions, such as
0218     %% ot_opt.Algorithm = 'dual-simplex'
0219     if isLP
0220         ot_opt = optimoptions('linprog');
0221         if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt)
0222             ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt);
0223         end
0224     else
0225         ot_opt = optimoptions('quadprog');
0226         if have_feature('quadprog_ls')
0227             ot_opt.Algorithm = 'interior-point-convex';
0228         else
0229             ot_opt.LargeScale = 'off';
0230         end
0231         if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt)
0232             ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt);
0233         end
0234     end
0235     ot_opt = optimoptions(ot_opt, 'Display', vrb);
0236 else                            %% need to use optimset()
0237     if isLP
0238         if matlab
0239             ot_opt = optimset('linprog');
0240         else
0241             ot_opt = optimset();
0242         end
0243         if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt)
0244             ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt);
0245         end
0246     else
0247         if matlab
0248             ot_opt = optimset('quadprog');
0249             if have_feature('quadprog_ls')
0250                 ot_opt = optimset(ot_opt, 'Algorithm', 'interior-point-convex');
0251             else
0252                 ot_opt = optimset(ot_opt, 'LargeScale', 'off');
0253             end
0254         else
0255             ot_opt = optimset();
0256         end
0257         if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt)
0258             ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt);
0259         end
0260     end
0261     ot_opt = optimset(ot_opt, 'Display', vrb);
0262 end
0263 
0264 %% call the solver
0265 if isLP
0266     if matlab
0267         [x, f, eflag, output, lam] = ...
0268             linprog(c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt);
0269     else
0270 % don't use linprog under Octave (using GLPK directly is recommended)
0271 %         [x, f] = linprog(c, Ai, bi, Ae, be, xmin, xmax);
0272 %         eflag = [];
0273 %         output = [];
0274 %         lam = [];
0275         [x, f, eflag, output, lam] = ...
0276             quadprog(sparse(nx,nx), c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt);
0277     end
0278 else
0279     [x, f, eflag, output, lam] = ...
0280         quadprog(H, c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt);
0281 end
0282 
0283 %% repackage lambdas
0284 if isempty(x)
0285     x = NaN(nx, 1);
0286 end
0287 if isempty(lam) || (isempty(lam.eqlin) && isempty(lam.ineqlin) && ...
0288                     isempty(lam.lower) && isempty(lam.upper))
0289     lambda = struct( ...
0290         'mu_l', NaN(nA, 1), ...
0291         'mu_u', NaN(nA, 1), ...
0292         'lower', NaN(nx, 1), ...
0293         'upper', NaN(nx, 1) ...
0294     );
0295 else
0296     kl = find(lam.eqlin < 0);   %% lower bound binding
0297     ku = find(lam.eqlin > 0);   %% upper bound binding
0298 
0299     mu_l = zeros(nA, 1);
0300 %     %% workaround for Octave optim 1.5.0 and earlier, which
0301 %     %% has opposite sign convention for equality multipliers
0302 %     if ~matlab && otver <= 1.005
0303 %         mu_l(ieq(ku)) = lam.eqlin(ku);
0304 %     else
0305         mu_l(ieq(kl)) = -lam.eqlin(kl);
0306 %     end
0307     mu_l(igt) = lam.ineqlin(nlt+(1:ngt));
0308     mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0309 
0310     mu_u = zeros(nA, 1);
0311 %     %% workaround for Octave optim 1.5.0 and earlier, which
0312 %     %% has opposite sign convention for equality multipliers
0313 %     if ~matlab && otver <= 1.005
0314 %         mu_u(ieq(kl)) = -lam.eqlin(kl);
0315 %     else
0316         mu_u(ieq(ku)) = lam.eqlin(ku);
0317 %     end
0318     mu_u(ilt) = lam.ineqlin(1:nlt);
0319     mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx));
0320 
0321     %% workaround for Octave optim 1.5.0 and earlier, which
0322     %% has opposite sign convention for equality multipliers
0323     % if ~matlab && otver <= 1.005
0324         %% there are also issues with variable bounds that are
0325         %% converted to equalities, and maybe other issues
0326     % end
0327 
0328     lambda = struct( ...
0329         'mu_l', mu_l, ...
0330         'mu_u', mu_u, ...
0331         'lower', lam.lower(1:nx), ...
0332         'upper', lam.upper(1:nx) ...
0333     );
0334 end

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005