Home > matpower5.0 > mips.m

mips

PURPOSE ^

MIPS MATLAB Interior Point Solver.

SYNOPSIS ^

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

DESCRIPTION ^

MIPS  MATLAB Interior Point Solver.
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
       MIPS(F_FCN, X0, A, L, U, XMIN, XMAX, GH_FCN, HESS_FCN, OPT)
   Primal-dual interior point method for NLP (nonlinear programming).
   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)
       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
           feastol (1e-6) - termination tolerance for feasibility
               condition
           gradtol (1e-6) - termination tolerance for gradient
               condition
           comptol (1e-6) - termination tolerance for complementarity
               condition
           costtol (1e-6) - termination tolerance for cost condition
           max_it (150) - maximum number of iterations
           step_control (0) - set to 1 to enable step-size control
           sc.red_it (20) - maximum number of step-size reductions if
               step-control is on
           cost_mult (1) - cost multiplier used to scale the objective
               function for improved conditioning. Note: This value is
               also passed as the 3rd argument to the Hessian evaluation
               function so that it can appropriately scale the
               objective function term in the Hessian of the
               Lagrangian.
           xi (0.99995) - constant used in alpha updates*
           sigma (0.1) - centering parameter*
           z0 (1) - used to initialize slack variables*
           alpha_min (1e-8) - returns EXITFLAG = -1 if either alpha
               parameter becomes smaller than this value*
           rho_min (0.95) - lower bound on rho_t*
           rho_max (1.05) - upper bound on rho_t*
           mu_threshold (1e-5) - KT multipliers smaller than this value
               for non-binding constraints are forced to zero
           max_stepsize (1e10) - returns EXITFLAG = -1 if the 2-norm of
               the reduced Newton step exceeds this value*
               * see the corresponding Appendix in the manual for details
       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 = first order optimality conditions satisfied
           0 = maximum number of iterations reached
           -1 = numerically failed
       OUTPUT : output struct with 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:
           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] = ...
           mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);

       x = mips(f_fcn, x0);
       x = mips(f_fcn, x0, A, l);
       x = mips(f_fcn, x0, A, l, u);
       x = mips(f_fcn, x0, A, l, u, xmin);
       x = mips(f_fcn, x0, A, l, u, xmin, xmax);
       x = mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn);
       x = mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn);
       x = mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);
       x = mips(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 = mips(...);
       [x, f] = mips(...);
       [x, f, exitflag] = mips(...);
       [x, f, exitflag, output] = mips(...);
       [x, f, exitflag, output, lambda] = mips(...);

   Example: (problem from http://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] = mips(problem);

   Ported by Ray Zimmerman from C code written by H. Wang for his
   PhD dissertation:
     "On the Computation and Application of Multi-period
     Security-Constrained Optimal Power Flow for Real-time
     Electricity Market Operations", Cornell University, May 2007.

   See also:
     H. Wang, C. E. Murillo-S‡nchez, R. D. Zimmerman, R. J. Thomas,
     "On Computational Issues of Market-Based Optimal Power Flow",
     IEEE Transactions on Power Systems, Vol. 22, No. 3, Aug. 2007,
     pp. 1185-1193.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x, f, eflag, output, lambda] = mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt)
