MIPS MATPOWER 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 https://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 MATPOWER 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 https://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-2018, Power Systems Engineering Research Center (PSERC) 0179 % by Ray Zimmerman, PSERC Cornell 0180 % 0181 % This file is part of MIPS. 0182 % Covered by the 3-clause BSD License (see LICENSE file for details). 0183 % See https://github.com/MATPOWER/mips for more info. 0184 0185 %%----- input argument handling ----- 0186 %% gather inputs 0187 if nargin == 1 && isstruct(f_fcn) %% problem struct 0188 p = f_fcn; 0189 f_fcn = p.f_fcn; 0190 x0 = p.x0; 0191 nx = size(x0, 1); %% number of optimization variables 0192 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0193 if isfield(p, 'hess_fcn'), hess_fcn = p.hess_fcn; else, hess_fcn = ''; end 0194 if isfield(p, 'gh_fcn'), gh_fcn = p.gh_fcn; else, gh_fcn = ''; end 0195 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0196 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0197 if isfield(p, 'u'), u = p.u; else, u = []; end 0198 if isfield(p, 'l'), l = p.l; else, l = []; end 0199 if isfield(p, 'A'), A = p.A; else, A=sparse(0,nx); end 0200 else %% individual args 0201 nx = size(x0, 1); %% number of optimization variables 0202 if nargin < 10 0203 opt = []; 0204 if nargin < 9 0205 hess_fcn = ''; 0206 if nargin < 8 0207 gh_fcn = ''; 0208 if nargin < 7 0209 xmax = []; 0210 if nargin < 6 0211 xmin = []; 0212 if nargin < 5 0213 u = []; 0214 if nargin < 4 0215 l = []; 0216 A = sparse(0,nx); 0217 end 0218 end 0219 end 0220 end 0221 end 0222 end 0223 end 0224 end 0225 %% set default argument values if missing 0226 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0227 (isempty(u) || all(u == Inf))) 0228 A = sparse(0,nx); %% no limits => no linear constraints 0229 end 0230 nA = size(A, 1); %% number of original linear constraints 0231 if isempty(u) %% By default, linear inequalities are ... 0232 u = Inf(nA, 1); %% ... unbounded above and ... 0233 end 0234 if isempty(l) 0235 l = -Inf(nA, 1); %% ... unbounded below. 0236 end 0237 if isempty(xmin) %% By default, optimization variables are ... 0238 xmin = -Inf(nx, 1); %% ... unbounded below and ... 0239 end 0240 if isempty(xmax) 0241 xmax = Inf(nx, 1); %% ... unbounded above. 0242 end 0243 if isempty(gh_fcn) 0244 nonlinear = false; %% no nonlinear constraints present 0245 gn = []; hn = []; 0246 else 0247 nonlinear = true; %% we have some nonlinear constraints 0248 end 0249 0250 %% default options 0251 if isempty(opt) 0252 opt = struct('verbose', 0); 0253 end 0254 if ~isfield(opt, 'linsolver') || isempty(opt.linsolver) 0255 opt.linsolver = ''; 0256 end 0257 if ~isfield(opt, 'feastol') || isempty(opt.feastol) || opt.feastol == 0 0258 opt.feastol = 1e-6; 0259 end 0260 if ~isfield(opt, 'gradtol') || isempty(opt.gradtol) || opt.gradtol == 0 0261 opt.gradtol = 1e-6; 0262 end 0263 if ~isfield(opt, 'comptol') || isempty(opt.comptol) || opt.comptol == 0 0264 opt.comptol = 1e-6; 0265 end 0266 if ~isfield(opt, 'costtol') || isempty(opt.costtol) || opt.costtol == 0 0267 opt.costtol = 1e-6; 0268 end 0269 if ~isfield(opt, 'max_it') || isempty(opt.max_it) 0270 opt.max_it = 150; 0271 end 0272 if ~isfield(opt, 'sc') || ~isfield(opt.sc, 'red_it') || isempty(opt.sc.red_it) 0273 opt.sc.red_it = 20; 0274 end 0275 if ~isfield(opt, 'step_control') || isempty(opt.step_control) 0276 opt.step_control = 0; 0277 end 0278 if ~isfield(opt, 'cost_mult') || isempty(opt.cost_mult) 0279 opt.cost_mult = 1; 0280 end 0281 if ~isfield(opt, 'verbose') || isempty(opt.verbose) 0282 opt.verbose = 0; 0283 end 0284 %% algorithm constants 0285 if ~isfield(opt, 'xi') || isempty(opt.xi) 0286 opt.xi = 0.99995; %% OPT_IPM_PHI 0287 end 0288 if ~isfield(opt, 'sigma') || isempty(opt.sigma) 0289 opt.sigma = 0.1; %% OPT_IPM_SIGMA, centering parameter 0290 end 0291 if ~isfield(opt, 'z0') || isempty(opt.z0) 0292 opt.z0 = 1; %% OPT_IPM_INIT_SLACK 0293 end 0294 if ~isfield(opt, 'alpha_min') || isempty(opt.alpha_min) 0295 opt.alpha_min = 1e-8; %% OPT_AP_AD_MIN 0296 end 0297 if ~isfield(opt, 'rho_min') || isempty(opt.rho_min) 0298 opt.rho_min = 0.95; %% OPT_IPM_QUAD_LOWTHRESH 0299 end 0300 if ~isfield(opt, 'rho_max') || isempty(opt.rho_max) 0301 opt.rho_max = 1.05; %% OPT_IPM_QUAD_HIGHTHRESH 0302 end 0303 if ~isfield(opt, 'mu_threshold') || isempty(opt.mu_threshold) 0304 opt.mu_threshold = 1e-5; %% SCOPF_MULTIPLIERS_FILTER_THRESH 0305 end 0306 if ~isfield(opt, 'max_stepsize') || isempty(opt.max_stepsize) 0307 opt.max_stepsize = 1e10; 0308 end 0309 0310 %% initialize history 0311 hist(opt.max_it+1) = struct('feascond', 0, 'gradcond', 0, 'compcond', 0, ... 0312 'costcond', 0, 'gamma', 0, 'stepsize', 0, 'obj', 0, ... 0313 'alphap', 0, 'alphad', 0); 0314 0315 %%----- set up problem ----- 0316 %% constants 0317 [xi, sigma, z0, alpha_min, rho_min, rho_max, mu_threshold, max_stepsize] = ... 0318 deal(opt.xi, opt.sigma, opt.z0, opt.alpha_min, ... 0319 opt.rho_min, opt.rho_max, opt.mu_threshold, opt.max_stepsize); 0320 if xi >= 1 || xi < 0.5 0321 error('mips: opt.xi (%g) must be a number slightly less than 1', opt.xi); 0322 end 0323 if sigma > 1 || sigma <= 0 0324 error('mips: opt.sigma (%g) must be a number between 0 and 1', opt.sigma); 0325 end 0326 0327 %% initialize 0328 i = 0; %% iteration counter 0329 converged = 0; %% flag 0330 eflag = 0; %% exit flag 0331 0332 %% add var limits to linear constraints 0333 AA = [speye(nx); A]; 0334 ll = [xmin; l]; 0335 uu = [xmax; u]; 0336 0337 %% split up linear constraints 0338 ieq = find( abs(uu-ll) <= eps ); %% equality 0339 igt = find( uu >= 1e10 & ll > -1e10 ); %% greater than, unbounded above 0340 ilt = find( ll <= -1e10 & uu < 1e10 ); %% less than, unbounded below 0341 ibx = find( (abs(uu-ll) > eps) & (uu < 1e10) & (ll > -1e10) ); 0342 Ae = AA(ieq, :); 0343 be = uu(ieq, 1); 0344 Ai = [ AA(ilt, :); -AA(igt, :); AA(ibx, :); -AA(ibx, :) ]; 0345 bi = [ uu(ilt, 1); -ll(igt, 1); uu(ibx, 1); -ll(ibx, 1) ]; 0346 0347 %% evaluate cost f(x0) and constraints g(x0), h(x0) 0348 x = x0; 0349 [f, df] = f_fcn(x); %% cost 0350 f = f * opt.cost_mult; 0351 df = df * opt.cost_mult; 0352 if nonlinear 0353 [hn, gn, dhn, dgn] = gh_fcn(x); %% nonlinear constraints 0354 h = [hn; Ai * x - bi]; %% inequality constraints 0355 g = [gn; Ae * x - be]; %% equality constraints 0356 dh = [dhn Ai']; %% 1st derivative of inequalities 0357 dg = [dgn Ae']; %% 1st derivative of equalities 0358 else 0359 h = Ai * x - bi; %% inequality constraints 0360 g = Ae * x - be; %% equality constraints 0361 dh = Ai'; %% 1st derivative of inequalities 0362 dg = Ae'; %% 1st derivative of equalities 0363 end 0364 0365 %% grab some dimensions 0366 neq = size(g, 1); %% number of equality constraints 0367 niq = size(h, 1); %% number of inequality constraints 0368 neqnln = size(gn, 1); %% number of nonlinear equality constraints 0369 niqnln = size(hn, 1); %% number of nonlinear inequality constraints 0370 nlt = length(ilt); %% number of upper bounded linear inequalities 0371 ngt = length(igt); %% number of lower bounded linear inequalities 0372 nbx = length(ibx); %% number of doubly bounded linear inequalities 0373 0374 %% initialize gamma, lam, mu, z, e 0375 gamma = 1; %% barrier coefficient, r in Harry's code 0376 lam = zeros(neq, 1); 0377 z = z0 * ones(niq, 1); 0378 mu = z; 0379 k = find(h < -z0); 0380 z(k) = -h(k); 0381 k = find(gamma / z > z0); %% (seems k is always empty if gamma = z0 = 1) 0382 if ~isempty(k) 0383 mu(k) = gamma / z(k); 0384 end 0385 e = ones(niq, 1); 0386 0387 %% check tolerance 0388 f0 = f; 0389 if opt.step_control 0390 L = f + lam' * g + mu' * (h+z) - gamma * sum(log(z)); 0391 end 0392 Lx = df + dg * lam + dh * mu; 0393 if isempty(h) 0394 maxh = zeros(1,0); 0395 else 0396 maxh = max(h); 0397 end 0398 feascond = max([norm(g, Inf), maxh]) / (1 + max([ norm(x, Inf), norm(z, Inf) ])); 0399 gradcond = norm(Lx, Inf) / (1 + max([ norm(lam, Inf), norm(mu, Inf) ])); 0400 compcond = (z' * mu) / (1 + norm(x, Inf)); 0401 costcond = abs(f - f0) / (1 + abs(f0)); 0402 %% save history 0403 hist(i+1) = struct('feascond', feascond, 'gradcond', gradcond, ... 0404 'compcond', compcond, 'costcond', costcond, 'gamma', gamma, ... 0405 'stepsize', 0, 'obj', f/opt.cost_mult, 'alphap', 0, 'alphad', 0); 0406 if strcmp(upper(opt.linsolver), 'PARDISO') 0407 ls = 'PARDISO'; 0408 mplinsolve_opt = struct('pardiso', ... 0409 struct('mtype', -2, ... 0410 'iparm', [8, 1; %% max it refinement steps 0411 10, 12;%% eps pivot 0412 11, 1; %% use scaling vectors 0413 13, 2; %% improved accuracy 0414 18, 0; %% do not determine nnz in LU 0415 21, 3; %% ? undocumented pivoting 0416 ] )); 0417 if ~have_fcn('pardiso') 0418 warning('mips: PARDISO linear solver not available, using default'); 0419 opt.linsolver = ''; 0420 ls = 'built-in'; 0421 end 0422 else 0423 ls = 'built-in'; 0424 mplinsolve_opt = []; 0425 end 0426 if opt.verbose 0427 if opt.step_control, s = '-sc'; else, s = ''; end 0428 v = mipsver('all'); 0429 fprintf('MATPOWER Interior Point Solver -- MIPS%s, Version %s, %s\n (using %s linear solver)', ... 0430 s, v.Version, v.Date, ls); 0431 if opt.verbose > 1 0432 fprintf('\n it objective step size feascond gradcond compcond costcond '); 0433 fprintf('\n---- ------------ --------- ------------ ------------ ------------ ------------'); 0434 fprintf('\n%3d %12.8g %10s %12g %12g %12g %12g', ... 0435 i, f/opt.cost_mult, '', feascond, gradcond, compcond, costcond); 0436 end 0437 end 0438 if feascond < opt.feastol && gradcond < opt.gradtol && ... 0439 compcond < opt.comptol && costcond < opt.costtol 0440 converged = 1; 0441 if opt.verbose 0442 fprintf('\nConverged!\n'); 0443 end 0444 end 0445 0446 %%----- do Newton iterations ----- 0447 while (~converged && i < opt.max_it) 0448 %% update iteration counter 0449 i = i + 1; 0450 0451 %% compute update step 0452 lambda = struct('eqnonlin', lam(1:neqnln), 'ineqnonlin', mu(1:niqnln)); 0453 if nonlinear 0454 if isempty(hess_fcn) 0455 fprintf('mips: Hessian evaluation via finite differences not yet implemented.\n Please provide your own hessian evaluation function.'); 0456 end 0457 Lxx = hess_fcn(x, lambda, opt.cost_mult); 0458 else 0459 [f_, df_, d2f] = f_fcn(x); %% cost 0460 Lxx = d2f * opt.cost_mult; 0461 end 0462 zinvdiag = sparse(1:niq, 1:niq, 1 ./ z, niq, niq); 0463 mudiag = sparse(1:niq, 1:niq, mu, niq, niq); 0464 dh_zinv = dh * zinvdiag; 0465 M = Lxx + dh_zinv * mudiag * dh'; 0466 N = Lx + dh_zinv * (mudiag * h + gamma * e); 0467 dxdlam = mplinsolve([M dg; dg' sparse(neq, neq)], [-N; -g], opt.linsolver, mplinsolve_opt); 0468 % AAA = [ 0469 % M dg; 0470 % dg' sparse(neq, neq) 0471 % ]; 0472 % rc = 1/condest(AAA); 0473 % if rc < 1e-22 0474 % fprintf('my RCOND = %g\n', rc); 0475 % n = size(AAA, 1); 0476 % AAA = AAA + 1e-3 * speye(n,n); 0477 % end 0478 % bbb = [-N; -g]; 0479 % dxdlam = AAA \ bbb; 0480 if any(isnan(dxdlam)) || norm(dxdlam) > max_stepsize 0481 if opt.verbose 0482 fprintf('\nNumerically Failed\n'); 0483 end 0484 eflag = -1; 0485 break; 0486 end 0487 dx = dxdlam(1:nx, 1); 0488 dlam = dxdlam(nx+(1:neq), 1); 0489 dz = -h - z - dh' * dx; 0490 dmu = -mu + zinvdiag *(gamma*e - mudiag * dz); 0491 0492 %% optional step-size control 0493 sc = 0; 0494 if opt.step_control 0495 x1 = x + dx; 0496 0497 %% evaluate cost, constraints, derivatives at x1 0498 [f1, df1] = f_fcn(x1); %% cost 0499 f1 = f1 * opt.cost_mult; 0500 df1 = df1 * opt.cost_mult; 0501 if nonlinear 0502 [hn1, gn1, dhn1, dgn1] = gh_fcn(x1); %% nonlinear constraints 0503 h1 = [hn1; Ai * x1 - bi]; %% inequality constraints 0504 g1 = [gn1; Ae * x1 - be]; %% equality constraints 0505 dh1 = [dhn1 Ai']; %% 1st derivative of inequalities 0506 dg1 = [dgn1 Ae']; %% 1st derivative of equalities 0507 else 0508 h1 = Ai * x1 - bi; %% inequality constraints 0509 g1 = Ae * x1 - be; %% equality constraints 0510 dh1 = dh; %% 1st derivative of inequalities 0511 dg1 = dg; %% 1st derivative of equalities 0512 end 0513 0514 %% check tolerance 0515 Lx1 = df1 + dg1 * lam + dh1 * mu; 0516 if isempty(h1) 0517 maxh1 = zeros(1,0); 0518 else 0519 maxh1 = max(h1); 0520 end 0521 feascond1 = max([norm(g1, Inf), maxh1]) / (1 + max([ norm(x1, Inf), norm(z, Inf) ])); 0522 gradcond1 = norm(Lx1, Inf) / (1 + max([ norm(lam, Inf), norm(mu, Inf) ])); 0523 0524 if feascond1 > feascond && gradcond1 > gradcond 0525 sc = 1; 0526 end 0527 end 0528 if sc 0529 alpha = 1; 0530 for j = 1:opt.sc.red_it 0531 dx1 = alpha * dx; 0532 x1 = x + dx1; 0533 f1 = f_fcn(x1); %% cost 0534 f1 = f1 * opt.cost_mult; 0535 if nonlinear 0536 [hn1, gn1] = gh_fcn(x1); %% nonlinear constraints 0537 h1 = [hn1; Ai * x1 - bi]; %% inequality constraints 0538 g1 = [gn1; Ae * x1 - be]; %% equality constraints 0539 else 0540 h1 = Ai * x1 - bi; %% inequality constraints 0541 g1 = Ae * x1 - be; %% equality constraints 0542 end 0543 L1 = f1 + lam' * g1 + mu' * (h1+z) - gamma * sum(log(z)); 0544 if opt.verbose > 2 0545 fprintf('\n %3d %10g', -j, norm(dx1)); 0546 end 0547 rho = (L1 - L) / (Lx' * dx1 + 0.5 * dx1' * Lxx * dx1); 0548 if rho > rho_min && rho < rho_max 0549 break; 0550 else 0551 alpha = alpha / 2; 0552 end 0553 end 0554 dx = alpha * dx; 0555 dz = alpha * dz; 0556 dlam = alpha * dlam; 0557 dmu = alpha * dmu; 0558 end 0559 0560 %% do the update 0561 k = find(dz < 0); 0562 if isempty(k) 0563 alphap = 1; 0564 else 0565 alphap = min( [xi * min(z(k) ./ -dz(k)) 1] ); 0566 end 0567 k = find(dmu < 0); 0568 if isempty(k) 0569 alphad = 1; 0570 else 0571 alphad = min( [xi * min(mu(k) ./ -dmu(k)) 1] ); 0572 end 0573 x = x + alphap * dx; 0574 z = z + alphap * dz; 0575 lam = lam + alphad * dlam; 0576 mu = mu + alphad * dmu; 0577 if niq > 0 0578 gamma = sigma * (z' * mu) / niq; 0579 end 0580 0581 %% evaluate cost, constraints, derivatives 0582 [f, df] = f_fcn(x); %% cost 0583 f = f * opt.cost_mult; 0584 df = df * opt.cost_mult; 0585 if nonlinear 0586 [hn, gn, dhn, dgn] = gh_fcn(x); %% nonlinear constraints 0587 h = [hn; Ai * x - bi]; %% inequality constraints 0588 g = [gn; Ae * x - be]; %% equality constraints 0589 dh = [dhn Ai']; %% 1st derivative of inequalities 0590 dg = [dgn Ae']; %% 1st derivative of equalities 0591 else 0592 h = Ai * x - bi; %% inequality constraints 0593 g = Ae * x - be; %% equality constraints 0594 %% 1st derivatives are constant, still dh = Ai', dg = Ae' 0595 end 0596 0597 %% check tolerance 0598 Lx = df + dg * lam + dh * mu; 0599 if isempty(h) 0600 maxh = zeros(1,0); 0601 else 0602 maxh = max(h); 0603 end 0604 feascond = max([norm(g, Inf), maxh]) / (1 + max([ norm(x, Inf), norm(z, Inf) ])); 0605 gradcond = norm(Lx, Inf) / (1 + max([ norm(lam, Inf), norm(mu, Inf) ])); 0606 compcond = (z' * mu) / (1 + norm(x, Inf)); 0607 costcond = abs(f - f0) / (1 + abs(f0)); 0608 %% save history 0609 hist(i+1) = struct('feascond', feascond, 'gradcond', gradcond, ... 0610 'compcond', compcond, 'costcond', costcond, 'gamma', gamma, ... 0611 'stepsize', norm(dx), 'obj', f/opt.cost_mult, ... 0612 'alphap', alphap, 'alphad', alphad); 0613 0614 if opt.verbose > 1 0615 fprintf('\n%3d %12.8g %10.5g %12g %12g %12g %12g', ... 0616 i, f/opt.cost_mult, norm(dx), feascond, gradcond, compcond, costcond); 0617 end 0618 if feascond < opt.feastol && gradcond < opt.gradtol && ... 0619 compcond < opt.comptol && costcond < opt.costtol 0620 converged = 1; 0621 if opt.verbose 0622 fprintf('\nConverged!\n'); 0623 end 0624 else 0625 if any(isnan(x)) || alphap < alpha_min || alphad < alpha_min || ... 0626 gamma < eps || gamma > 1/eps 0627 if opt.verbose 0628 fprintf('\nNumerically Failed\n'); 0629 end 0630 eflag = -1; 0631 break; 0632 end 0633 f0 = f; 0634 if opt.step_control 0635 L = f + lam' * g + mu' * (h+z) - gamma * sum(log(z)); 0636 end 0637 end 0638 end 0639 0640 if opt.verbose 0641 if ~converged 0642 fprintf('\nDid not converge in %d iterations.\n', i); 0643 end 0644 end 0645 0646 %%----- package up results ----- 0647 hist = hist(1:i+1); 0648 if eflag ~= -1 0649 eflag = converged; 0650 end 0651 output = struct('iterations', i, 'hist', hist, 'message', ''); 0652 if eflag == 0 0653 output.message = 'Did not converge'; 0654 elseif eflag == 1 0655 output.message = 'Converged'; 0656 elseif eflag == -1 0657 output.message = 'Numerically failed'; 0658 else 0659 output.message = 'Please hang up and dial again'; 0660 end 0661 0662 %% zero out multipliers on non-binding constraints 0663 mu(h < -opt.feastol & mu < mu_threshold) = 0; 0664 0665 %% un-scale cost and prices 0666 f = f / opt.cost_mult; 0667 lam = lam / opt.cost_mult; 0668 mu = mu / opt.cost_mult; 0669 0670 %% re-package multipliers into struct 0671 lam_lin = lam((neqnln+1):neq); %% lambda for linear constraints 0672 mu_lin = mu((niqnln+1):niq); %% mu for linear constraints 0673 kl = find(lam_lin < 0); %% lower bound binding 0674 ku = find(lam_lin > 0); %% upper bound binding 0675 0676 mu_l = zeros(nx+nA, 1); 0677 mu_l(ieq(kl)) = -lam_lin(kl); 0678 mu_l(igt) = mu_lin(nlt+(1:ngt)); 0679 mu_l(ibx) = mu_lin(nlt+ngt+nbx+(1:nbx)); 0680 0681 mu_u = zeros(nx+nA, 1); 0682 mu_u(ieq(ku)) = lam_lin(ku); 0683 mu_u(ilt) = mu_lin(1:nlt); 0684 mu_u(ibx) = mu_lin(nlt+ngt+(1:nbx)); 0685 0686 fields = { ... 0687 'mu_l', mu_l((nx+1):end), ... 0688 'mu_u', mu_u((nx+1):end), ... 0689 'lower', mu_l(1:nx), ... 0690 'upper', mu_u(1:nx) ... 0691 }; 0692 0693 if niqnln > 0 0694 fields = { ... 0695 'ineqnonlin', mu(1:niqnln), ... 0696 fields{:} ... 0697 }; 0698 end 0699 if neqnln > 0 0700 fields = { ... 0701 'eqnonlin', lam(1:neqnln), ... 0702 fields{:} ... 0703 }; 0704 end 0705 0706 lambda = struct(fields{:}); 0707 0708 % lambda = struct( ... 0709 % 'eqnonlin', lam(1:neqnln), ... 0710 % 'ineqnonlin', mu(1:niqnln), ... 0711 % 'mu_l', mu_l((nx+1):end), ... 0712 % 'mu_u', mu_u((nx+1):end), ... 0713 % 'lower', mu_l(1:nx), ... 0714 % 'upper', mu_u(1:nx) );