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