Home > matpower5.1 > 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)
   Uses 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 = first order optimality conditions satisfied
           0 = maximum number of iterations reached
           -1 = numerically failed
       OUTPUT : output struct with the following fields:
           iterations - number of iterations performed
           hist - struct array with trajectories of the following:
                   feascond, gradcond, compcond, costcond, gamma,
                   stepsize, obj, alphap, alphad
           message - exit message
       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 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_ipopt(H, c, A, l, u, xmin, [], x0, opt);

   See also IPOPT, IPOPT_OPTIONS.
   http://www.coin-or.org/projects/Ipopt.xml.

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 %   Uses IPOPT to solve the following QP (quadratic programming) problem:
0006 %
0007 %       min 1/2 X'*H*X + C'*X
0008 %        X
0009 %
0010 %   subject to
0011 %
0012 %       L <= A*X <= U       (linear constraints)
0013 %       XMIN <= X <= XMAX   (variable bounds)
0014 %
0015 %   Inputs (all optional except H, C, A and L):
0016 %       H : matrix (possibly sparse) of quadratic cost coefficients
0017 %       C : vector of linear cost coefficients
0018 %       A, L, U : define the optional linear constraints. Default
0019 %           values for the elements of L and U are -Inf and Inf,
0020 %           respectively.
0021 %       XMIN, XMAX : optional lower and upper bounds on the
0022 %           X variables, defaults are -Inf and Inf, respectively.
0023 %       X0 : optional starting value of optimization vector X
0024 %       OPT : optional options structure with the following fields,
0025 %           all of which are also optional (default values shown in
0026 %           parentheses)
0027 %           verbose (0) - controls level of progress output displayed
0028 %               0 = no progress output
0029 %               1 = some progress output
0030 %               2 = verbose progress output
0031 %           ipopt_opt - options struct for IPOPT, value in verbose
0032 %               overrides these options
0033 %       PROBLEM : The inputs can alternatively be supplied in a single
0034 %           PROBLEM struct with fields corresponding to the input arguments
0035 %           described above: H, c, A, l, u, xmin, xmax, x0, opt
0036 %
0037 %   Outputs:
0038 %       X : solution vector
0039 %       F : final objective function value
0040 %       EXITFLAG : exit flag
0041 %           1 = first order optimality conditions satisfied
0042 %           0 = maximum number of iterations reached
0043 %           -1 = numerically failed
0044 %       OUTPUT : output struct with the following fields:
0045 %           iterations - number of iterations performed
0046 %           hist - struct array with trajectories of the following:
0047 %                   feascond, gradcond, compcond, costcond, gamma,
0048 %                   stepsize, obj, alphap, alphad
0049 %           message - exit message
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_ipopt(H, c, A, l, u, xmin, xmax, x0, opt)
0065 %
0066 %       x = qps_ipopt(H, c, A, l, u)
0067 %       x = qps_ipopt(H, c, A, l, u, xmin, xmax)
0068 %       x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0)
0069 %       x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt)
0070 %       x = qps_ipopt(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_ipopt(...)
0074 %       [x, f] = qps_ipopt(...)
0075 %       [x, f, exitflag] = qps_ipopt(...)
0076 %       [x, f, exitflag, output] = qps_ipopt(...)
0077 %       [x, f, exitflag, output, lambda] = qps_ipopt(...)
0078 %
0079 %   Example: (problem from from http://www.jmu.edu/docs/sasdoc/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_ipopt(H, c, A, l, u, xmin, [], x0, opt);
0093 %
0094 %   See also IPOPT, IPOPT_OPTIONS.
0095 %   http://www.coin-or.org/projects/Ipopt.xml.
0096 
0097 %   MATPOWER
0098 %   Copyright (c) 2010-2015 by Power System Engineering Research Center (PSERC)
0099 %   by Ray Zimmerman, PSERC Cornell
0100 %
0101 %   $Id: qps_ipopt.m 2661 2015-03-20 17:02:46Z ray $
0102 %
0103 %   This file is part of MATPOWER.
0104 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0105 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0106 
0107 %% check for IPOPT
0108 % if ~have_fcn('ipopt')
0109 %     error('qps_ipopt: requires IPOPT (http://www.coin-or.org/projects/Ipopt.xml)');
0110 % end
0111 
0112 %%----- input argument handling  -----
0113 %% gather inputs
0114 if nargin == 1 && isstruct(H)       %% problem struct
0115     p = H;
0116     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0117     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0118     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0119     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0120     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0121     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0122     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0123     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0124     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0125 else                                %% individual args
0126     if nargin < 9
0127         opt = [];
0128         if nargin < 8
0129             x0 = [];
0130             if nargin < 7
0131                 xmax = [];
0132                 if nargin < 6
0133                     xmin = [];
0134                 end
0135             end
0136         end
0137     end
0138 end
0139 
0140 %% define nx, set default values for missing optional inputs
0141 if isempty(H) || ~any(any(H))
0142     if isempty(A) && isempty(xmin) && isempty(xmax)
0143         error('qps_ipopt: LP problem must include constraints or variable bounds');
0144     else
0145         if ~isempty(A)
0146             nx = size(A, 2);
0147         elseif ~isempty(xmin)
0148             nx = length(xmin);
0149         else    % if ~isempty(xmax)
0150             nx = length(xmax);
0151         end
0152     end
0153     H = sparse(nx,nx);
0154 else
0155     nx = size(H, 1);
0156 end
0157 if isempty(c)
0158     c = zeros(nx, 1);
0159 end
0160 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0161                                  (isempty(u) || all(u == Inf)))
0162     A = sparse(0,nx);           %% no limits => no linear constraints
0163 end
0164 nA = size(A, 1);                %% number of original linear constraints
0165 if nA
0166     if isempty(u)               %% By default, linear inequalities are ...
0167         u = Inf(nA, 1);             %% ... unbounded above and ...
0168     end
0169     if isempty(l)
0170         l = -Inf(nA, 1);            %% ... unbounded below.
0171     end
0172 end
0173 if isempty(x0)
0174     x0 = zeros(nx, 1);
0175 end
0176 
0177 %% default options
0178 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0179     verbose = opt.verbose;
0180 else
0181     verbose = 0;
0182 end
0183 
0184 %% make sure args are sparse/full as expected by IPOPT
0185 if ~isempty(H)
0186     if ~issparse(H)
0187         H = sparse(H);
0188     end
0189 end
0190 if ~issparse(A)
0191     A = sparse(A);
0192 end
0193 
0194 
0195 %%-----  run optimization  -----
0196 %% set options struct for IPOPT
0197 if ~isempty(opt) && isfield(opt, 'ipopt_opt') && ~isempty(opt.ipopt_opt)
0198     options.ipopt = ipopt_options(opt.ipopt_opt);
0199 else
0200     options.ipopt = ipopt_options;
0201 end
0202 options.ipopt.jac_c_constant    = 'yes';
0203 options.ipopt.jac_d_constant    = 'yes';
0204 options.ipopt.hessian_constant  = 'yes';
0205 options.ipopt.least_square_init_primal  = 'yes';
0206 options.ipopt.least_square_init_duals   = 'yes';
0207 % options.ipopt.mehrotra_algorithm        = 'yes';        %% default 'no'
0208 if verbose
0209     options.ipopt.print_level = min(12, verbose*2+1);
0210 else
0211     options.ipopt.print_level = 0;
0212 end
0213 
0214 %% define variable and constraint bounds, if given
0215 if nA
0216     options.cu = u;
0217     options.cl = l;
0218 end
0219 if ~isempty(xmin)
0220     options.lb = xmin;
0221 end
0222 if ~isempty(xmax)
0223     options.ub = xmax;
0224 end
0225 
0226 %% assign function handles
0227 funcs.objective         = @(x) 0.5 * x' * H * x + c' * x;
0228 funcs.gradient          = @(x) H * x + c;
0229 funcs.constraints       = @(x) A * x;
0230 funcs.jacobian          = @(x) A;
0231 funcs.jacobianstructure = @() A;
0232 funcs.hessian           = @(x, sigma, lambda) tril(H);
0233 funcs.hessianstructure  = @() tril(H);
0234 
0235 %% run the optimization
0236 [x, info] = ipopt(x0,funcs,options);
0237 
0238 if info.status == 0 || info.status == 1
0239     eflag = 1;
0240 else
0241     eflag = 0;
0242 end
0243 if isfield(info, 'iter')
0244     output.iterations = info.iter;
0245 end
0246 output.info       = info.status;
0247 f = funcs.objective(x);
0248 
0249 %% check for empty results (in case optimization failed)
0250 if isempty(info.lambda)
0251     lam = NaN(nA, 1);
0252 else
0253     lam = info.lambda;
0254 end
0255 if isempty(info.zl) && ~isempty(xmin)
0256     zl = NaN(nx, 1);
0257 else
0258     zl = info.zl;
0259 end
0260 if isempty(info.zu) && ~isempty(xmax)
0261     zu = NaN(nx, 1);
0262 else
0263     zu = info.zu;
0264 end
0265 
0266 %% repackage lambdas
0267 kl = find(lam < 0);                     %% lower bound binding
0268 ku = find(lam > 0);                     %% upper bound binding
0269 mu_l = zeros(nA, 1);
0270 mu_l(kl) = -lam(kl);
0271 mu_u = zeros(nA, 1);
0272 mu_u(ku) = lam(ku);
0273 
0274 lambda = struct( ...
0275     'mu_l', mu_l, ...
0276     'mu_u', mu_u, ...
0277     'lower', zl, ...
0278     'upper', zu    );

Generated on Fri 20-Mar-2015 18:23:34 by m2html © 2005