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