Home > matpower4.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
           max_it (0) - maximum number of iterations allowed
               0 = use algorithm default
           ipopt_opt - options struct for IPOPT, values in
               verbose and max_it override 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 %           max_it (0) - maximum number of iterations allowed
0032 %               0 = use algorithm default
0033 %           ipopt_opt - options struct for IPOPT, values in
0034 %               verbose and max_it override 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 = first order optimality conditions satisfied
0044 %           0 = maximum number of iterations reached
0045 %           -1 = numerically failed
0046 %       OUTPUT : output struct with the following fields:
0047 %           iterations - number of iterations performed
0048 %           hist - struct array with trajectories of the following:
0049 %                   feascond, gradcond, compcond, costcond, gamma,
0050 %                   stepsize, obj, alphap, alphad
0051 %           message - exit message
0052 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0053 %           multipliers on the constraints, with fields:
0054 %           mu_l - lower (left-hand) limit on linear constraints
0055 %           mu_u - upper (right-hand) limit on linear constraints
0056 %           lower - lower bound on optimization variables
0057 %           upper - upper bound on optimization variables
0058 %
0059 %   Note the calling syntax is almost identical to that of QUADPROG
0060 %   from MathWorks' Optimization Toolbox. The main difference is that
0061 %   the linear constraints are specified with A, L, U instead of
0062 %   A, B, Aeq, Beq.
0063 %
0064 %   Calling syntax options:
0065 %       [x, f, exitflag, output, lambda] = ...
0066 %           qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt)
0067 %
0068 %       x = qps_ipopt(H, c, A, l, u)
0069 %       x = qps_ipopt(H, c, A, l, u, xmin, xmax)
0070 %       x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0)
0071 %       x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt)
0072 %       x = qps_ipopt(problem), where problem is a struct with fields:
0073 %                       H, c, A, l, u, xmin, xmax, x0, opt
0074 %                       all fields except 'c', 'A' and 'l' or 'u' are optional
0075 %       x = qps_ipopt(...)
0076 %       [x, f] = qps_ipopt(...)
0077 %       [x, f, exitflag] = qps_ipopt(...)
0078 %       [x, f, exitflag, output] = qps_ipopt(...)
0079 %       [x, f, exitflag, output, lambda] = qps_ipopt(...)
0080 %
0081 %   Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm)
0082 %       H = [   1003.1  4.3     6.3     5.9;
0083 %               4.3     2.2     2.1     3.9;
0084 %               6.3     2.1     3.5     4.8;
0085 %               5.9     3.9     4.8     10  ];
0086 %       c = zeros(4,1);
0087 %       A = [   1       1       1       1;
0088 %               0.17    0.11    0.10    0.18    ];
0089 %       l = [1; 0.10];
0090 %       u = [1; Inf];
0091 %       xmin = zeros(4,1);
0092 %       x0 = [1; 0; 0; 1];
0093 %       opt = struct('verbose', 2);
0094 %       [x, f, s, out, lambda] = qps_ipopt(H, c, A, l, u, xmin, [], x0, opt);
0095 %
0096 %   See also IPOPT, IPOPT_OPTIONS.
0097 %   https://projects.coin-or.org/Ipopt/.
0098 
0099 %   MATPOWER
0100 %   $Id: qps_ipopt.m,v 1.6 2011/09/09 15:26:08 cvs Exp $
0101 %   by Ray Zimmerman, PSERC Cornell
0102 %   Copyright (c) 2010 by Power System Engineering Research Center (PSERC)
0103 %
0104 %   This file is part of MATPOWER.
0105 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0106 %
0107 %   MATPOWER is free software: you can redistribute it and/or modify
0108 %   it under the terms of the GNU General Public License as published
0109 %   by the Free Software Foundation, either version 3 of the License,
0110 %   or (at your option) any later version.
0111 %
0112 %   MATPOWER is distributed in the hope that it will be useful,
0113 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0114 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0115 %   GNU General Public License for more details.
0116 %
0117 %   You should have received a copy of the GNU General Public License
0118 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0119 %
0120 %   Additional permission under GNU GPL version 3 section 7
0121 %
0122 %   If you modify MATPOWER, or any covered work, to interface with
0123 %   other modules (such as MATLAB code and MEX-files) available in a
0124 %   MATLAB(R) or comparable environment containing parts covered
0125 %   under other licensing terms, the licensors of MATPOWER grant
0126 %   you additional permission to convey the resulting work.
0127 
0128 %% check for IPOPT
0129 % if ~have_fcn('ipopt')
0130 %     error('qps_ipopt: requires IPOPT (https://projects.coin-or.org/Ipopt/)');
0131 % end
0132 
0133 %%----- input argument handling  -----
0134 %% gather inputs
0135 if nargin == 1 && isstruct(H)       %% problem struct
0136     p = H;
0137     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0138     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0139     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0140     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0141     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0142     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0143     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0144     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0145     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0146 else                                %% individual args
0147     if nargin < 9
0148         opt = [];
0149         if nargin < 8
0150             x0 = [];
0151             if nargin < 7
0152                 xmax = [];
0153                 if nargin < 6
0154                     xmin = [];
0155                 end
0156             end
0157         end
0158     end
0159 end
0160 
0161 %% define nx, set default values for missing optional inputs
0162 if isempty(H) || ~any(any(H))
0163     if isempty(A) && isempty(xmin) && isempty(xmax)
0164         error('qps_ipopt: LP problem must include constraints or variable bounds');
0165     else
0166         if ~isempty(A)
0167             nx = size(A, 2);
0168         elseif ~isempty(xmin)
0169             nx = length(xmin);
0170         else    % if ~isempty(xmax)
0171             nx = length(xmax);
0172         end
0173     end
0174     H = sparse(nx,nx);
0175 else
0176     nx = size(H, 1);
0177 end
0178 if isempty(c)
0179     c = zeros(nx, 1);
0180 end
0181 if  ~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0182                    (isempty(u) || all(u == Inf))
0183     A = sparse(0,nx);           %% no limits => no linear constraints
0184 end
0185 nA = size(A, 1);                %% number of original linear constraints
0186 if nA
0187     if isempty(u)               %% By default, linear inequalities are ...
0188         u = Inf * ones(nA, 1);      %% ... unbounded above and ...
0189     end
0190     if isempty(l)
0191         l = -Inf * ones(nA, 1);     %% ... unbounded below.
0192     end
0193 end
0194 if isempty(x0)
0195     x0 = zeros(nx, 1);
0196 end
0197 
0198 %% default options
0199 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0200     verbose = opt.verbose;
0201 else
0202     verbose = 0;
0203 end
0204 if ~isempty(opt) && isfield(opt, 'max_it') && ~isempty(opt.max_it)
0205     max_it = opt.max_it;
0206 else
0207     max_it = 0;
0208 end
0209 
0210 %% make sure args are sparse/full as expected by IPOPT
0211 if ~isempty(H)
0212     if ~issparse(H)
0213         H = sparse(H);
0214     end
0215 end
0216 if ~issparse(A)
0217     A = sparse(A);
0218 end
0219 
0220 
0221 %%-----  run optimization  -----
0222 %% set options struct for IPOPT
0223 if ~isempty(opt) && isfield(opt, 'ipopt_opt') && ~isempty(opt.ipopt_opt)
0224     options.ipopt = ipopt_options(opt.ipopt_opt);
0225 else
0226     options.ipopt = ipopt_options;
0227 end
0228 options.ipopt.jac_c_constant    = 'yes';
0229 options.ipopt.jac_d_constant    = 'yes';
0230 options.ipopt.hessian_constant  = 'yes';
0231 options.ipopt.least_square_init_primal  = 'yes';
0232 options.ipopt.least_square_init_duals   = 'yes';
0233 % options.ipopt.mehrotra_algorithm        = 'yes';        %% default 'no'
0234 if verbose
0235     options.ipopt.print_level = min(12, verbose*2+1);
0236 else
0237     options.ipopt.print_level = 0;
0238 end
0239 if max_it
0240     options.ipopt.max_iter = max_it;
0241 end
0242 
0243 %% define variable and constraint bounds, if given
0244 if nA
0245     options.cu = u;
0246     options.cl = l;
0247 end
0248 if ~isempty(xmin)
0249     options.lb = xmin;
0250 end
0251 if ~isempty(xmax)
0252     options.ub = xmax;
0253 end
0254 
0255 %% assign function handles
0256 funcs.objective         = @(x) 0.5 * x' * H * x + c' * x;
0257 funcs.gradient          = @(x) H * x + c;
0258 funcs.constraints       = @(x) A * x;
0259 funcs.jacobian          = @(x) A;
0260 funcs.jacobianstructure = @() A;
0261 funcs.hessian           = @(x, sigma, lambda) tril(H);
0262 funcs.hessianstructure  = @() tril(H);
0263 
0264 %% run the optimization
0265 [x, info] = ipopt(x0,funcs,options);
0266 
0267 if info.status == 0 || info.status == 1
0268     eflag = 1;
0269 else
0270     eflag = 0;
0271 end
0272 if isfield(info, 'iter')
0273     output.iterations = info.iter;
0274 end
0275 output.info       = info.status;
0276 f = funcs.objective(x);
0277 
0278 %% repackage lambdas
0279 kl = find(info.lambda < 0);                     %% lower bound binding
0280 ku = find(info.lambda > 0);                     %% upper bound binding
0281 mu_l = zeros(nA, 1);
0282 mu_l(kl) = -info.lambda(kl);
0283 mu_u = zeros(nA, 1);
0284 mu_u(ku) = info.lambda(ku);
0285 
0286 lambda = struct( ...
0287     'mu_l', mu_l, ...
0288     'mu_u', mu_u, ...
0289     'lower', info.zl, ...
0290     'upper', info.zu    );

Generated on Mon 26-Jan-2015 15:00:13 by m2html © 2005