Home > matpower5.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 %   $Id: qps_ot.m 2385 2014-10-11 15:33:45Z ray $
0097 %   by Ray Zimmerman, PSERC Cornell
0098 %   Copyright (c) 2010-2011 by Power System Engineering Research Center (PSERC)
0099 %
0100 %   This file is part of MATPOWER.
0101 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0102 %
0103 %   MATPOWER is free software: you can redistribute it and/or modify
0104 %   it under the terms of the GNU General Public License as published
0105 %   by the Free Software Foundation, either version 3 of the License,
0106 %   or (at your option) any later version.
0107 %
0108 %   MATPOWER is distributed in the hope that it will be useful,
0109 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0110 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0111 %   GNU General Public License for more details.
0112 %
0113 %   You should have received a copy of the GNU General Public License
0114 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0115 %
0116 %   Additional permission under GNU GPL version 3 section 7
0117 %
0118 %   If you modify MATPOWER, or any covered work, to interface with
0119 %   other modules (such as MATLAB code and MEX-files) available in a
0120 %   MATLAB(R) or comparable environment containing parts covered
0121 %   under other licensing terms, the licensors of MATPOWER grant
0122 %   you additional permission to convey the resulting work.
0123 
0124 %% check for Optimization Toolbox
0125 % if ~have_fcn('quadprog')
0126 %     error('qps_ot: requires the Optimization Toolbox');
0127 % end
0128 
0129 %%----- input argument handling  -----
0130 %% gather inputs
0131 if nargin == 1 && isstruct(H)       %% problem struct
0132     p = H;
0133     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0134     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0135     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0136     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0137     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0138     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0139     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0140     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0141     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0142 else                                %% individual args
0143     if nargin < 9
0144         opt = [];
0145         if nargin < 8
0146             x0 = [];
0147             if nargin < 7
0148                 xmax = [];
0149                 if nargin < 6
0150                     xmin = [];
0151                 end
0152             end
0153         end
0154     end
0155 end
0156 
0157 %% define nx, set default values for missing optional inputs
0158 if isempty(H) || ~any(any(H))
0159     if isempty(A) && isempty(xmin) && isempty(xmax)
0160         error('qps_ot: LP problem must include constraints or variable bounds');
0161     else
0162         if ~isempty(A)
0163             nx = size(A, 2);
0164         elseif ~isempty(xmin)
0165             nx = length(xmin);
0166         else    % if ~isempty(xmax)
0167             nx = length(xmax);
0168         end
0169     end
0170 else
0171     nx = size(H, 1);
0172 end
0173 if isempty(c)
0174     c = zeros(nx, 1);
0175 end
0176 if  ~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0177                    (isempty(u) || all(u == Inf))
0178     A = sparse(0,nx);           %% no limits => no linear constraints
0179 end
0180 nA = size(A, 1);                %% number of original linear constraints
0181 if isempty(u)                   %% By default, linear inequalities are ...
0182     u = Inf * ones(nA, 1);      %% ... unbounded above and ...
0183 end
0184 if isempty(l)
0185     l = -Inf * ones(nA, 1);     %% ... unbounded below.
0186 end
0187 if isempty(xmin)                %% By default, optimization variables are ...
0188     xmin = -Inf * ones(nx, 1);  %% ... unbounded below and ...
0189 end
0190 if isempty(xmax)
0191     xmax = Inf * ones(nx, 1);   %% ... unbounded above.
0192 end
0193 if isempty(x0)
0194     x0 = zeros(nx, 1);
0195 end
0196 if isempty(H) || ~any(any(H))
0197     isLP = 1;   %% it's an LP
0198 else
0199     isLP = 0;   %% nope, it's a QP
0200 end
0201 
0202 %% default options
0203 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0204     verbose = opt.verbose;
0205 else
0206     verbose = 0;
0207 end
0208 
0209 %% split up linear constraints
0210 ieq = find( abs(u-l) <= eps );          %% equality
0211 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0212 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0213 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0214 Ae = A(ieq, :);
0215 be = u(ieq);
0216 Ai  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0217 bi  = [ u(ilt);    -l(igt);    u(ibx);    -l(ibx)];
0218 
0219 %% grab some dimensions
0220 nlt = length(ilt);      %% number of upper bounded linear inequalities
0221 ngt = length(igt);      %% number of lower bounded linear inequalities
0222 nbx = length(ibx);      %% number of doubly bounded linear inequalities
0223 
0224 %% set up options
0225 if verbose > 1
0226     vrb = 'iter';       %% seems to be same as 'final'
0227 elseif verbose == 1
0228     vrb = 'final';
0229 else
0230     vrb = 'off';
0231 end
0232 if have_fcn('optimoptions')     %% Optimization Tbx 6.3 + (R2013a +)
0233     %% could use optimset for everything, except some options are not
0234     %% recognized by optimset, only optimoptions, such as
0235     %% ot_opt.Algorithm = 'dual-simplex'
0236     if isLP
0237         ot_opt = optimoptions('linprog');
0238         if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt)
0239             ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt);
0240         end
0241     else
0242         ot_opt = optimoptions('quadprog');
0243         if have_fcn('quadprog_ls')
0244             ot_opt.Algorithm = 'interior-point-convex';
0245         else
0246             ot_opt.LargeScale = 'off';
0247         end
0248         if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt)
0249             ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt);
0250         end
0251     end
0252     ot_opt = optimoptions(ot_opt, 'Display', vrb);
0253 else                            %% need to use optimset()
0254     if isLP
0255         ot_opt = optimset('linprog');
0256         if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt)
0257             ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt);
0258         end
0259     else
0260         ot_opt = optimset('quadprog');
0261         if have_fcn('quadprog_ls')
0262             ot_opt = optimset(ot_opt, 'Algorithm', 'interior-point-convex');
0263         else
0264             ot_opt = optimset(ot_opt, 'LargeScale', 'off');
0265         end
0266         if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt)
0267             ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt);
0268         end
0269     end
0270     ot_opt = optimset(ot_opt, 'Display', vrb);
0271 end
0272 
0273 %% call the solver
0274 if isLP
0275     [x, f, eflag, output, lam] = ...
0276         linprog(c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt);
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     mu_l(ieq(kl)) = -lam.eqlin(kl);
0300     mu_l(igt) = lam.ineqlin(nlt+(1:ngt));
0301     mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0302 
0303     mu_u = zeros(nA, 1);
0304     mu_u(ieq(ku)) = lam.eqlin(ku);
0305     mu_u(ilt) = lam.ineqlin(1:nlt);
0306     mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx));
0307 
0308     lambda = struct( ...
0309         'mu_l', mu_l, ...
0310         'mu_u', mu_u, ...
0311         'lower', lam.lower(1:nx), ...
0312         'upper', lam.upper(1:nx) ...
0313     );
0314 end

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