Home > matpower4.0 > mips6.m

mips6

PURPOSE ^

------------------------------ deprecated ------------------------------

SYNOPSIS ^

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

DESCRIPTION ^

------------------------------  deprecated  ------------------------------
   MATLAB 6.x support to be removed in a future version.
--------------------------------------------------------------------------
MIPS6  MATLAB Interior Point Solver (for MATLAB 6.x).
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
       MIPS6(F_FCN, X0, A, L, U, XMIN, XMAX, GH_FCN, HESS_FCN, OPT, VARARGIN)
   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] = ...
           mips6(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, ...
                   opt, varargin);

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

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

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