Home > matpower6.0 > 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-2016, Power Systems Engineering Research Center (PSERC)
0099 %   by Ray Zimmerman, PSERC Cornell
0100 %
0101 %   This file is part of MATPOWER.
0102 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0103 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0104 
0105 %% check for IPOPT
0106 % if ~have_fcn('ipopt')
0107 %     error('qps_ipopt: requires IPOPT (http://www.coin-or.org/projects/Ipopt.xml)');
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     if isempty(A) && isempty(xmin) && isempty(xmax)
0141         error('qps_ipopt: LP problem must include constraints or variable bounds');
0142     else
0143         if ~isempty(A)
0144             nx = size(A, 2);
0145         elseif ~isempty(xmin)
0146             nx = length(xmin);
0147         else    % if ~isempty(xmax)
0148             nx = length(xmax);
0149         end
0150     end
0151     H = sparse(nx,nx);
0152 else
0153     nx = size(H, 1);
0154 end
0155 if isempty(c)
0156     c = zeros(nx, 1);
0157 end
0158 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0159                                  (isempty(u) || all(u == Inf)))
0160     A = sparse(0,nx);           %% no limits => no linear constraints
0161 end
0162 nA = size(A, 1);                %% number of original linear constraints
0163 if nA
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 end
0171 if isempty(x0)
0172     x0 = zeros(nx, 1);
0173 end
0174 
0175 %% default options
0176 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0177     verbose = opt.verbose;
0178 else
0179     verbose = 0;
0180 end
0181 
0182 %% make sure args are sparse/full as expected by IPOPT
0183 if ~isempty(H)
0184     if ~issparse(H)
0185         H = sparse(H);
0186     end
0187 end
0188 if ~issparse(A)
0189     A = sparse(A);
0190 end
0191 
0192 
0193 %%-----  run optimization  -----
0194 %% set options struct for IPOPT
0195 if ~isempty(opt) && isfield(opt, 'ipopt_opt') && ~isempty(opt.ipopt_opt)
0196     options.ipopt = ipopt_options(opt.ipopt_opt);
0197 else
0198     options.ipopt = ipopt_options;
0199 end
0200 options.ipopt.jac_c_constant    = 'yes';
0201 options.ipopt.jac_d_constant    = 'yes';
0202 options.ipopt.hessian_constant  = 'yes';
0203 options.ipopt.least_square_init_primal  = 'yes';
0204 options.ipopt.least_square_init_duals   = 'yes';
0205 % options.ipopt.mehrotra_algorithm        = 'yes';        %% default 'no'
0206 if verbose
0207     options.ipopt.print_level = min(12, verbose*2+1);
0208 else
0209     options.ipopt.print_level = 0;
0210 end
0211 
0212 %% define variable and constraint bounds, if given
0213 if nA
0214     options.cu = u;
0215     options.cl = l;
0216 end
0217 if ~isempty(xmin)
0218     options.lb = xmin;
0219 end
0220 if ~isempty(xmax)
0221     options.ub = xmax;
0222 end
0223 
0224 %% assign function handles
0225 funcs.objective         = @(x) 0.5 * x' * H * x + c' * x;
0226 funcs.gradient          = @(x) H * x + c;
0227 funcs.constraints       = @(x) A * x;
0228 funcs.jacobian          = @(x) A;
0229 funcs.jacobianstructure = @() A;
0230 funcs.hessian           = @(x, sigma, lambda) tril(H);
0231 funcs.hessianstructure  = @() tril(H);
0232 
0233 %% run the optimization
0234 [x, info] = ipopt(x0,funcs,options);
0235 
0236 if info.status == 0 || info.status == 1
0237     eflag = 1;
0238 else
0239     eflag = 0;
0240 end
0241 if isfield(info, 'iter')
0242     output.iterations = info.iter;
0243 end
0244 output.info       = info.status;
0245 f = funcs.objective(x);
0246 
0247 %% check for empty results (in case optimization failed)
0248 if isempty(info.lambda)
0249     lam = NaN(nA, 1);
0250 else
0251     lam = info.lambda;
0252 end
0253 if isempty(info.zl) && ~isempty(xmin)
0254     zl = NaN(nx, 1);
0255 else
0256     zl = info.zl;
0257 end
0258 if isempty(info.zu) && ~isempty(xmax)
0259     zu = NaN(nx, 1);
0260 else
0261     zu = info.zu;
0262 end
0263 
0264 %% repackage lambdas
0265 kl = find(lam < 0);                     %% lower bound binding
0266 ku = find(lam > 0);                     %% upper bound binding
0267 mu_l = zeros(nA, 1);
0268 mu_l(kl) = -lam(kl);
0269 mu_u = zeros(nA, 1);
0270 mu_u(ku) = lam(ku);
0271 
0272 lambda = struct( ...
0273     'mu_l', mu_l, ...
0274     'mu_u', mu_u, ...
0275     'lower', zl, ...
0276     'upper', zu    );

Generated on Fri 16-Dec-2016 12:45:37 by m2html © 2005