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

qps_ipopt

PURPOSE ^

QPS_IPOPT Quadratic Program Solver based on IPOPT.

SYNOPSIS ^

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

DESCRIPTION ^

QPS_IPOPT  Quadratic Program Solver based on IPOPT.
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
       QPS_IPOPT(H, C, A, L, U, XMIN, XMAX, X0, OPT)
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_IPOPT(PROBLEM)
   A wrapper function providing a standardized interface for using
   IPOPT 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
           ipopt_opt - options struct for IPOPT, 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 = converged
           0 = failed to converge
       OUTPUT : output struct with the following fields:
           status - see IPOPT documentation for INFO.status
             https://coin-or.github.io/Ipopt/IpReturnCodes__inc_8h_source.html
           iterations - number of iterations performed (INFO.iter)
           cpu - see IPOPT documentation for INFO.cpu
           eval - see IPOPT documentation for INFO.eval
       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_ipopt(H, c, A, l, u, xmin, xmax, x0, opt)

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

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

   See also QPS_MASTER, IPOPT, IPOPT_OPTIONS.
   https://github.com/coin-or/Ipopt.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x, f, eflag, output, lambda] = qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt)
0002 %QPS_IPOPT  Quadratic Program Solver based on IPOPT.
0003 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
0004 %       QPS_IPOPT(H, C, A, L, U, XMIN, XMAX, X0, OPT)
0005 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_IPOPT(PROBLEM)
0006 %   A wrapper function providing a standardized interface for using
0007 %   IPOPT 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 %           ipopt_opt - options struct for IPOPT, 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
0043 %           1 = converged
0044 %           0 = failed to converge
0045 %       OUTPUT : output struct with the following fields:
0046 %           status - see IPOPT documentation for INFO.status
0047 %             https://coin-or.github.io/Ipopt/IpReturnCodes__inc_8h_source.html
0048 %           iterations - number of iterations performed (INFO.iter)
0049 %           cpu - see IPOPT documentation for INFO.cpu
0050 %           eval - see IPOPT documentation for INFO.eval
0051 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0052 %           multipliers on the constraints, with fields:
0053 %           mu_l - lower (left-hand) limit on linear constraints
0054 %           mu_u - upper (right-hand) limit on linear constraints
0055 %           lower - lower bound on optimization variables
0056 %           upper - upper bound on optimization variables
0057 %
0058 %   Note the calling syntax is almost identical to that of QUADPROG
0059 %   from MathWorks' Optimization Toolbox. The main difference is that
0060 %   the linear constraints are specified with A, L, U instead of
0061 %   A, B, Aeq, Beq.
0062 %
0063 %   Calling syntax options:
0064 %       [x, f, exitflag, output, lambda] = ...
0065 %           qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt)
0066 %
0067 %       x = qps_ipopt(H, c, A, l, u)
0068 %       x = qps_ipopt(H, c, A, l, u, xmin, xmax)
0069 %       x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0)
0070 %       x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt)
0071 %       x = qps_ipopt(problem), where problem is a struct with fields:
0072 %                       H, c, A, l, u, xmin, xmax, x0, opt
0073 %                       all fields except 'c', 'A' and 'l' or 'u' are optional
0074 %       x = qps_ipopt(...)
0075 %       [x, f] = qps_ipopt(...)
0076 %       [x, f, exitflag] = qps_ipopt(...)
0077 %       [x, f, exitflag, output] = qps_ipopt(...)
0078 %       [x, f, exitflag, output, lambda] = qps_ipopt(...)
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_ipopt(H, c, A, l, u, xmin, [], x0, opt);
0094 %
0095 %   See also QPS_MASTER, IPOPT, IPOPT_OPTIONS.
0096 %   https://github.com/coin-or/Ipopt.
0097 
0098 %   MP-Opt-Model
0099 %   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
0100 %   by Ray Zimmerman, PSERC Cornell
0101 %
0102 %   This file is part of MP-Opt-Model.
0103 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0104 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0105 
0106 %% check for IPOPT
0107 % if ~have_feature('ipopt')
0108 %     error('qps_ipopt: requires IPOPT (https://github.com/coin-or/Ipopt)');
0109 % end
0110 
0111 %%----- input argument handling  -----
0112 %% gather inputs
0113 if nargin == 1 && isstruct(H)       %% problem struct
0114     p = H;
0115     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0116     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0117     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0118     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0119     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0120     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0121     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0122     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0123     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0124 else                                %% individual args
0125     if nargin < 9
0126         opt = [];
0127         if nargin < 8
0128             x0 = [];
0129             if nargin < 7
0130                 xmax = [];
0131                 if nargin < 6
0132                     xmin = [];
0133                 end
0134             end
0135         end
0136     end
0137 end
0138 
0139 %% define nx, set default values for missing optional inputs
0140 if isempty(H) || ~any(any(H))
0141     if isempty(A) && isempty(xmin) && isempty(xmax)
0142         error('qps_ipopt: 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     H = sparse(nx,nx);
0153 else
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 nA
0165     if isempty(u)               %% By default, linear inequalities are ...
0166         u = Inf(nA, 1);             %% ... unbounded above and ...
0167     end
0168     if isempty(l)
0169         l = -Inf(nA, 1);            %% ... unbounded below.
0170     end
0171 end
0172 if isempty(x0)
0173     x0 = zeros(nx, 1);
0174 end
0175 
0176 %% default options
0177 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0178     verbose = opt.verbose;
0179 else
0180     verbose = 0;
0181 end
0182 
0183 %% make sure args are sparse/full as expected by IPOPT
0184 if ~isempty(H)
0185     if ~issparse(H)
0186         H = sparse(H);
0187     end
0188 end
0189 if ~issparse(A)
0190     A = sparse(A);
0191 end
0192 
0193 
0194 %%-----  run optimization  -----
0195 %% set options struct for IPOPT
0196 if ~isempty(opt) && isfield(opt, 'ipopt_opt') && ~isempty(opt.ipopt_opt)
0197     options.ipopt = ipopt_options(opt.ipopt_opt);
0198 else
0199     options.ipopt = ipopt_options;
0200 end
0201 options.ipopt.jac_c_constant    = 'yes';
0202 options.ipopt.jac_d_constant    = 'yes';
0203 options.ipopt.hessian_constant  = 'yes';
0204 options.ipopt.least_square_init_primal  = 'yes';
0205 options.ipopt.least_square_init_duals   = 'yes';
0206 % options.ipopt.mehrotra_algorithm        = 'yes';        %% default 'no'
0207 if verbose
0208     options.ipopt.print_level = min(12, verbose*2+1);
0209 else
0210     options.ipopt.print_level = 0;
0211 end
0212 
0213 %% define variable and constraint bounds, if given
0214 if nA
0215     options.cu = u;
0216     options.cl = l;
0217 end
0218 if ~isempty(xmin)
0219     options.lb = xmin;
0220 end
0221 if ~isempty(xmax)
0222     options.ub = xmax;
0223 end
0224 
0225 %% assign function handles
0226 funcs.objective         = @(x) 0.5 * x' * H * x + c' * x;
0227 funcs.gradient          = @(x) H * x + c;
0228 funcs.constraints       = @(x) A * x;
0229 funcs.jacobian          = @(x) A;
0230 funcs.jacobianstructure = @() A;
0231 funcs.hessian           = @(x, sigma, lambda) tril(H);
0232 funcs.hessianstructure  = @() tril(H);
0233 
0234 %% run the optimization
0235 [x, info] = ipopt(x0,funcs,options);
0236 
0237 if info.status == 0 || info.status == 1
0238     eflag = 1;
0239 else
0240     eflag = 0;
0241 end
0242 output = struct('status', info.status);
0243 if isfield(info, 'iter')
0244     output.iterations = info.iter;
0245 else
0246     output.iterations = [];
0247 end
0248 if isfield(info, 'cpu')
0249     output.cpu = info.cpu;
0250 end
0251 if isfield(info, 'eval')
0252     output.eval = info.eval;
0253 end
0254 f = funcs.objective(x);
0255 
0256 %% check for empty results (in case optimization failed)
0257 if isempty(info.lambda)
0258     lam = NaN(nA, 1);
0259 else
0260     lam = info.lambda;
0261 end
0262 if isempty(info.zl) && ~isempty(xmin)
0263     zl = NaN(nx, 1);
0264 else
0265     zl = info.zl;
0266 end
0267 if isempty(info.zu) && ~isempty(xmax)
0268     zu = NaN(nx, 1);
0269 else
0270     zu = info.zu;
0271 end
0272 
0273 %% repackage lambdas
0274 kl = find(lam < 0);                     %% lower bound binding
0275 ku = find(lam > 0);                     %% upper bound binding
0276 mu_l = zeros(nA, 1);
0277 mu_l(kl) = -lam(kl);
0278 mu_u = zeros(nA, 1);
0279 mu_u(ku) = lam(ku);
0280 
0281 lambda = struct( ...
0282     'mu_l', mu_l, ...
0283     'mu_u', mu_u, ...
0284     'lower', zl, ...
0285     'upper', zu    );

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