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

qps_osqp

PURPOSE ^

QPS_OSQP Quadratic Program Solver based on OSQP.

SYNOPSIS ^

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

DESCRIPTION ^

QPS_OSQP  Quadratic Program Solver based on OSQP.
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
       QPS_OSQP(H, C, A, L, U, XMIN, XMAX, X0, OPT)
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_OSQP(PROBLEM)
   A wrapper function providing a standardized interface for using
   OSQP 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
               3 = even more verbose progress output
           osqp_opt - options struct for OSQP (see
               https://osqp.org/docs/interfaces/solver_settings.html),
               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 : OSQP exit flag
           1 = converged
           0 or negative values OSQP status value
           (see OSQP documentation for details)
       OUTPUT : OSQP results struct
           (see OSQP 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_osqp(H, c, A, l, u, xmin, xmax, x0, opt)

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


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

   See also QPS_MASTER, OSQP.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x, f, eflag, output, lambda] = qps_osqp(H, c, A, l, u, xmin, xmax, x0, opt)
0002 %QPS_OSQP  Quadratic Program Solver based on OSQP.
0003 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
0004 %       QPS_OSQP(H, C, A, L, U, XMIN, XMAX, X0, OPT)
0005 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_OSQP(PROBLEM)
0006 %   A wrapper function providing a standardized interface for using
0007 %   OSQP 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
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 %               3 = even more verbose progress output
0034 %           osqp_opt - options struct for OSQP (see
0035 %               https://osqp.org/docs/interfaces/solver_settings.html),
0036 %               value in 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 : OSQP exit flag
0045 %           1 = converged
0046 %           0 or negative values OSQP status value
0047 %           (see OSQP documentation for details)
0048 %       OUTPUT : OSQP results struct
0049 %           (see OSQP documentation for details)
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 QUADPROG
0058 %   from MathWorks' Optimization Toolbox. The main difference is that
0059 %   the linear constraints are specified with A, L, U instead of
0060 %   A, B, Aeq, Beq.
0061 %
0062 %   Calling syntax options:
0063 %       [x, f, exitflag, output, lambda] = ...
0064 %           qps_osqp(H, c, A, l, u, xmin, xmax, x0, opt)
0065 %
0066 %       x = qps_osqp(H, c, A, l, u)
0067 %       x = qps_osqp(H, c, A, l, u, xmin, xmax)
0068 %       x = qps_osqp(H, c, A, l, u, xmin, xmax, x0)
0069 %       x = qps_osqp(H, c, A, l, u, xmin, xmax, x0, opt)
0070 %       x = qps_osqp(problem), where problem is a struct with fields:
0071 %                       H, c, A, l, u, xmin, xmax, x0, opt
0072 %                       all fields except 'c', 'A' and 'l' or 'u' are optional
0073 %       x = qps_osqp(...)
0074 %       [x, f] = qps_osqp(...)
0075 %       [x, f, exitflag] = qps_osqp(...)
0076 %       [x, f, exitflag, output] = qps_osqp(...)
0077 %       [x, f, exitflag, output, lambda] = qps_osqp(...)
0078 %
0079 %
0080 %   Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm)
0081 %       H = [   1003.1  4.3     6.3     5.9;
0082 %               4.3     2.2     2.1     3.9;
0083 %               6.3     2.1     3.5     4.8;
0084 %               5.9     3.9     4.8     10  ];
0085 %       c = zeros(4,1);
0086 %       A = [   1       1       1       1;
0087 %               0.17    0.11    0.10    0.18    ];
0088 %       l = [1; 0.10];
0089 %       u = [1; Inf];
0090 %       xmin = zeros(4,1);
0091 %       x0 = [1; 0; 0; 1];
0092 %       opt = struct('verbose', 2);
0093 %       [x, f, s, out, lambda] = qps_osqp(H, c, A, l, u, xmin, [], x0, opt);
0094 %
0095 %   See also QPS_MASTER, OSQP.
0096 
0097 %   MP-Opt-Model
0098 %   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
0099 %   by Ray Zimmerman, PSERC Cornell
0100 %
0101 %   This file is part of MP-Opt-Model.
0102 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0103 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0104 
0105 %% check for OSQP
0106 % if ~have_feature('osqp')
0107 %     error('qps_osqp: requires OSQP (https://osqp.org)');
0108 % end
0109 
0110 %%----- input argument handling  -----
0111 %% gather inputs
0112 if nargin == 1 && isstruct(H)       %% problem struct
0113     p = H;
0114     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0115     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0116     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0117     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0118     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0119     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0120     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0121     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0122     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0123 else                                %% individual args
0124     if nargin < 9
0125         opt = [];
0126         if nargin < 8
0127             x0 = [];
0128             if nargin < 7
0129                 xmax = [];
0130                 if nargin < 6
0131                     xmin = [];
0132                 end
0133             end
0134         end
0135     end
0136 end
0137 
0138 %% define nx, set default values for missing optional inputs
0139 if isempty(H) || ~any(any(H))
0140     lpqp = 'LP';
0141     if isempty(A) && isempty(xmin) && isempty(xmax)
0142         error('qps_osqp: LP problem must include constraints or variable bounds');
0143     else
0144         if ~isempty(A)
0145             nx = size(A, 2);
0146         elseif ~isempty(xmin)
0147             nx = length(xmin);
0148         else    % if ~isempty(xmax)
0149             nx = length(xmax);
0150         end
0151     end
0152 else
0153     lpqp = 'QP';
0154     nx = size(H, 1);
0155 end
0156 if isempty(c)
0157     c = zeros(nx, 1);
0158 end
0159 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0160                                  (isempty(u) || all(u == Inf)))
0161     A = sparse(0,nx);           %% no limits => no linear constraints
0162 end
0163 nA = size(A, 1);                %% number of original linear constraints
0164 if isempty(u)                   %% By default, linear inequalities are ...
0165     u = Inf(nA, 1);             %% ... unbounded above and ...
0166 end
0167 if isempty(l)
0168     l = -Inf(nA, 1);            %% ... unbounded below.
0169 end
0170 
0171 %% default options
0172 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0173     verbose = opt.verbose;
0174 else
0175     verbose = 0;
0176 end
0177 
0178 %% set up options struct for OSQP
0179 osqp_opt = struct();
0180 if ~isempty(opt) && isfield(opt, 'osqp_opt') && ~isempty(opt.osqp_opt)
0181     osqp_opt = nested_struct_copy(osqp_opt, opt.osqp_opt);
0182 end
0183 if verbose > 1
0184     osqp_opt.verbose = 1;
0185 %     if verbose > 2
0186 %         osqp_opt.
0187 %     else
0188 %         osqp_opt.
0189 %     end
0190 else
0191     osqp_opt.verbose = 0;
0192 end
0193 
0194 %% add variable bounds to linear constraints
0195 if isempty(xmin)
0196     if isempty(xmax)
0197         k = [];
0198     else
0199         k = find(xmax < Inf);
0200     end
0201 else
0202     if isempty(xmax)
0203         k = find(xmin > -Inf);
0204     else
0205         k = find(xmin > -Inf | xmax < Inf);
0206     end
0207 end
0208 nv = length(k);
0209 if nv
0210     Av = sparse(1:nv, k, 1, nv, nx);
0211     if isempty(xmin)
0212         lv = -Inf(nv,1);
0213     else
0214         lv = xmin(k);
0215     end
0216     if isempty(xmax)
0217         uv = Inf(nv,1);
0218     else
0219         uv = xmax(k);
0220     end
0221     AA = [A; Av];
0222     ll = [l; lv];
0223     uu = [u; uv];
0224 else
0225     AA = A;
0226     ll = l;
0227     uu = u;
0228 end
0229 
0230 %% set up model
0231 o = osqp();
0232 o.setup(H, c, AA, ll, uu, osqp_opt);
0233 if ~isempty(x0)
0234     o.warm_start('x', x0);
0235 end
0236 
0237 %% solve model
0238 if verbose
0239     vn = osqpver;
0240     fprintf('OSQP Version %s -- %s solver\n', vn, lpqp);
0241 end
0242 res = o.solve();
0243 if verbose
0244     fprintf('OSQP solution status: %s\n', res.info.status);
0245 end
0246 
0247 %% extract results
0248 x = res.x;
0249 f = res.info.obj_val;
0250 eflag = res.info.status_val;
0251 if eflag > 1
0252     eflag = -10 * eflag;
0253 end
0254 output = res;
0255 
0256 %% separate duals into binding lower & upper bounds
0257 mu_l = -res.y;
0258 mu_u =  res.y;
0259 mu_l(mu_l < 0) = 0;
0260 mu_u(mu_u < 0) = 0;
0261 
0262 %% variable bounds
0263 lam_lower = zeros(nx, 1);
0264 lam_lower(k) = mu_l(nA+1:end);
0265 lam_upper = zeros(nx, 1);
0266 lam_upper(k) = mu_u(nA+1:end);
0267 
0268 %% package lambdas
0269 lambda = struct( ...
0270     'mu_l', mu_l(1:nA), ...
0271     'mu_u', mu_u(1:nA), ...
0272     'lower', lam_lower, ...
0273     'upper', lam_upper ...
0274 );

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