Home > matpower7.1 > mp-opt-model > lib > qps_clp.m

qps_clp

PURPOSE ^

QPS_CLP Quadratic Program Solver based on CLP - COIN-OR Linear Programming.

SYNOPSIS ^

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

DESCRIPTION ^

QPS_CLP  Quadratic Program Solver based on CLP - COIN-OR Linear Programming.
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
       QPS_CLP(H, C, A, L, U, XMIN, XMAX, X0, OPT)
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_CLP(PROBLEM)
   A wrapper function providing a standardized interface for using
   CLP 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 (NOT USED)
       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
           clp_opt - options struct for CLP, 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 : exit flag, 1 - optimal, -1 - infeasible, -2 - unbounded
                   -3 - max iterations/time exceeded
       OUTPUT : struct with fields
           exitflag - raw CLP exit flag: 0 - optimal, 1 - infeasible,
               2 - unbounded, 3 - max iterations/time exceeded
           status - string with explanation of exitflag
           (iter - depending on build of solver this may contain
                   the number of iterations)
       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 CLP. 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_clp(H, c, A, l, u, xmin, xmax, x0, opt)

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


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

   See also QPS_MASTER, CLP.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x, f, eflag, output, lambda] = qps_clp(H, c, A, l, u, xmin, xmax, x0, opt)
