Home > matpower4.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
           max_it (0) - maximum number of iterations allowed
               0 = use algorithm default
           ot_opt - options struct for QUADPROG/LINPROG, values in
               verbose and max_it override 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, lam] = 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 %           max_it (0) - maximum number of iterations allowed
0034 %               0 = use algorithm default
0035 %           ot_opt - options struct for QUADPROG/LINPROG, values in
0036 %               verbose and max_it override 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, lam] = qps_ot(H, c, A, l, u, xmin, [], x0, opt);
0092 %
0093 %   See also QUADPROG, LINPROG.
0094 
0095 %   MATPOWER
0096 %   $Id: qps_ot.m,v 1.12 2011/01/18 15:23:36 cvs Exp $
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 
0197 %% default options
0198 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0199     verbose = opt.verbose;
0200 else
0201     verbose = 0;
0202 end
0203 if ~isempty(opt) && isfield(opt, 'max_it') && ~isempty(opt.max_it)
0204     max_it = opt.max_it;
0205 else
0206     max_it = 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 ~isempty(opt) && isfield(opt, 'ot_opt') && ~isempty(opt.ot_opt)
0226     ot_opt = opt.ot_opt;
0227 else
0228     if isempty(H) || ~any(any(H))
0229         ot_opt = optimset('linprog');
0230     else
0231         ot_opt = optimset('quadprog');
0232         if have_fcn('quadprog_ls')
0233             ot_opt = optimset(ot_opt, 'Algorithm', 'interior-point-convex');
0234         else
0235             ot_opt = optimset(ot_opt, 'LargeScale', 'off');
0236         end
0237     end
0238 end
0239 if max_it
0240     ot_opt = optimset(ot_opt, 'MaxIter', max_it);
0241 end
0242 if verbose > 1
0243     ot_opt = optimset(ot_opt, 'Display', 'iter');   %% seems to be same as 'final'
0244 elseif verbose == 1
0245     ot_opt = optimset(ot_opt, 'Display', 'final');
0246 else
0247     ot_opt = optimset(ot_opt, 'Display', 'off');
0248 end
0249 
0250 %% call the solver
0251 if isempty(H) || ~any(any(H))
0252     [x, f, eflag, output, lam] = ...
0253         linprog(c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt);
0254 else
0255     [x, f, eflag, output, lam] = ...
0256         quadprog(H, c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt);
0257 end
0258 
0259 %% repackage lambdas
0260 kl = find(lam.eqlin < 0);   %% lower bound binding
0261 ku = find(lam.eqlin > 0);   %% upper bound binding
0262 
0263 mu_l = zeros(nA, 1);
0264 mu_l(ieq(kl)) = -lam.eqlin(kl);
0265 mu_l(igt) = lam.ineqlin(nlt+(1:ngt));
0266 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0267 
0268 mu_u = zeros(nA, 1);
0269 mu_u(ieq(ku)) = lam.eqlin(ku);
0270 mu_u(ilt) = lam.ineqlin(1:nlt);
0271 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx));
0272 
0273 lambda = struct( ...
0274     'mu_l', mu_l, ...
0275     'mu_u', mu_u, ...
0276     'lower', lam.lower(1:nx), ...
0277     'upper', lam.upper(1:nx) ...
0278 );

Generated on Mon 26-Jan-2015 14:56:45 by m2html © 2005