0002 %MIPS  MATLAB Interior Point Solver.
0003 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
0004 %       MIPS(F_FCN, X0, A, L, U, XMIN, XMAX, GH_FCN, HESS_FCN, OPT)
0005 %   Primal-dual interior point method for NLP (nonlinear programming).
0006 %   Minimize a function F(X) beginning from a starting point X0, subject to
0007 %   optional linear and nonlinear constraints and variable bounds.
0008 %
0009 %       min F(X)
0010 %        X
0011 %
0012 %   subject to
0013 %
0014 %       G(X) = 0            (nonlinear equalities)
0015 %       H(X) <= 0           (nonlinear inequalities)
0016 %       L <= A*X <= U       (linear constraints)
0017 %       XMIN <= X <= XMAX   (variable bounds)
0018 %
0019 %   Inputs (all optional except F_FCN and X0):
0020 %       F_FCN : handle to function that evaluates the objective function,
0021 %           its gradients and Hessian for a given value of X. If there
0022 %           are nonlinear constraints, the Hessian information is
0023 %           provided by the HESS_FCN function passed in the 9th argument
0024 %           and is not required here. Calling syntax for this function:
0025 %               [F, DF, D2F] = F_FCN(X)
0026 %       X0 : starting value of optimization vector X
0027 %       A, L, U : define the optional linear constraints. Default
0028 %           values for the elements of L and U are -Inf and Inf,
0029 %           respectively.
0030 %       XMIN, XMAX : optional lower and upper bounds on the
0031 %           X variables, defaults are -Inf and Inf, respectively.
0032 %       GH_FCN : handle to function that evaluates the optional
0033 %           nonlinear constraints and their gradients for a given
0034 %           value of X. Calling syntax for this function is:
0035 %               [H, G, DH, DG] = GH_FCN(X)
0036 %       HESS_FCN : handle to function that computes the Hessian of the
0037 %           Lagrangian for given values of X, lambda and mu, where
0038 %           lambda and mu are the multipliers on the equality and
0039 %           inequality constraints, g and h, respectively. The calling
0040 %           syntax for this function is:
0041 %               LXX = HESS_FCN(X, LAM)
0042 %           where lambda = LAM.eqnonlin and mu = LAM.ineqnonlin.
0043 %       OPT : optional options structure with the following fields,
0044 %           all of which are also optional (default values shown in
0045 %           parentheses)
0046 %           verbose (0) - controls level of progress output displayed
0047 %           feastol (1e-6) - termination tolerance for feasibility
0048 %               condition
0049 %           gradtol (1e-6) - termination tolerance for gradient
0050 %               condition
0051 %           comptol (1e-6) - termination tolerance for complementarity
0052 %               condition
0053 %           costtol (1e-6) - termination tolerance for cost condition
0054 %           max_it (150) - maximum number of iterations
0055 %           step_control (0) - set to 1 to enable step-size control
0056 %           sc.red_it (20) - maximum number of step-size reductions if
0057 %               step-control is on
0058 %           cost_mult (1) - cost multiplier used to scale the objective
0059 %               function for improved conditioning. Note: This value is
0060 %               also passed as the 3rd argument to the Hessian evaluation
0061 %               function so that it can appropriately scale the
0062 %               objective function term in the Hessian of the
0063 %               Lagrangian.
0064 %           xi (0.99995) - constant used in alpha updates*
0065 %           sigma (0.1) - centering parameter*
0066 %           z0 (1) - used to initialize slack variables*
0067 %           alpha_min (1e-8) - returns EXITFLAG = -1 if either alpha
0068 %               parameter becomes smaller than this value*
0069 %           rho_min (0.95) - lower bound on rho_t*
0070 %           rho_max (1.05) - upper bound on rho_t*
0071 %           mu_threshold (1e-5) - KT multipliers smaller than this value
0072 %               for non-binding constraints are forced to zero
0073 %           max_stepsize (1e10) - returns EXITFLAG = -1 if the 2-norm of
0074 %               the reduced Newton step exceeds this value*
0075 %               * see the corresponding Appendix in the manual for details
0076 %       PROBLEM : The inputs can alternatively be supplied in a single
0077 %           PROBLEM struct with fields corresponding to the input arguments
0078 %           described above: f_fcn, x0, A, l, u, xmin, xmax,
0079 %                            gh_fcn, hess_fcn, opt
0080 %
0081 %   Outputs:
0082 %       X : solution vector
0083 %       F : final objective function value
0084 %       EXITFLAG : exit flag
0085 %           1 = first order optimality conditions satisfied
0086 %           0 = maximum number of iterations reached
0087 %           -1 = numerically failed
0088 %       OUTPUT : output struct with fields:
0089 %           iterations - number of iterations performed
0090 %           hist - struct array with trajectories of the following:
0091 %                   feascond, gradcond, compcond, costcond, gamma,
0092 %                   stepsize, obj, alphap, alphad
0093 %           message - exit message
0094 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0095 %           multipliers on the constraints, with fields:
0096 %           eqnonlin - nonlinear equality constraints
0097 %           ineqnonlin - nonlinear inequality constraints
0098 %           mu_l - lower (left-hand) limit on linear constraints
0099 %           mu_u - upper (right-hand) limit on linear constraints
0100 %           lower - lower bound on optimization variables
0101 %           upper - upper bound on optimization variables
0102 %
0103 %   Note the calling syntax is almost identical to that of FMINCON
0104 %   from MathWorks' Optimization Toolbox. The main difference is that
0105 %   the linear constraints are specified with A, L, U instead of
0106 %   A, B, Aeq, Beq. The functions for evaluating the objective
0107 %   function, constraints and Hessian are identical.
0108 %
0109 %   Calling syntax options:
0110 %       [x, f, exitflag, output, lambda] = ...
0111 %           mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);
0112 %
0113 %       x = mips(f_fcn, x0);
0114 %       x = mips(f_fcn, x0, A, l);
0115 %       x = mips(f_fcn, x0, A, l, u);
0116 %       x = mips(f_fcn, x0, A, l, u, xmin);
0117 %       x = mips(f_fcn, x0, A, l, u, xmin, xmax);
0118 %       x = mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn);
0119 %       x = mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn);
0120 %       x = mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);
0121 %       x = mips(problem);
0122 %               where problem is a struct with fields:
0123 %                   f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt
0124 %                   all fields except 'f_fcn' and 'x0' are optional
0125 %       x = mips(...);
0126 %       [x, f] = mips(...);
0127 %       [x, f, exitflag] = mips(...);
0128 %       [x, f, exitflag, output] = mips(...);
0129 %       [x, f, exitflag, output, lambda] = mips(...);
0130 %
0131 %   Example: (problem from http://en.wikipedia.org/wiki/Nonlinear_programming)
0132 %       function [f, df, d2f] = f2(x)
0133 %       f = -x(1)*x(2) - x(2)*x(3);
0134 %       if nargout > 1           %% gradient is required
0135 %           df = -[x(2); x(1)+x(3); x(2)];
0136 %           if nargout > 2       %% Hessian is required
0137 %               d2f = -[0 1 0; 1 0 1; 0 1 0];   %% actually not used since
0138 %           end                                 %% 'hess_fcn' is provided
0139 %       end
0140 %
0141 %       function [h, g, dh, dg] = gh2(x)
0142 %       h = [ 1 -1 1; 1 1 1] * x.^2 + [-2; -10];
0143 %       dh = 2 * [x(1) x(1); -x(2) x(2); x(3) x(3)];
0144 %       g = []; dg = [];
0145 %
0146 %       function Lxx = hess2(x, lam, cost_mult)
0147 %       if nargin < 3, cost_mult = 1; end
0148 %       mu = lam.ineqnonlin;
0149 %       Lxx = cost_mult * [0 -1 0; -1 0 -1; 0 -1 0] + ...
0150 %               [2*[1 1]*mu 0 0; 0 2*[-1 1]*mu 0; 0 0 2*[1 1]*mu];
0151 %
0152 %       problem = struct( ...
0153 %           'f_fcn',    @(x)f2(x), ...
0154 %           'gh_fcn',   @(x)gh2(x), ...
0155 %           'hess_fcn', @(x, lam, cost_mult)hess2(x, lam, cost_mult), ...
0156 %           'x0',       [1; 1; 0], ...
0157 %           'opt',      struct('verbose', 2) ...
0158 %       );
0159 %       [x, f, exitflag, output, lambda] = mips(problem);
0160 %
0161 %   Ported by Ray Zimmerman from C code written by H. Wang for his
0162 %   PhD dissertation:
0163 %     "On the Computation and Application of Multi-period
0164 %     Security-Constrained Optimal Power Flow for Real-time
0165 %     Electricity Market Operations", Cornell University, May 2007.
0166 %
0167 %   See also:
0168 %     H. Wang, C. E. Murillo-S‡nchez, R. D. Zimmerman, R. J. Thomas,
0169 %     "On Computational Issues of Market-Based Optimal Power Flow",
0170 %     IEEE Transactions on Power Systems, Vol. 22, No. 3, Aug. 2007,
0171 %     pp. 1185-1193.
0172 
0173 %   MIPS
0174 %   $Id: mips.m 2437 2014-12-04 19:13:54Z ray $
0175 %   by Ray Zimmerman, PSERC Cornell
0176 %   Copyright (c) 2009-2013 by Power System Engineering Research Center (PSERC)
0177 %
0178 %   This file is part of MIPS.
0179 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0180 %
0181 %   MIPS is free software: you can redistribute it and/or modify
0182 %   it under the terms of the GNU General Public License as published
0183 %   by the Free Software Foundation, either version 3 of the License,
0184 %   or (at your option) any later version.
0185 %
0186 %   MIPS is distributed in the hope that it will be useful,
0187 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0188 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0189 %   GNU General Public License for more details.
0190 %
0191 %   You should have received a copy of the GNU General Public License
0192 %   along with MIPS. If not, see <http://www.gnu.org/licenses/>.
0193 %
0194 %   Additional permission under GNU GPL version 3 section 7
0195 %
0196 %   If you modify MIPS, or any covered work, to interface with
0197 %   other modules (such as MATLAB code and MEX-files) available in a
0198 %   MATLAB(R) or comparable environment containing parts covered
0199 %   under other licensing terms, the licensors of MIPS grant
0200 %   you additional permission to convey the resulting work.
0201 
0202 %%----- input argument handling  -----
0203 %% gather inputs
0204 if nargin == 1 && isstruct(f_fcn)       %% problem struct
0205     p = f_fcn;
0206     f_fcn = p.f_fcn;
0207     x0 = p.x0;
0208     nx = size(x0, 1);       %% number of optimization variables
0209     if isfield(p, 'opt'),       opt = p.opt;            else,   opt = [];       end
0210     if isfield(p, 'hess_fcn'),  hess_fcn = p.hess_fcn;  else,   hess_fcn = '';  end
0211     if isfield(p, 'gh_fcn'),    gh_fcn = p.gh_fcn;      else,   gh_fcn = '';    end
0212     if isfield(p, 'xmax'),      xmax = p.xmax;          else,   xmax = [];      end
0213     if isfield(p, 'xmin'),      xmin = p.xmin;          else,   xmin = [];      end
0214     if isfield(p, 'u'),         u = p.u;                else,   u = [];         end
0215     if isfield(p, 'l'),         l = p.l;                else,   l = [];         end
0216     if isfield(p, 'A'),         A = p.A;                else,   A=sparse(0,nx); end
0217 else                                    %% individual args
0218     nx = size(x0, 1);       %% number of optimization variables
0219     if nargin < 10
0220         opt = [];
0221         if nargin < 9
0222             hess_fcn = '';
0223             if nargin < 8
0224                 gh_fcn = '';
0225                 if nargin < 7
0226                     xmax = [];
0227                     if nargin < 6
0228                         xmin = [];
0229                         if nargin < 5
0230                             u = [];
0231                             if nargin < 4
0232                                 l = [];
0233                                 A = sparse(0,nx);
0234                             end
0235                         end
0236                     end
0237                 end
0238             end
0239         end
0240     end
0241 end
0242 %% set default argument values if missing
0243 if  ~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0244                    (isempty(u) || all(u == Inf))
0245     A = sparse(0,nx);           %% no limits => no linear constraints
0246 end
0247 nA = size(A, 1);                %% number of original linear constraints
0248 if isempty(u)                   %% By default, linear inequalities are ...
0249     u = Inf * ones(nA, 1);      %% ... unbounded above and ...
0250 end
0251 if isempty(l)
0252     l = -Inf * ones(nA, 1);     %% ... unbounded below.
0253 end
0254 if isempty(xmin)                %% By default, optimization variables are ...
0255     xmin = -Inf * ones(nx, 1);  %% ... unbounded below and ...
0256 end
0257 if isempty(xmax)
0258     xmax = Inf * ones(nx, 1);   %% ... unbounded above.
0259 end
0260 if isempty(gh_fcn)
0261     nonlinear = false;          %% no nonlinear constraints present
0262     gn = []; hn = [];
0263 else
0264     nonlinear = true;           %% we have some nonlinear constraints
0265 end
0266 
0267 %% default options
0268 if isempty(opt)
0269     opt = struct('verbose', 0);
0270 end
0271 if ~isfield(opt, 'feastol') || isempty(opt.feastol)
0272     opt.feastol = 1e-6;
0273 end
0274 if ~isfield(opt, 'gradtol') || isempty(opt.gradtol)
0275     opt.gradtol = 1e-6;
0276 end
0277 if ~isfield(opt, 'comptol') || isempty(opt.comptol)
0278     opt.comptol = 1e-6;
0279 end
0280 if ~isfield(opt, 'costtol') || isempty(opt.costtol)
0281     opt.costtol = 1e-6;
0282 end
0283 if ~isfield(opt, 'max_it') || isempty(opt.max_it)
0284     opt.max_it = 150;
0285 end
0286 if ~isfield(opt, 'sc') || ~isfield(opt.sc, 'red_it') || isempty(opt.sc.red_it)
0287     opt.sc.red_it = 20;
0288 end
0289 if ~isfield(opt, 'step_control') || isempty(opt.step_control)
0290     opt.step_control = 0;
0291 end
0292 if ~isfield(opt, 'cost_mult') || isempty(opt.cost_mult)
0293     opt.cost_mult = 1;
0294 end
0295 if ~isfield(opt, 'verbose') || isempty(opt.verbose)
0296     opt.verbose = 0;
0297 end
0298 %% algorithm constants
0299 if ~isfield(opt, 'xi') || isempty(opt.xi)
0300     opt.xi = 0.99995;           %% OPT_IPM_PHI
0301 end
0302 if ~isfield(opt, 'sigma') || isempty(opt.sigma)
0303     opt.sigma = 0.1;            %% OPT_IPM_SIGMA, centering parameter
0304 end
0305 if ~isfield(opt, 'z0') || isempty(opt.z0)
0306     opt.z0 = 1;                 %% OPT_IPM_INIT_SLACK
0307 end
0308 if ~isfield(opt, 'alpha_min') || isempty(opt.alpha_min)
0309     opt.alpha_min = 1e-8;       %% OPT_AP_AD_MIN
0310 end
0311 if ~isfield(opt, 'rho_min') || isempty(opt.rho_min)
0312     opt.rho_min = 0.95;         %% OPT_IPM_QUAD_LOWTHRESH
0313 end
0314 if ~isfield(opt, 'rho_max') || isempty(opt.rho_max)
0315     opt.rho_max = 1.05;         %% OPT_IPM_QUAD_HIGHTHRESH
0316 end
0317 if ~isfield(opt, 'mu_threshold') || isempty(opt.mu_threshold)
0318     opt.mu_threshold = 1e-5;    %% SCOPF_MULTIPLIERS_FILTER_THRESH
0319 end
0320 if ~isfield(opt, 'max_stepsize') || isempty(opt.max_stepsize)
0321     opt.max_stepsize = 1e10;
0322 end
0323 
0324 %% initialize history
0325 hist(opt.max_it+1) = struct('feascond', 0, 'gradcond', 0, 'compcond', 0, ...
0326     'costcond', 0, 'gamma', 0, 'stepsize', 0, 'obj', 0, ...
0327     'alphap', 0, 'alphad', 0);
0328 
0329 %%-----  set up problem  -----
0330 %% constants
0331 [xi, sigma, z0, alpha_min, rho_min, rho_max, mu_threshold, max_stepsize] = ...
0332     deal(opt.xi, opt.sigma, opt.z0, opt.alpha_min, ...
0333         opt.rho_min, opt.rho_max, opt.mu_threshold, opt.max_stepsize);
0334 if xi >= 1 || xi < 0.5
0335     error('mips: opt.xi (%g) must be a number slightly less than 1', opt.xi);
0336 end
0337 if sigma > 1 || sigma <= 0
0338     error('mips: opt.sigma (%g) must be a number between 0 and 1', opt.sigma);
0339 end
0340 
0341 %% initialize
0342 i = 0;                      %% iteration counter
0343 converged = 0;              %% flag
0344 eflag = 0;                  %% exit flag
0345 
0346 %% add var limits to linear constraints
0347 AA = [speye(nx); A];
0348 ll = [xmin; l];
0349 uu = [xmax; u];
0350 
0351 %% split up linear constraints
0352 ieq = find( abs(uu-ll) <= eps );            %% equality
0353 igt = find( uu >=  1e10 & ll > -1e10 );     %% greater than, unbounded above
0354 ilt = find( ll <= -1e10 & uu <  1e10 );     %% less than, unbounded below
0355 ibx = find( (abs(uu-ll) > eps) & (uu < 1e10) & (ll > -1e10) );
0356 Ae = AA(ieq, :);
0357 be = uu(ieq, 1);
0358 Ai  = [ AA(ilt, :); -AA(igt, :); AA(ibx, :); -AA(ibx, :) ];
0359 bi  = [ uu(ilt, 1); -ll(igt, 1); uu(ibx, 1); -ll(ibx, 1) ];
0360 
0361 %% evaluate cost f(x0) and constraints g(x0), h(x0)
0362 x = x0;
0363 [f, df] = f_fcn(x);             %% cost
0364 f = f * opt.cost_mult;
0365 df = df * opt.cost_mult;
0366 if nonlinear
0367     [hn, gn, dhn, dgn] = gh_fcn(x); %% nonlinear constraints
0368     h = [hn; Ai * x - bi];          %% inequality constraints
0369     g = [gn; Ae * x - be];          %% equality constraints
0370     dh = [dhn Ai'];                 %% 1st derivative of inequalities
0371     dg = [dgn Ae'];                 %% 1st derivative of equalities
0372 else
0373     h = Ai * x - bi;                %% inequality constraints
0374     g = Ae * x - be;                %% equality constraints
0375     dh = Ai';                       %% 1st derivative of inequalities
0376     dg = Ae';                       %% 1st derivative of equalities
0377 end
0378 
0379 %% grab some dimensions
0380 neq = size(g, 1);           %% number of equality constraints
0381 niq = size(h, 1);           %% number of inequality constraints
0382 neqnln = size(gn, 1);       %% number of nonlinear equality constraints
0383 niqnln = size(hn, 1);       %% number of nonlinear inequality constraints
0384 nlt = length(ilt);          %% number of upper bounded linear inequalities
0385 ngt = length(igt);          %% number of lower bounded linear inequalities
0386 nbx = length(ibx);          %% number of doubly bounded linear inequalities
0387 
0388 %% initialize gamma, lam, mu, z, e
0389 gamma = 1;                  %% barrier coefficient, r in Harry's code
0390 lam = zeros(neq, 1);
0391 z   = z0 * ones(niq, 1);
0392 mu  = z;
0393 k = find(h < -z0);
0394 z(k) = -h(k);
0395 k = find(gamma / z > z0);   %% (seems k is always empty if gamma = z0 = 1)
0396 if ~isempty(k)
0397     mu(k) = gamma / z(k);
0398 end
0399 e = ones(niq, 1);
0400 
0401 %% check tolerance
0402 f0 = f;
0403 if opt.step_control
0404     L = f + lam' * g + mu' * (h+z) - gamma * sum(log(z));
0405 end
0406 Lx = df + dg * lam + dh * mu;
0407 if isempty(h)
0408     maxh = zeros(1,0);
0409 else
0410     maxh = max(h);
0411 end
0412 feascond = max([norm(g, Inf), maxh]) / (1 + max([ norm(x, Inf), norm(z, Inf) ]));
0413 gradcond = norm(Lx, Inf) / (1 + max([ norm(lam, Inf), norm(mu, Inf) ]));
0414 compcond = (z' * mu) / (1 + norm(x, Inf));
0415 costcond = abs(f - f0) / (1 + abs(f0));
0416 %% save history
0417 hist(i+1) = struct('feascond', feascond, 'gradcond', gradcond, ...
0418     'compcond', compcond, 'costcond', costcond, 'gamma', gamma, ...
0419     'stepsize', 0, 'obj', f/opt.cost_mult, 'alphap', 0, 'alphad', 0);
0420 if opt.verbose
0421     if opt.step_control, s = '-sc'; else, s = ''; end
0422     v = mipsver('all');
0423     fprintf('MATLAB Interior Point Solver -- MIPS%s, Version %s, %s', ...
0424         s, v.Version, v.Date);
0425     if opt.verbose > 1
0426         fprintf('\n it    objective   step size   feascond     gradcond     compcond     costcond  ');
0427         fprintf('\n----  ------------ --------- ------------ ------------ ------------ ------------');
0428         fprintf('\n%3d  %12.8g %10s %12g %12g %12g %12g', ...
0429             i, f/opt.cost_mult, '', feascond, gradcond, compcond, costcond);
0430     end
0431 end
0432 if feascond < opt.feastol && gradcond < opt.gradtol && ...
0433                 compcond < opt.comptol && costcond < opt.costtol
0434     converged = 1;
0435     if opt.verbose
0436         fprintf('\nConverged!\n');
0437     end
0438 end
0439 
0440 %%-----  do Newton iterations  -----
0441 while (~converged && i < opt.max_it)
0442     %% update iteration counter
0443     i = i + 1;
0444 
0445     %% compute update step
0446     lambda = struct('eqnonlin', lam(1:neqnln), 'ineqnonlin', mu(1:niqnln));
0447     if nonlinear
0448         if isempty(hess_fcn)
0449             fprintf('mips: Hessian evaluation via finite differences not yet implemented.\n       Please provide your own hessian evaluation function.');
0450         end
0451         Lxx = hess_fcn(x, lambda, opt.cost_mult);
0452     else
0453         [f_, df_, d2f] = f_fcn(x);      %% cost
0454         Lxx = d2f * opt.cost_mult;
0455     end
0456     zinvdiag = sparse(1:niq, 1:niq, 1 ./ z, niq, niq);
0457     mudiag = sparse(1:niq, 1:niq, mu, niq, niq);
0458     dh_zinv = dh * zinvdiag;
0459     M = Lxx + dh_zinv * mudiag * dh';
0460     N = Lx + dh_zinv * (mudiag * h + gamma * e);
0461     dxdlam = [M dg; dg' sparse(neq, neq)] \ [-N; -g];
0462 %     AAA = [
0463 %         M  dg;
0464 %         dg'  sparse(neq, neq)
0465 %     ];
0466 %     rc = 1/condest(AAA);
0467 %     if rc < 1e-22
0468 %         fprintf('my RCOND = %g\n', rc);
0469 %         n = size(AAA, 1);
0470 %         AAA = AAA + 1e-3 * speye(n,n);
0471 %     end
0472 %     bbb = [-N; -g];
0473 %     dxdlam = AAA \ bbb;
0474     if any(isnan(dxdlam)) || norm(dxdlam) > max_stepsize
0475         if opt.verbose
0476             fprintf('\nNumerically Failed\n');
0477         end
0478         eflag = -1;
0479         break;
0480     end
0481     dx = dxdlam(1:nx, 1);
0482     dlam = dxdlam(nx+(1:neq), 1);
0483     dz = -h - z - dh' * dx;
0484     dmu = -mu + zinvdiag *(gamma*e - mudiag * dz);
0485 
0486     %% optional step-size control
0487     sc = 0;
0488     if opt.step_control
0489         x1 = x + dx;
0490 
0491         %% evaluate cost, constraints, derivatives at x1
0492         [f1, df1] = f_fcn(x1);          %% cost
0493         f1 = f1 * opt.cost_mult;
0494         df1 = df1 * opt.cost_mult;
0495         if nonlinear
0496             [hn1, gn1, dhn1, dgn1] = gh_fcn(x1); %% nonlinear constraints
0497             h1 = [hn1; Ai * x1 - bi];       %% inequality constraints
0498             g1 = [gn1; Ae * x1 - be];       %% equality constraints
0499             dh1 = [dhn1 Ai'];               %% 1st derivative of inequalities
0500             dg1 = [dgn1 Ae'];               %% 1st derivative of equalities
0501         else
0502             h1 = Ai * x1 - bi;              %% inequality constraints
0503             g1 = Ae * x1 - be;              %% equality constraints
0504             dh1 = dh;                       %% 1st derivative of inequalities
0505             dg1 = dg;                       %% 1st derivative of equalities
0506         end
0507 
0508         %% check tolerance
0509         Lx1 = df1 + dg1 * lam + dh1 * mu;
0510         if isempty(h1)
0511             maxh1 = zeros(1,0);
0512         else
0513             maxh1 = max(h1);
0514         end
0515         feascond1 = max([norm(g1, Inf), maxh1]) / (1 + max([ norm(x1, Inf), norm(z, Inf) ]));
0516         gradcond1 = norm(Lx1, Inf) / (1 + max([ norm(lam, Inf), norm(mu, Inf) ]));
0517 
0518         if feascond1 > feascond && gradcond1 > gradcond
0519             sc = 1;
0520         end
0521     end
0522     if sc
0523         alpha = 1;
0524         for j = 1:opt.sc.red_it
0525             dx1 = alpha * dx;
0526             x1 = x + dx1;
0527             f1 = f_fcn(x1);                 %% cost
0528             f1 = f1 * opt.cost_mult;
0529             if nonlinear
0530                 [hn1, gn1] = gh_fcn(x1);    %% nonlinear constraints
0531                 h1 = [hn1; Ai * x1 - bi];   %% inequality constraints
0532                 g1 = [gn1; Ae * x1 - be];   %% equality constraints
0533             else
0534                 h1 = Ai * x1 - bi;          %% inequality constraints
0535                 g1 = Ae * x1 - be;          %% equality constraints
0536             end
0537             L1 = f1 + lam' * g1 + mu' * (h1+z) - gamma * sum(log(z));
0538             if opt.verbose > 2
0539                 fprintf('\n   %3d            %10g', -j, norm(dx1));
0540             end
0541             rho = (L1 - L) / (Lx' * dx1 + 0.5 * dx1' * Lxx * dx1);
0542             if rho > rho_min && rho < rho_max
0543                 break;
0544             else
0545                 alpha = alpha / 2;
0546             end
0547         end
0548         dx = alpha * dx;
0549         dz = alpha * dz;
0550         dlam = alpha * dlam;
0551         dmu = alpha * dmu;
0552     end
0553 
0554     %% do the update
0555     k = find(dz < 0);
0556     if isempty(k)
0557         alphap = 1;
0558     else
0559         alphap = min( [xi * min(z(k) ./ -dz(k)) 1] );
0560     end
0561     k = find(dmu < 0);
0562     if isempty(k)
0563         alphad = 1;
0564     else
0565         alphad = min( [xi * min(mu(k) ./ -dmu(k)) 1] );
0566     end
0567     x = x + alphap * dx;
0568     z = z + alphap * dz;
0569     lam = lam + alphad * dlam;
0570     mu  = mu  + alphad * dmu;
0571     if niq > 0
0572         gamma = sigma * (z' * mu) / niq;
0573     end
0574 
0575     %% evaluate cost, constraints, derivatives
0576     [f, df] = f_fcn(x);                 %% cost
0577     f = f * opt.cost_mult;
0578     df = df * opt.cost_mult;
0579     if nonlinear
0580         [hn, gn, dhn, dgn] = gh_fcn(x); %% nonlinear constraints
0581         h = [hn; Ai * x - bi];          %% inequality constraints
0582         g = [gn; Ae * x - be];          %% equality constraints
0583         dh = [dhn Ai'];                 %% 1st derivative of inequalities
0584         dg = [dgn Ae'];                 %% 1st derivative of equalities
0585     else
0586         h = Ai * x - bi;                %% inequality constraints
0587         g = Ae * x - be;                %% equality constraints
0588         %% 1st derivatives are constant, still dh = Ai', dg = Ae'
0589     end
0590 
0591     %% check tolerance
0592     Lx = df + dg * lam + dh * mu;
0593     if isempty(h)
0594         maxh = zeros(1,0);
0595     else
0596         maxh = max(h);
0597     end
0598     feascond = max([norm(g, Inf), maxh]) / (1 + max([ norm(x, Inf), norm(z, Inf) ]));
0599     gradcond = norm(Lx, Inf) / (1 + max([ norm(lam, Inf), norm(mu, Inf) ]));
0600     compcond = (z' * mu) / (1 + norm(x, Inf));
0601     costcond = abs(f - f0) / (1 + abs(f0));
0602     %% save history
0603     hist(i+1) = struct('feascond', feascond, 'gradcond', gradcond, ...
0604         'compcond', compcond, 'costcond', costcond, 'gamma', gamma, ...
0605         'stepsize', norm(dx), 'obj', f/opt.cost_mult, ...
0606         'alphap', alphap, 'alphad', alphad);
0607 
0608     if opt.verbose > 1
0609         fprintf('\n%3d  %12.8g %10.5g %12g %12g %12g %12g', ...
0610             i, f/opt.cost_mult, norm(dx), feascond, gradcond, compcond, costcond);
0611     end
0612     if feascond < opt.feastol && gradcond < opt.gradtol && ...
0613                     compcond < opt.comptol && costcond < opt.costtol
0614         converged = 1;
0615         if opt.verbose
0616             fprintf('\nConverged!\n');
0617         end
0618     else
0619         if any(isnan(x)) || alphap < alpha_min || alphad < alpha_min || ...
0620                 gamma < eps || gamma > 1/eps
0621             if opt.verbose
0622                 fprintf('\nNumerically Failed\n');
0623             end
0624             eflag = -1;
0625             break;
0626         end
0627         f0 = f;
0628         if opt.step_control
0629             L = f + lam' * g + mu' * (h+z) - gamma * sum(log(z));
0630         end
0631     end
0632 end
0633 
0634 if opt.verbose
0635     if ~converged
0636         fprintf('\nDid not converge in %d iterations.\n', i);
0637     end
0638 end
0639 
0640 %%-----  package up results  -----
0641 hist = hist(1:i+1);
0642 if eflag ~= -1
0643     eflag = converged;
0644 end
0645 output = struct('iterations', i, 'hist', hist, 'message', '');
0646 if eflag == 0
0647     output.message = 'Did not converge';
0648 elseif eflag == 1
0649     output.message = 'Converged';
0650 elseif eflag == -1
0651     output.message = 'Numerically failed';
0652 else
0653     output.message = 'Please hang up and dial again';
0654 end
0655 
0656 %% zero out multipliers on non-binding constraints
0657 mu(h < -opt.feastol & mu < mu_threshold) = 0;
0658 
0659 %% un-scale cost and prices
0660 f   = f   / opt.cost_mult;
0661 lam = lam / opt.cost_mult;
0662 mu  = mu  / opt.cost_mult;
0663 
0664 %% re-package multipliers into struct
0665 lam_lin = lam((neqnln+1):neq);              %% lambda for linear constraints
0666 mu_lin  = mu((niqnln+1):niq);               %% mu for linear constraints
0667 kl = find(lam_lin < 0);                     %% lower bound binding
0668 ku = find(lam_lin > 0);                     %% upper bound binding
0669 
0670 mu_l = zeros(nx+nA, 1);
0671 mu_l(ieq(kl)) = -lam_lin(kl);
0672 mu_l(igt) = mu_lin(nlt+(1:ngt));
0673 mu_l(ibx) = mu_lin(nlt+ngt+nbx+(1:nbx));
0674 
0675 mu_u = zeros(nx+nA, 1);
0676 mu_u(ieq(ku)) = lam_lin(ku);
0677 mu_u(ilt) = mu_lin(1:nlt);
0678 mu_u(ibx) = mu_lin(nlt+ngt+(1:nbx));
0679 
0680 fields = { ...
0681     'mu_l', mu_l((nx+1):end), ...
0682     'mu_u', mu_u((nx+1):end), ...
0683     'lower', mu_l(1:nx), ...
0684     'upper', mu_u(1:nx) ...
0685 };
0686 
0687 if niqnln > 0
0688     fields = { ...
0689         'ineqnonlin', mu(1:niqnln), ...
0690         fields{:} ...
0691     };
0692 end
0693 if neqnln > 0
0694     fields = { ...
0695         'eqnonlin', lam(1:neqnln), ...
0696         fields{:} ...
0697     };
0698 end
0699 
0700 lambda = struct(fields{:});
0701 
0702 % lambda = struct( ...
0703 %     'eqnonlin', lam(1:neqnln), ...
0704 %     'ineqnonlin', mu(1:niqnln), ...
0705 %     'mu_l', mu_l((nx+1):end), ...
0706 %     'mu_u', mu_u((nx+1):end), ...
0707 %     'lower', mu_l(1:nx), ...
0708 %     'upper', mu_u(1:nx) );

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