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

nlps_ipopt

PURPOSE ^

NLPS_IPOPT Nonlinear programming (NLP) Solver based on IPOPT.

SYNOPSIS ^

function [x, f, eflag, output, lambda] = nlps_ipopt(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt)

DESCRIPTION ^

NLPS_IPOPT  Nonlinear programming (NLP) Solver based on IPOPT.
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
       NLPS_IPOPT(F_FCN, X0, A, L, U, XMIN, XMAX, GH_FCN, HESS_FCN, OPT)
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = NLPS_IPOPT(PROBLEM)
   A wrapper function providing a standardized interface for using
   IPOPT to solve the following NLP (nonlinear programming) problem:

   Minimize a function F(X) beginning from a starting point X0, subject to
   optional linear and nonlinear constraints and variable bounds.

       min F(X)
        X

   subject to

       G(X) = 0            (nonlinear equalities)
       H(X) <= 0           (nonlinear inequalities)
       L <= A*X <= U       (linear constraints)
       XMIN <= X <= XMAX   (variable bounds)

   Inputs (all optional except F_FCN and X0):
       F_FCN : handle to function that evaluates the objective function,
           its gradients and Hessian for a given value of X. If there
           are nonlinear constraints, the Hessian information is
           provided by the HESS_FCN function passed in the 9th argument
           and is not required here. Calling syntax for this function:
               [F, DF, D2F] = F_FCN(X)
       X0 : starting value of optimization vector X
       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.
       GH_FCN : handle to function that evaluates the optional
           nonlinear constraints and their gradients for a given
           value of X. Calling syntax for this function is:
               [H, G, DH, DG] = GH_FCN(X)
           where the columns of DH and DG are the gradients of the
           corresponding elements of H and G, i.e. DH and DG are
           transposes of the Jacobians of H and G, respectively.
       HESS_FCN : handle to function that computes the Hessian of the
           Lagrangian for given values of X, lambda and mu, where
           lambda and mu are the multipliers on the equality and
           inequality constraints, g and h, respectively. The calling
           syntax for this function is:
               LXX = HESS_FCN(X, LAM)
           where lambda = LAM.eqnonlin and mu = LAM.ineqnonlin.
       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: f_fcn, x0, A, l, u, xmin, xmax,
                            gh_fcn, hess_fcn, 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:
           eqnonlin - nonlinear equality constraints
           ineqnonlin - nonlinear inequality constraints
           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 FMINCON
   from MathWorks' Optimization Toolbox. The main difference is that
   the linear constraints are specified with A, L, U instead of
   A, B, Aeq, Beq. The functions for evaluating the objective
   function, constraints and Hessian are identical.

   Calling syntax options:
       [x, f, exitflag, output, lambda] = ...
           nlps_ipopt(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);

       x = nlps_ipopt(f_fcn, x0);
       x = nlps_ipopt(f_fcn, x0, A, l);
       x = nlps_ipopt(f_fcn, x0, A, l, u);
       x = nlps_ipopt(f_fcn, x0, A, l, u, xmin);
       x = nlps_ipopt(f_fcn, x0, A, l, u, xmin, xmax);
       x = nlps_ipopt(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn);
       x = nlps_ipopt(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn);
       x = nlps_ipopt(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);
       x = nlps_ipopt(problem);
               where problem is a struct with fields:
                   f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt
                   all fields except 'f_fcn' and 'x0' are optional
       x = nlps_ipopt(...);
       [x, f] = nlps_ipopt(...);
       [x, f, exitflag] = nlps_ipopt(...);
       [x, f, exitflag, output] = nlps_ipopt(...);
       [x, f, exitflag, output, lambda] = nlps_ipopt(...);

   Example: (problem from https://en.wikipedia.org/wiki/Nonlinear_programming)
       function [f, df, d2f] = f2(x)
       f = -x(1)*x(2) - x(2)*x(3);
       if nargout > 1           %% gradient is required
           df = -[x(2); x(1)+x(3); x(2)];
           if nargout > 2       %% Hessian is required
               d2f = -[0 1 0; 1 0 1; 0 1 0];   %% actually not used since
           end                                 %% 'hess_fcn' is provided
       end
       
       function [h, g, dh, dg] = gh2(x)
       h = [ 1 -1 1; 1 1 1] * x.^2 + [-2; -10];
       dh = 2 * [x(1) x(1); -x(2) x(2); x(3) x(3)];
       g = []; dg = [];
       
       function Lxx = hess2(x, lam, cost_mult)
       if nargin < 3, cost_mult = 1; end
       mu = lam.ineqnonlin;
       Lxx = cost_mult * [0 -1 0; -1 0 -1; 0 -1 0] + ...
               [2*[1 1]*mu 0 0; 0 2*[-1 1]*mu 0; 0 0 2*[1 1]*mu];
       
       problem = struct( ...
           'f_fcn',    @(x)f2(x), ...
           'gh_fcn',   @(x)gh2(x), ...
           'hess_fcn', @(x, lam, cost_mult)hess2(x, lam, cost_mult), ...
           'x0',       [1; 1; 0], ...
           'opt',      struct('verbose', 2) ...
       );
       [x, f, exitflag, output, lambda] = nlps_ipopt(problem);

   See also NLPS_MASTER, IPOPT.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [x, f, eflag, output, lambda] = nlps_ipopt(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt)
0002 %NLPS_IPOPT  Nonlinear programming (NLP) Solver based on IPOPT.
0003 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
0004 %       NLPS_IPOPT(F_FCN, X0, A, L, U, XMIN, XMAX, GH_FCN, HESS_FCN, OPT)
0005 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = NLPS_IPOPT(PROBLEM)
0006 %   A wrapper function providing a standardized interface for using
0007 %   IPOPT to solve the following NLP (nonlinear programming) problem:
0008 %
0009 %   Minimize a function F(X) beginning from a starting point X0, subject to
0010 %   optional linear and nonlinear constraints and variable bounds.
0011 %
0012 %       min F(X)
0013 %        X
0014 %
0015 %   subject to
0016 %
0017 %       G(X) = 0            (nonlinear equalities)
0018 %       H(X) <= 0           (nonlinear inequalities)
0019 %       L <= A*X <= U       (linear constraints)
0020 %       XMIN <= X <= XMAX   (variable bounds)
0021 %
0022 %   Inputs (all optional except F_FCN and X0):
0023 %       F_FCN : handle to function that evaluates the objective function,
0024 %           its gradients and Hessian for a given value of X. If there
0025 %           are nonlinear constraints, the Hessian information is
0026 %           provided by the HESS_FCN function passed in the 9th argument
0027 %           and is not required here. Calling syntax for this function:
0028 %               [F, DF, D2F] = F_FCN(X)
0029 %       X0 : starting value of optimization vector X
0030 %       A, L, U : define the optional linear constraints. Default
0031 %           values for the elements of L and U are -Inf and Inf,
0032 %           respectively.
0033 %       XMIN, XMAX : optional lower and upper bounds on the
0034 %           X variables, defaults are -Inf and Inf, respectively.
0035 %       GH_FCN : handle to function that evaluates the optional
0036 %           nonlinear constraints and their gradients for a given
0037 %           value of X. Calling syntax for this function is:
0038 %               [H, G, DH, DG] = GH_FCN(X)
0039 %           where the columns of DH and DG are the gradients of the
0040 %           corresponding elements of H and G, i.e. DH and DG are
0041 %           transposes of the Jacobians of H and G, respectively.
0042 %       HESS_FCN : handle to function that computes the Hessian of the
0043 %           Lagrangian for given values of X, lambda and mu, where
0044 %           lambda and mu are the multipliers on the equality and
0045 %           inequality constraints, g and h, respectively. The calling
0046 %           syntax for this function is:
0047 %               LXX = HESS_FCN(X, LAM)
0048 %           where lambda = LAM.eqnonlin and mu = LAM.ineqnonlin.
0049 %       OPT : optional options structure with the following fields,
0050 %           all of which are also optional (default values shown in
0051 %           parentheses)
0052 %           verbose (0) - controls level of progress output displayed
0053 %               0 = no progress output
0054 %               1 = some progress output
0055 %               2 = verbose progress output
0056 %           ipopt_opt   - options struct for IPOPT, value in verbose
0057 %                   overrides these options
0058 %       PROBLEM : The inputs can alternatively be supplied in a single
0059 %           PROBLEM struct with fields corresponding to the input arguments
0060 %           described above: f_fcn, x0, A, l, u, xmin, xmax,
0061 %                            gh_fcn, hess_fcn, opt
0062 %
0063 %   Outputs:
0064 %       X : solution vector
0065 %       F : final objective function value
0066 %       EXITFLAG : exit flag
0067 %           1 = converged
0068 %           0 = failed to converge
0069 %       OUTPUT : output struct with the following fields:
0070 %           status - see IPOPT documentation for INFO.status
0071 %             https://coin-or.github.io/Ipopt/IpReturnCodes__inc_8h_source.html
0072 %           iterations - number of iterations performed (INFO.iter)
0073 %           cpu - see IPOPT documentation for INFO.cpu
0074 %           eval - see IPOPT documentation for INFO.eval
0075 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0076 %           multipliers on the constraints, with fields:
0077 %           eqnonlin - nonlinear equality constraints
0078 %           ineqnonlin - nonlinear inequality constraints
0079 %           mu_l - lower (left-hand) limit on linear constraints
0080 %           mu_u - upper (right-hand) limit on linear constraints
0081 %           lower - lower bound on optimization variables
0082 %           upper - upper bound on optimization variables
0083 %
0084 %   Note the calling syntax is almost identical to that of FMINCON
0085 %   from MathWorks' Optimization Toolbox. The main difference is that
0086 %   the linear constraints are specified with A, L, U instead of
0087 %   A, B, Aeq, Beq. The functions for evaluating the objective
0088 %   function, constraints and Hessian are identical.
0089 %
0090 %   Calling syntax options:
0091 %       [x, f, exitflag, output, lambda] = ...
0092 %           nlps_ipopt(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);
0093 %
0094 %       x = nlps_ipopt(f_fcn, x0);
0095 %       x = nlps_ipopt(f_fcn, x0, A, l);
0096 %       x = nlps_ipopt(f_fcn, x0, A, l, u);
0097 %       x = nlps_ipopt(f_fcn, x0, A, l, u, xmin);
0098 %       x = nlps_ipopt(f_fcn, x0, A, l, u, xmin, xmax);
0099 %       x = nlps_ipopt(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn);
0100 %       x = nlps_ipopt(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn);
0101 %       x = nlps_ipopt(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);
0102 %       x = nlps_ipopt(problem);
0103 %               where problem is a struct with fields:
0104 %                   f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt
0105 %                   all fields except 'f_fcn' and 'x0' are optional
0106 %       x = nlps_ipopt(...);
0107 %       [x, f] = nlps_ipopt(...);
0108 %       [x, f, exitflag] = nlps_ipopt(...);
0109 %       [x, f, exitflag, output] = nlps_ipopt(...);
0110 %       [x, f, exitflag, output, lambda] = nlps_ipopt(...);
0111 %
0112 %   Example: (problem from https://en.wikipedia.org/wiki/Nonlinear_programming)
0113 %       function [f, df, d2f] = f2(x)
0114 %       f = -x(1)*x(2) - x(2)*x(3);
0115 %       if nargout > 1           %% gradient is required
0116 %           df = -[x(2); x(1)+x(3); x(2)];
0117 %           if nargout > 2       %% Hessian is required
0118 %               d2f = -[0 1 0; 1 0 1; 0 1 0];   %% actually not used since
0119 %           end                                 %% 'hess_fcn' is provided
0120 %       end
0121 %
0122 %       function [h, g, dh, dg] = gh2(x)
0123 %       h = [ 1 -1 1; 1 1 1] * x.^2 + [-2; -10];
0124 %       dh = 2 * [x(1) x(1); -x(2) x(2); x(3) x(3)];
0125 %       g = []; dg = [];
0126 %
0127 %       function Lxx = hess2(x, lam, cost_mult)
0128 %       if nargin < 3, cost_mult = 1; end
0129 %       mu = lam.ineqnonlin;
0130 %       Lxx = cost_mult * [0 -1 0; -1 0 -1; 0 -1 0] + ...
0131 %               [2*[1 1]*mu 0 0; 0 2*[-1 1]*mu 0; 0 0 2*[1 1]*mu];
0132 %
0133 %       problem = struct( ...
0134 %           'f_fcn',    @(x)f2(x), ...
0135 %           'gh_fcn',   @(x)gh2(x), ...
0136 %           'hess_fcn', @(x, lam, cost_mult)hess2(x, lam, cost_mult), ...
0137 %           'x0',       [1; 1; 0], ...
0138 %           'opt',      struct('verbose', 2) ...
0139 %       );
0140 %       [x, f, exitflag, output, lambda] = nlps_ipopt(problem);
0141 %
0142 %   See also NLPS_MASTER, IPOPT.
0143 
0144 %   MP-Opt-Model
0145 %   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
0146 %   by Ray Zimmerman, PSERC Cornell
0147 %
0148 %   This file is part of MP-Opt-Model.
0149 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0150 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0151 
0152 %%----- input argument handling  -----
0153 %% gather inputs
0154 if nargin == 1 && isstruct(f_fcn)       %% problem struct
0155     p = f_fcn;
0156     f_fcn = p.f_fcn;
0157     x0 = p.x0;
0158     nx = size(x0, 1);       %% number of optimization variables
0159     if isfield(p, 'opt'),       opt = p.opt;            else,   opt = [];       end
0160     if isfield(p, 'hess_fcn'),  hess_fcn = p.hess_fcn;  else,   hess_fcn = '';  end
0161     if isfield(p, 'gh_fcn'),    gh_fcn = p.gh_fcn;      else,   gh_fcn = '';    end
0162     if isfield(p, 'xmax'),      xmax = p.xmax;          else,   xmax = [];      end
0163     if isfield(p, 'xmin'),      xmin = p.xmin;          else,   xmin = [];      end
0164     if isfield(p, 'u'),         u = p.u;                else,   u = [];         end
0165     if isfield(p, 'l'),         l = p.l;                else,   l = [];         end
0166     if isfield(p, 'A'),         A = p.A;                else,   A=sparse(0,nx); end
0167 else                                    %% individual args
0168     nx = size(x0, 1);       %% number of optimization variables
0169     if nargin < 10
0170         opt = [];
0171         if nargin < 9
0172             hess_fcn = '';
0173             if nargin < 8
0174                 gh_fcn = '';
0175                 if nargin < 7
0176                     xmax = [];
0177                     if nargin < 6
0178                         xmin = [];
0179                         if nargin < 5
0180                             u = [];
0181                             if nargin < 4
0182                                 l = [];
0183                                 A = sparse(0,nx);
0184                             end
0185                         end
0186                     end
0187                 end
0188             end
0189         end
0190     end
0191 end
0192 
0193 %% set default argument values if missing
0194 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0195                                  (isempty(u) || all(u == Inf)))
0196     A = sparse(0,nx);           %% no limits => no linear constraints
0197 end
0198 nA = size(A, 1);                %% number of original linear constraints
0199 if isempty(u)                   %% By default, linear inequalities are ...
0200     u = Inf(nA, 1);             %% ... unbounded above and ...
0201 end
0202 if isempty(l)
0203     l = -Inf(nA, 1);            %% ... unbounded below.
0204 end
0205 if isempty(xmin)                %% By default, optimization variables are ...
0206     xmin = -Inf(nx, 1);         %% ... unbounded below and ...
0207 end
0208 if isempty(xmax)
0209     xmax = Inf(nx, 1);          %% ... unbounded above.
0210 end
0211 if isempty(gh_fcn)
0212     nonlinear = false;          %% no nonlinear constraints present
0213 else
0214     nonlinear = true;           %% we have some nonlinear constraints
0215 end
0216 
0217 %% default options
0218 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0219     verbose = opt.verbose;
0220 else
0221     verbose = 0;
0222 end
0223 
0224 %% make sure A is sparse
0225 if ~issparse(A)
0226     A = sparse(A);
0227 end
0228 
0229 %% replace equality variable bounds with an equality constraint
0230 %% (since IPOPT does not return shadow prices on variables that it eliminates)
0231 kk = find(xmin == xmax);
0232 nk = length(kk);
0233 if nk
0234     A = [ A; sparse((1:nk)', kk, 1, nk, nx) ];
0235     l = [ l; xmin(kk) ];
0236     u = [ u; xmax(kk) ];
0237     xmin(kk) = -Inf;
0238     xmax(kk) = Inf;
0239     nA = size(A, 1);            %% updated number of linear constraints
0240 end
0241 
0242 %%-----  set up problem  -----
0243 %% build Jacobian and Hessian structure
0244 randx = rand(size(x0));
0245 nonz = 1e-20;
0246 if nonlinear
0247     [h, g, dhs, dgs] = gh_fcn(randx);
0248     dhs(dhs ~= 0) = nonz;   %% set non-zero entries to tiny value (for adding later)
0249     dgs(dgs ~= 0) = nonz;   %% set non-zero entries to tiny value (for adding later)
0250     Js = [dgs'; dhs'; A];
0251 else
0252     g = []; h = [];
0253     dhs = sparse(0, nx); dgs = dhs; Js = A;
0254 end
0255 neq = length(g);
0256 niq = length(h);
0257 if isempty(hess_fcn)
0258     [f_, df_, Hs] = f_fcn(randx);   %% cost
0259 else
0260     lam = struct('eqnonlin', rand(neq, 1), 'ineqnonlin', rand(niq, 1));
0261     Hs = hess_fcn(randx, lam, 1);
0262 end
0263 if ~issparse(Hs), Hs = sparse(Hs); end      %% convert to sparse if necessary
0264 Hs(Hs ~= 0) = nonz;     %% set non-zero entries to tiny value (for adding later)
0265 
0266 %% set options struct for IPOPT
0267 if ~isempty(opt) && isfield(opt, 'ipopt_opt') && ~isempty(opt.ipopt_opt)
0268     options.ipopt = ipopt_options(opt.ipopt_opt);
0269 else
0270     options.ipopt = ipopt_options;
0271 end
0272 if verbose
0273     options.ipopt.print_level = min(12, verbose*2+1);
0274 else
0275     options.ipopt.print_level = 0;
0276 end
0277 
0278 %% extra data to pass to functions
0279 options.auxdata = struct( ...
0280     'f_fcn',    f_fcn, ...
0281     'gh_fcn',   gh_fcn, ...
0282     'hess_fcn', hess_fcn, ...
0283     'A',        A, ...
0284     'nx',       nx, ...
0285     'nA',       nA, ...
0286     'neqnln',   neq, ...
0287     'niqnln',   niq, ...
0288     'dgs',      dgs, ...
0289     'dhs',      dhs, ...
0290     'Hs',       Hs    );
0291 
0292 %% define variable and constraint bounds
0293 options.lb = xmin;
0294 options.ub = xmax;
0295 options.cl = [zeros(neq, 1);  -Inf(niq, 1); l];
0296 options.cu = [zeros(neq, 1); zeros(niq, 1); u];
0297 
0298 %% assign function handles
0299 funcs.objective         = @objective;
0300 funcs.gradient          = @gradient;
0301 funcs.constraints       = @constraints;
0302 funcs.jacobian          = @jacobian;
0303 funcs.hessian           = @hessian;
0304 funcs.jacobianstructure = @(d) Js;
0305 funcs.hessianstructure  = @(d) tril(Hs);
0306 
0307 %%-----  run solver  -----
0308 %% run the optimization
0309 if have_feature('ipopt_auxdata')
0310     [x, info] = ipopt_auxdata(x0,funcs,options);
0311 else
0312     [x, info] = ipopt(x0,funcs,options);
0313 end
0314 
0315 if info.status == 0 || info.status == 1
0316     eflag = 1;
0317 else
0318     eflag = 0;
0319 end
0320 output = struct('status', info.status);
0321 if isfield(info, 'iter')
0322     output.iterations = info.iter;
0323 else
0324     output.iterations = [];
0325 end
0326 if isfield(info, 'cpu')
0327     output.cpu = info.cpu;
0328 end
0329 if isfield(info, 'eval')
0330     output.eval = info.eval;
0331 end
0332 f = f_fcn(x);
0333 
0334 %% check for empty results (in case optimization failed)
0335 if isempty(info.lambda)
0336     lam = NaN(nA, 1);
0337 else
0338     lam = info.lambda;
0339 end
0340 if isempty(info.zl) && ~isempty(xmin)
0341     zl = NaN(nx, 1);
0342 else
0343     zl = info.zl;
0344 end
0345 if isempty(info.zu) && ~isempty(xmax)
0346     zu = NaN(nx, 1);
0347 else
0348     zu = info.zu;
0349 end
0350 
0351 %% extract shadow prices for equality var bounds converted to eq constraints
0352 %% (since IPOPT does not return shadow prices on variables that it eliminates)
0353 if nk
0354     offset = neq + niq + nA - nk;
0355     lam_tmp = lam(offset+(1:nk));
0356     kl = find(lam_tmp < 0);         %% lower bound binding
0357     ku = find(lam_tmp > 0);         %% upper bound binding
0358     zl(kk(kl)) = -lam_tmp(kl);
0359     zu(kk(ku)) =  lam_tmp(ku);
0360     lam(offset+(1:nk)) = [];        %% remove these shadow prices
0361     nA = nA - nk;                   %% reduce dimension accordingly
0362 end
0363 
0364 %% extract multipliers for linear constraints
0365 lam_lin = lam(neq+niq+(1:nA));      %% lambda for linear constraints
0366 kl = find(lam_lin < 0);             %% lower bound binding
0367 ku = find(lam_lin > 0);             %% upper bound binding
0368 mu_l = zeros(nA, 1);
0369 mu_l(kl) = -lam_lin(kl);
0370 mu_u = zeros(nA, 1);
0371 mu_u(ku) = lam_lin(ku);
0372 
0373 lambda = struct( ...
0374     'lower', zl, ...
0375     'upper', zu, ...
0376     'eqnonlin', lam(1:neq), ...
0377     'ineqnonlin', lam(neq+(1:niq)), ...
0378     'mu_l', mu_l, ...
0379     'mu_u', mu_u );
0380 
0381 
0382 %-----  callback functions  -----
0383 function f = objective(x, d)
0384 f = d.f_fcn(x);
0385 
0386 function df = gradient(x, d)
0387 [f, df] = d.f_fcn(x);
0388 
0389 function c = constraints(x, d)
0390 if isempty(d.gh_fcn)
0391     c = d.A*x;
0392 else
0393     [h, g] = d.gh_fcn(x);
0394     c = [g; h; d.A*x];
0395 end
0396 
0397 function J = jacobian(x, d)
0398 if isempty(d.gh_fcn)
0399     J = d.A;
0400 else
0401     [h, g, dh, dg] = d.gh_fcn(x);
0402     %% add sparse structure (with tiny values) to current matrices to
0403     %% ensure that sparsity structure matches that supplied
0404     J = [(dg + d.dgs)'; (dh + d.dhs)'; d.A];
0405 end
0406 
0407 function H = hessian(x, sigma, lambda, d)
0408 if isempty(d.hess_fcn)
0409     [f_, df_, d2f] = d.f_fcn(x);    %% cost
0410     if ~issparse(d2f), d2f = sparse(d2f); end   %% convert to sparse if necessary
0411     H = tril(d2f * sigma);
0412 else
0413     lam.eqnonlin   = lambda(1:d.neqnln);
0414     lam.ineqnonlin = lambda(d.neqnln+(1:d.niqnln));
0415     H = d.hess_fcn(x, lam, sigma);
0416     if ~issparse(H), H = sparse(H); end     %% convert to sparse if necessary
0417     %% add sparse structure (with tiny values) to current matrices to
0418     %% ensure that sparsity structure matches that supplied
0419     H = tril(H + d.Hs);
0420 end

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