0002 %QPS_CLP  Quadratic Program Solver based on CLP - COIN-OR Linear Programming.
0003 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
0004 %       QPS_CLP(H, C, A, L, U, XMIN, XMAX, X0, OPT)
0005 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_CLP(PROBLEM)
0006 %   A wrapper function providing a standardized interface for using
0007 %   CLP to solve the 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 (NOT USED)
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 %           clp_opt - options struct for CLP, value in verbose
0034 %                   overrides these options
0035 %       PROBLEM : The inputs can alternatively be supplied in a single
0036 %           PROBLEM struct with fields corresponding to the input arguments
0037 %           described above: H, c, A, l, u, xmin, xmax, x0, opt
0038 %
0039 %   Outputs:
0040 %       X : solution vector
0041 %       F : final objective function value
0042 %       EXITFLAG : exit flag, 1 - optimal, -1 - infeasible, -2 - unbounded
0043 %                   -3 - max iterations/time exceeded
0044 %       OUTPUT : struct with fields
0045 %           exitflag - raw CLP exit flag: 0 - optimal, 1 - infeasible,
0046 %               2 - unbounded, 3 - max iterations/time exceeded
0047 %           status - string with explanation of exitflag
0048 %           (iter - depending on build of solver this may contain
0049 %                   the number of iterations)
0050 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0051 %           multipliers on the constraints, with fields:
0052 %           mu_l - lower (left-hand) limit on linear constraints
0053 %           mu_u - upper (right-hand) limit on linear constraints
0054 %           lower - lower bound on optimization variables
0055 %           upper - upper bound on optimization variables
0056 %
0057 %   Note the calling syntax is almost identical to that of CLP. The main
0058 %   difference is that the linear constraints are specified with A, L, U
0059 %   instead of A, B, Aeq, Beq.
0060 %
0061 %   Calling syntax options:
0062 %       [x, f, exitflag, output, lambda] = ...
0063 %           qps_clp(H, c, A, l, u, xmin, xmax, x0, opt)
0064 %
0065 %       x = qps_clp(H, c, A, l, u)
0066 %       x = qps_clp(H, c, A, l, u, xmin, xmax)
0067 %       x = qps_clp(H, c, A, l, u, xmin, xmax, x0)
0068 %       x = qps_clp(H, c, A, l, u, xmin, xmax, x0, opt)
0069 %       x = qps_clp(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_clp(...)
0073 %       [x, f] = qps_clp(...)
0074 %       [x, f, exitflag] = qps_clp(...)
0075 %       [x, f, exitflag, output] = qps_clp(...)
0076 %       [x, f, exitflag, output, lambda] = qps_clp(...)
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_clp(H, c, A, l, u, xmin, [], x0, opt);
0093 %
0094 %   See also QPS_MASTER, CLP.
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_clp: requires the MEX interface to CLP');
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_clp: 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 ~issparse(A)
0177     A = sparse(A);
0178 end
0179 if ~issparse(H)
0180     H = sparse(H);
0181 end
0182 
0183 
0184 %% default options
0185 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0186     verbose = opt.verbose;
0187 else
0188     verbose = 0;
0189 end
0190 
0191 %% set up options struct for CLP
0192 if ~isempty(opt) && isfield(opt, 'clp_opt') && ~isempty(opt.clp_opt)
0193     clp_opt = clp_options(opt.clp_opt);
0194 else
0195     clp_opt = clp_options;
0196 end
0197 
0198 if have_feature('opti_clp')     %% use OPTI Toolbox verision's MEX interface
0199     clp_opt.display = verbose;
0200 
0201     [x, f, exitflag, iter, lam] = clp(tril(H), c, A, l, u, xmin, xmax, clp_opt);
0202     
0203     output.iter = iter;
0204 
0205     %% repackage lambdas
0206 %     if isempty(x)
0207 %         x = NaN(nx, 1);
0208 %         f = NaN;
0209 %     end
0210 %     if isempty(lam)
0211 %         lambda = struct( ...
0212 %             'mu_l', zeros(nA, 1), ...
0213 %             'mu_u', zeros(nA, 1), ...
0214 %             'lower', zeros(nx, 1), ...
0215 %             'upper', zeros(nx, 1) ...
0216 %         );
0217 %     else
0218         mu_l = lam.dual_row;
0219         mu_u = -lam.dual_row;
0220         lower = lam.dual_col;
0221         upper = -lam.dual_col;
0222         mu_l(mu_l < 0) = 0;
0223         mu_u(mu_u < 0) = 0;
0224         lower(lower < 0) = 0;
0225         upper(upper < 0) = 0;
0226         
0227         lambda = struct( ...
0228             'mu_l', mu_l, ...
0229             'mu_u', mu_u, ...
0230             'lower', lower, ...
0231             'upper', upper ...
0232         );
0233 %     end
0234 else
0235     clp_opt.verbose = verbose;
0236 
0237     %% split up linear constraints
0238     ieq = find( abs(u-l) <= eps );          %% equality
0239     igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0240     ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0241     ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0242     Ae = A(ieq, :);
0243     be = u(ieq);
0244     Ai  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0245     bi  = [ u(ilt);    -l(igt);    u(ibx);    -l(ibx)];
0246 
0247     %% grab some dimensions
0248     nlt = length(ilt);      %% number of upper bounded linear inequalities
0249     ngt = length(igt);      %% number of lower bounded linear inequalities
0250     nbx = length(ibx);      %% number of doubly bounded linear inequalities
0251 
0252     %% call the solver
0253     [x, z, exitflag] = ...
0254         clp(H, c, Ai, bi, Ae, be, xmin, xmax, clp_opt);
0255 
0256     %% repackage lambdas
0257     if isempty(x)
0258         x = NaN(nx, 1);
0259         f = NaN;
0260     else
0261         if isempty(H) || ~any(any(H))
0262             f = c'*x;
0263         else
0264             f = 0.5 * x'*H*x + c'*x;
0265         end
0266     end
0267     if isempty(z)
0268         lambda = struct( ...
0269             'mu_l', zeros(nA, 1), ...
0270             'mu_u', zeros(nA, 1), ...
0271             'lower', zeros(nx, 1), ...
0272             'upper', zeros(nx, 1) ...
0273         );
0274     else
0275         neq = length(be);
0276         nie = length(bi);
0277         lam.eqlin = z(1:neq);
0278         lam.ineqlin = z(neq+(1:nie));
0279     %%-----  MEXCLP DOES NOT RETURN MULTIPLIERS ON VARIABLE BOUNDS :-/  -----
0280         lam.lower = NaN(nx, 1);
0281         lam.upper = NaN(nx, 1);
0282         kl = find(lam.eqlin > 0);   %% lower bound binding
0283         ku = find(lam.eqlin < 0);   %% upper bound binding
0284 
0285         mu_l = zeros(nA, 1);
0286         mu_l(ieq(kl)) = lam.eqlin(kl);
0287         mu_l(igt) = -lam.ineqlin(nlt+(1:ngt));
0288         mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0289 
0290         mu_u = zeros(nA, 1);
0291         mu_u(ieq(ku)) = -lam.eqlin(ku);
0292         mu_u(ilt) = -lam.ineqlin(1:nlt);
0293         mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx));
0294 
0295         lambda = struct( ...
0296             'mu_l', mu_l, ...
0297             'mu_u', mu_u, ...
0298             'lower', lam.lower, ...
0299             'upper', lam.upper ...
0300         );
0301     end
0302 end
0303 
0304 %% set eflag
0305 eflag = -exitflag;
0306 if eflag == 0       %% success
0307     eflag = 1;
0308 end
0309 
0310 %% set status
0311 output.exitflag = exitflag;
0312 switch exitflag
0313     case 0
0314         output.status = 'optimal';
0315     case 1
0316         output.status = 'primal infeasible';
0317     case 2
0318         output.status = 'dual infeasible';
0319     case 3
0320         output.status = 'max iterations or time exceeded';
0321     otherwise
0322         output.status = 'unknown exit code';
0323 end

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