Home > matpower5.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.
   https://projects.coin-or.org/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 %   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 %   https://projects.coin-or.org/Ipopt/.
0096 
0097 %   MATPOWER
0098 %   $Id: qps_ipopt.m 2340 2014-06-27 19:15:10Z ray $
0099 %   by Ray Zimmerman, PSERC Cornell
0100 %   Copyright (c) 2010 by Power System Engineering Research Center (PSERC)
0101 %
0102 %   This file is part of MATPOWER.
0103 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0104 %
0105 %   MATPOWER is free software: you can redistribute it and/or modify
0106 %   it under the terms of the GNU General Public License as published
0107 %   by the Free Software Foundation, either version 3 of the License,
0108 %   or (at your option) any later version.
0109 %
0110 %   MATPOWER is distributed in the hope that it will be useful,
0111 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0112 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0113 %   GNU General Public License for more details.
0114 %
0115 %   You should have received a copy of the GNU General Public License
0116 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0117 %
0118 %   Additional permission under GNU GPL version 3 section 7
0119 %
0120 %   If you modify MATPOWER, or any covered work, to interface with
0121 %   other modules (such as MATLAB code and MEX-files) available in a
0122 %   MATLAB(R) or comparable environment containing parts covered
0123 %   under other licensing terms, the licensors of MATPOWER grant
0124 %   you additional permission to convey the resulting work.
0125 
0126 %% check for IPOPT
0127 % if ~have_fcn('ipopt')
0128 %     error('qps_ipopt: requires IPOPT (https://projects.coin-or.org/Ipopt/)');
0129 % end
0130 
0131 %%----- input argument handling  -----
0132 %% gather inputs
0133 if nargin == 1 && isstruct(H)       %% problem struct
0134     p = H;
0135     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0136     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0137     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0138     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0139     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0140     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0141     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0142     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0143     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0144 else                                %% individual args
0145     if nargin < 9
0146         opt = [];
0147         if nargin < 8
0148             x0 = [];
0149             if nargin < 7
0150                 xmax = [];
0151                 if nargin < 6
0152                     xmin = [];
0153                 end
0154             end
0155         end
0156     end
0157 end
0158 
0159 %% define nx, set default values for missing optional inputs
0160 if isempty(H) || ~any(any(H))
0161     if isempty(A) && isempty(xmin) && isempty(xmax)
0162         error('qps_ipopt: LP problem must include constraints or variable bounds');
0163     else
0164         if ~isempty(A)
0165             nx = size(A, 2);
0166         elseif ~isempty(xmin)
0167             nx = length(xmin);
0168         else    % if ~isempty(xmax)
0169             nx = length(xmax);
0170         end
0171     end
0172     H = sparse(nx,nx);
0173 else
0174     nx = size(H, 1);
0175 end
0176 if isempty(c)
0177     c = zeros(nx, 1);
0178 end
0179 if  ~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0180                    (isempty(u) || all(u == Inf))
0181     A = sparse(0,nx);           %% no limits => no linear constraints
0182 end
0183 nA = size(A, 1);                %% number of original linear constraints
0184 if nA
0185     if isempty(u)               %% By default, linear inequalities are ...
0186         u = Inf * ones(nA, 1);      %% ... unbounded above and ...
0187     end
0188     if isempty(l)
0189         l = -Inf * ones(nA, 1);     %% ... unbounded below.
0190     end
0191 end
0192 if isempty(x0)
0193     x0 = zeros(nx, 1);
0194 end
0195 
0196 %% default options
0197 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0198     verbose = opt.verbose;
0199 else
0200     verbose = 0;
0201 end
0202 
0203 %% make sure args are sparse/full as expected by IPOPT
0204 if ~isempty(H)
0205     if ~issparse(H)
0206         H = sparse(H);
0207     end
0208 end
0209 if ~issparse(A)
0210     A = sparse(A);
0211 end
0212 
0213 
0214 %%-----  run optimization  -----
0215 %% set options struct for IPOPT
0216 if ~isempty(opt) && isfield(opt, 'ipopt_opt') && ~isempty(opt.ipopt_opt)
0217     options.ipopt = ipopt_options(opt.ipopt_opt);
0218 else
0219     options.ipopt = ipopt_options;
0220 end
0221 options.ipopt.jac_c_constant    = 'yes';
0222 options.ipopt.jac_d_constant    = 'yes';
0223 options.ipopt.hessian_constant  = 'yes';
0224 options.ipopt.least_square_init_primal  = 'yes';
0225 options.ipopt.least_square_init_duals   = 'yes';
0226 % options.ipopt.mehrotra_algorithm        = 'yes';        %% default 'no'
0227 if verbose
0228     options.ipopt.print_level = min(12, verbose*2+1);
0229 else
0230     options.ipopt.print_level = 0;
0231 end
0232 
0233 %% define variable and constraint bounds, if given
0234 if nA
0235     options.cu = u;
0236     options.cl = l;
0237 end
0238 if ~isempty(xmin)
0239     options.lb = xmin;
0240 end
0241 if ~isempty(xmax)
0242     options.ub = xmax;
0243 end
0244 
0245 %% assign function handles
0246 funcs.objective         = @(x) 0.5 * x' * H * x + c' * x;
0247 funcs.gradient          = @(x) H * x + c;
0248 funcs.constraints       = @(x) A * x;
0249 funcs.jacobian          = @(x) A;
0250 funcs.jacobianstructure = @() A;
0251 funcs.hessian           = @(x, sigma, lambda) tril(H);
0252 funcs.hessianstructure  = @() tril(H);
0253 
0254 %% run the optimization
0255 [x, info] = ipopt(x0,funcs,options);
0256 
0257 if info.status == 0 || info.status == 1
0258     eflag = 1;
0259 else
0260     eflag = 0;
0261 end
0262 if isfield(info, 'iter')
0263     output.iterations = info.iter;
0264 end
0265 output.info       = info.status;
0266 f = funcs.objective(x);
0267 
0268 %% check for empty results (in case optimization failed)
0269 if isempty(info.lambda)
0270     lam = NaN(nA, 1);
0271 else
0272     lam = info.lambda;
0273 end
0274 if isempty(info.zl) && ~isempty(xmin)
0275     zl = NaN(nx, 1);
0276 else
0277     zl = info.zl;
0278 end
0279 if isempty(info.zu) && ~isempty(xmax)
0280     zu = NaN(nx, 1);
0281 else
0282     zu = info.zu;
0283 end
0284 
0285 %% repackage lambdas
0286 kl = find(lam < 0);                     %% lower bound binding
0287 ku = find(lam > 0);                     %% upper bound binding
0288 mu_l = zeros(nA, 1);
0289 mu_l(kl) = -lam(kl);
0290 mu_u = zeros(nA, 1);
0291 mu_u(ku) = lam(ku);
0292 
0293 lambda = struct( ...
0294     'mu_l', mu_l, ...
0295     'mu_u', mu_u, ...
0296     'lower', zl, ...
0297     'upper', zu    );

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005