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 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 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 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-2016, 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 http://www.pserc.cornell.edu/matpower/ 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 if ~have_fcn('pardiso') 0409 warning('mips: PARDISO linear solver not available, using default'); 0410 opt.linsolver = ''; 0411 ls = 'built-in'; 0412 end 0413 else 0414 ls = 'built-in'; 0415 end 0416 if opt.verbose 0417 if opt.step_control, s = '-sc'; else, s = ''; end 0418 v = mipsver('all'); 0419 fprintf('MATPOWER Interior Point Solver -- MIPS%s, Version %s, %s\n (using %s linear solver)', ... 0420 s, v.Version, v.Date, ls); 0421 if opt.verbose > 1 0422 fprintf('\n it objective step size feascond gradcond compcond costcond '); 0423 fprintf('\n---- ------------ --------- ------------ ------------ ------------ ------------'); 0424 fprintf('\n%3d %12.8g %10s %12g %12g %12g %12g', ... 0425 i, f/opt.cost_mult, '', feascond, gradcond, compcond, costcond); 0426 end 0427 end 0428 if feascond < opt.feastol && gradcond < opt.gradtol && ... 0429 compcond < opt.comptol && costcond < opt.costtol 0430 converged = 1; 0431 if opt.verbose 0432 fprintf('\nConverged!\n'); 0433 end 0434 end 0435 0436 %%----- do Newton iterations ----- 0437 while (~converged && i < opt.max_it) 0438 %% update iteration counter 0439 i = i + 1; 0440 0441 %% compute update step 0442 lambda = struct('eqnonlin', lam(1:neqnln), 'ineqnonlin', mu(1:niqnln)); 0443 if nonlinear 0444 if isempty(hess_fcn) 0445 fprintf('mips: Hessian evaluation via finite differences not yet implemented.\n Please provide your own hessian evaluation function.'); 0446 end 0447 Lxx = hess_fcn(x, lambda, opt.cost_mult); 0448 else 0449 [f_, df_, d2f] = f_fcn(x); %% cost 0450 Lxx = d2f * opt.cost_mult; 0451 end 0452 zinvdiag = sparse(1:niq, 1:niq, 1 ./ z, niq, niq); 0453 mudiag = sparse(1:niq, 1:niq, mu, niq, niq); 0454 dh_zinv = dh * zinvdiag; 0455 M = Lxx + dh_zinv * mudiag * dh'; 0456 N = Lx + dh_zinv * (mudiag * h + gamma * e); 0457 dxdlam = mplinsolve([M dg; dg' sparse(neq, neq)], [-N; -g], opt.linsolver, []); 0458 % AAA = [ 0459 % M dg; 0460 % dg' sparse(neq, neq) 0461 % ]; 0462 % rc = 1/condest(AAA); 0463 % if rc < 1e-22 0464 % fprintf('my RCOND = %g\n', rc); 0465 % n = size(AAA, 1); 0466 % AAA = AAA + 1e-3 * speye(n,n); 0467 % end 0468 % bbb = [-N; -g]; 0469 % dxdlam = AAA \ bbb; 0470 if any(isnan(dxdlam)) || norm(dxdlam) > max_stepsize 0471 if opt.verbose 0472 fprintf('\nNumerically Failed\n'); 0473 end 0474 eflag = -1; 0475 break; 0476 end 0477 dx = dxdlam(1:nx, 1); 0478 dlam = dxdlam(nx+(1:neq), 1); 0479 dz = -h - z - dh' * dx; 0480 dmu = -mu + zinvdiag *(gamma*e - mudiag * dz); 0481 0482 %% optional step-size control 0483 sc = 0; 0484 if opt.step_control 0485 x1 = x + dx; 0486 0487 %% evaluate cost, constraints, derivatives at x1 0488 [f1, df1] = f_fcn(x1); %% cost 0489 f1 = f1 * opt.cost_mult; 0490 df1 = df1 * opt.cost_mult; 0491 if nonlinear 0492 [hn1, gn1, dhn1, dgn1] = gh_fcn(x1); %% nonlinear constraints 0493 h1 = [hn1; Ai * x1 - bi]; %% inequality constraints 0494 g1 = [gn1; Ae * x1 - be]; %% equality constraints 0495 dh1 = [dhn1 Ai']; %% 1st derivative of inequalities 0496 dg1 = [dgn1 Ae']; %% 1st derivative of equalities 0497 else 0498 h1 = Ai * x1 - bi; %% inequality constraints 0499 g1 = Ae * x1 - be; %% equality constraints 0500 dh1 = dh; %% 1st derivative of inequalities 0501 dg1 = dg; %% 1st derivative of equalities 0502 end 0503 0504 %% check tolerance 0505 Lx1 = df1 + dg1 * lam + dh1 * mu; 0506 if isempty(h1) 0507 maxh1 = zeros(1,0); 0508 else 0509 maxh1 = max(h1); 0510 end 0511 feascond1 = max([norm(g1, Inf), maxh1]) / (1 + max([ norm(x1, Inf), norm(z, Inf) ])); 0512 gradcond1 = norm(Lx1, Inf) / (1 + max([ norm(lam, Inf), norm(mu, Inf) ])); 0513 0514 if feascond1 > feascond && gradcond1 > gradcond 0515 sc = 1; 0516 end 0517 end 0518 if sc 0519 alpha = 1; 0520 for j = 1:opt.sc.red_it 0521 dx1 = alpha * dx; 0522 x1 = x + dx1; 0523 f1 = f_fcn(x1); %% cost 0524 f1 = f1 * opt.cost_mult; 0525 if nonlinear 0526 [hn1, gn1] = gh_fcn(x1); %% nonlinear constraints 0527 h1 = [hn1; Ai * x1 - bi]; %% inequality constraints 0528 g1 = [gn1; Ae * x1 - be]; %% equality constraints 0529 else 0530 h1 = Ai * x1 - bi; %% inequality constraints 0531 g1 = Ae * x1 - be; %% equality constraints 0532 end 0533 L1 = f1 + lam' * g1 + mu' * (h1+z) - gamma * sum(log(z)); 0534 if opt.verbose > 2 0535 fprintf('\n %3d %10g', -j, norm(dx1)); 0536 end 0537 rho = (L1 - L) / (Lx' * dx1 + 0.5 * dx1' * Lxx * dx1); 0538 if rho > rho_min && rho < rho_max 0539 break; 0540 else 0541 alpha = alpha / 2; 0542 end 0543 end 0544 dx = alpha * dx; 0545 dz = alpha * dz; 0546 dlam = alpha * dlam; 0547 dmu = alpha * dmu; 0548 end 0549 0550 %% do the update 0551 k = find(dz < 0); 0552 if isempty(k) 0553 alphap = 1; 0554 else 0555 alphap = min( [xi * min(z(k) ./ -dz(k)) 1] ); 0556 end 0557 k = find(dmu < 0); 0558 if isempty(k) 0559 alphad = 1; 0560 else 0561 alphad = min( [xi * min(mu(k) ./ -dmu(k)) 1] ); 0562 end 0563 x = x + alphap * dx; 0564 z = z + alphap * dz; 0565 lam = lam + alphad * dlam; 0566 mu = mu + alphad * dmu; 0567 if niq > 0 0568 gamma = sigma * (z' * mu) / niq; 0569 end 0570 0571 %% evaluate cost, constraints, derivatives 0572 [f, df] = f_fcn(x); %% cost 0573 f = f * opt.cost_mult; 0574 df = df * opt.cost_mult; 0575 if nonlinear 0576 [hn, gn, dhn, dgn] = gh_fcn(x); %% nonlinear constraints 0577 h = [hn; Ai * x - bi]; %% inequality constraints 0578 g = [gn; Ae * x - be]; %% equality constraints 0579 dh = [dhn Ai']; %% 1st derivative of inequalities 0580 dg = [dgn Ae']; %% 1st derivative of equalities 0581 else 0582 h = Ai * x - bi; %% inequality constraints 0583 g = Ae * x - be; %% equality constraints 0584 %% 1st derivatives are constant, still dh = Ai', dg = Ae' 0585 end 0586 0587 %% check tolerance 0588 Lx = df + dg * lam + dh * mu; 0589 if isempty(h) 0590 maxh = zeros(1,0); 0591 else 0592 maxh = max(h); 0593 end 0594 feascond = max([norm(g, Inf), maxh]) / (1 + max([ norm(x, Inf), norm(z, Inf) ])); 0595 gradcond = norm(Lx, Inf) / (1 + max([ norm(lam, Inf), norm(mu, Inf) ])); 0596 compcond = (z' * mu) / (1 + norm(x, Inf)); 0597 costcond = abs(f - f0) / (1 + abs(f0)); 0598 %% save history 0599 hist(i+1) = struct('feascond', feascond, 'gradcond', gradcond, ... 0600 'compcond', compcond, 'costcond', costcond, 'gamma', gamma, ... 0601 'stepsize', norm(dx), 'obj', f/opt.cost_mult, ... 0602 'alphap', alphap, 'alphad', alphad); 0603 0604 if opt.verbose > 1 0605 fprintf('\n%3d %12.8g %10.5g %12g %12g %12g %12g', ... 0606 i, f/opt.cost_mult, norm(dx), feascond, gradcond, compcond, costcond); 0607 end 0608 if feascond < opt.feastol && gradcond < opt.gradtol && ... 0609 compcond < opt.comptol && costcond < opt.costtol 0610 converged = 1; 0611 if opt.verbose 0612 fprintf('\nConverged!\n'); 0613 end 0614 else 0615 if any(isnan(x)) || alphap < alpha_min || alphad < alpha_min || ... 0616 gamma < eps || gamma > 1/eps 0617 if opt.verbose 0618 fprintf('\nNumerically Failed\n'); 0619 end 0620 eflag = -1; 0621 break; 0622 end 0623 f0 = f; 0624 if opt.step_control 0625 L = f + lam' * g + mu' * (h+z) - gamma * sum(log(z)); 0626 end 0627 end 0628 end 0629 0630 if opt.verbose 0631 if ~converged 0632 fprintf('\nDid not converge in %d iterations.\n', i); 0633 end 0634 end 0635 0636 %%----- package up results ----- 0637 hist = hist(1:i+1); 0638 if eflag ~= -1 0639 eflag = converged; 0640 end 0641 output = struct('iterations', i, 'hist', hist, 'message', ''); 0642 if eflag == 0 0643 output.message = 'Did not converge'; 0644 elseif eflag == 1 0645 output.message = 'Converged'; 0646 elseif eflag == -1 0647 output.message = 'Numerically failed'; 0648 else 0649 output.message = 'Please hang up and dial again'; 0650 end 0651 0652 %% zero out multipliers on non-binding constraints 0653 mu(h < -opt.feastol & mu < mu_threshold) = 0; 0654 0655 %% un-scale cost and prices 0656 f = f / opt.cost_mult; 0657 lam = lam / opt.cost_mult; 0658 mu = mu / opt.cost_mult; 0659 0660 %% re-package multipliers into struct 0661 lam_lin = lam((neqnln+1):neq); %% lambda for linear constraints 0662 mu_lin = mu((niqnln+1):niq); %% mu for linear constraints 0663 kl = find(lam_lin < 0); %% lower bound binding 0664 ku = find(lam_lin > 0); %% upper bound binding 0665 0666 mu_l = zeros(nx+nA, 1); 0667 mu_l(ieq(kl)) = -lam_lin(kl); 0668 mu_l(igt) = mu_lin(nlt+(1:ngt)); 0669 mu_l(ibx) = mu_lin(nlt+ngt+nbx+(1:nbx)); 0670 0671 mu_u = zeros(nx+nA, 1); 0672 mu_u(ieq(ku)) = lam_lin(ku); 0673 mu_u(ilt) = mu_lin(1:nlt); 0674 mu_u(ibx) = mu_lin(nlt+ngt+(1:nbx)); 0675 0676 fields = { ... 0677 'mu_l', mu_l((nx+1):end), ... 0678 'mu_u', mu_u((nx+1):end), ... 0679 'lower', mu_l(1:nx), ... 0680 'upper', mu_u(1:nx) ... 0681 }; 0682 0683 if niqnln > 0 0684 fields = { ... 0685 'ineqnonlin', mu(1:niqnln), ... 0686 fields{:} ... 0687 }; 0688 end 0689 if neqnln > 0 0690 fields = { ... 0691 'eqnonlin', lam(1:neqnln), ... 0692 fields{:} ... 0693 }; 0694 end 0695 0696 lambda = struct(fields{:}); 0697 0698 % lambda = struct( ... 0699 % 'eqnonlin', lam(1:neqnln), ... 0700 % 'ineqnonlin', mu(1:niqnln), ... 0701 % 'mu_l', mu_l((nx+1):end), ... 0702 % 'mu_u', mu_u((nx+1):end), ... 0703 % 'lower', mu_l(1:nx), ... 0704 % 'upper', mu_u(1:nx) );