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

Generated on Mon 26-Jan-2015 14:56:45 by m2html © 2005