------------------------------ 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.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) );