Page 1

	

Optimal	Power	Flow	Paper	4	

A	staff	paper	by

Anya	Castillo

Richard	P.	O’Neill	

March 2013

The views presented are the personal views of the authors and not the

Federal Energy Regulatory Commission or any of its Commissioners.
Page 2

Survey	of	Approaches	to	Solving	the	ACOPF	

Optimal	Power	Flow	Paper	4	

Anya	Castillo	and	Richard	P	O’Neill	

anya.castillo@ferc.gov;	richard.oneill@ferc.gov	

March	2013			

	

Abstract:	

We	present	a	background	on	approaches	historically	applied	to	solve	the	ACOPF,	

many	which	are	used	in	our	following	companion	study	on	testing	solution	

techniques	ሺCastillo,	2013ሻ.	In	this	paper	we	present	an	introduction	on	the	

associated	theory	in	nonlinear	optimization,	and	then	discuss	the	solvers	and	

published	algorithms	that	have	been	applied	to	the	ACOPF,	dating	initially	from	

Carpentier	in	1962	to	current	day	approaches.	We	provide	insight	into	the	major	

contributions	in	solution	methods	applied	to	the	ACOPF	to	date.			

	

Disclaimer:	The	views	presented	are	the	personal	views	of	the	authors	and	not	the	

Federal	Energy	Regulatory	Commission	or	any	of	its	Commissioners.	 	
Page 3

Table	of	Contents	

	

1. Introduction		 	 	 	 	 	 	 	 	 4	

2. Definitions	 	 	 	 	 	 	 	 	 4	

3. Unconstrained	Nonlinear	Optimization		 	 	 	 	 5	

4. Constrained	Nonlinear	Optimization	 	 	 	 	 	 14	

5. Decomposition	Techniques	 	 	 	 	 	 	 23	

6. Commercial	Solvers	 	 	 	 	 	 	 	 28	

7. History	of	ACOPF	Solution	Techniques	 	 	 	 	 	 29	

8. Summary	and	Conclusions	 	 	 	 	 	 	 40	

References	 	 	 	 	 	 	 	 	 	 41	

	 	
Page 4

1.		Introduction	

	

	 Recently,	mathematical	optimization	has	become	a	powerful	tool	in	many	

real‐life	applications	along	with	the	more	traditional	uses	in	engineering	and	

science.	The	history	is	driven,	and	will	continue	to	be	driven,	by	the	demand	to	solve	

increasingly	larger	nonlinear	optimization	problems	under	the	limitations	of	

storage	and	computing	time.	

	 Here,	we	present	background	on	approaches	that	are	used	in	a	companion	

study	and	are	often	applied	to	solve	the	ACOPF.	We	are	focused	on	understanding	

the	reliability	of	such	approaches	and	whether	there	is	good	reason	to	believe	that	

convergence	of	the	ACOPF	would	be	guaranteed	at	a	reasonable	rate.	Therefore,	we	

introduce	and	discuss	some	of	the	associated	theory.	

	 In	section	2,	we	present	definitions.	In	section	3,	we	review	unconstrained	

optimization	for	nonlinear	programming	followed	by	a	review	of	constrained	

optimization	for	nonlinear	programming	in	section	4;	there	are	numerous	nonlinear	

programming	resources	and,	for	example,	the	reader	can	refer	to	ሺNocedal,	1999ሻ,	

ሺPanos,	2006ሻ	for	a	more	detailed	discussion.	In	section	5,	we	review	problem	

decomposition	techniques	and	popular	approximations	to	the	ACOPF.	In	section	6,	

we	review	the	nonlinear	commercial	solvers	applied	to	the	ACOPF	in	our	companion	

computational	study	ሺCastillo,	2013ሻ.	Then	in	section	7,	we	review	the	literature	on	

solving	the	ACOPF	followed	by	a	brief	discussion	in	section	8.	

2.		Definitions	

	

Definition	1.	If	x	є	X		Rn,	x	is	called	a	feasible	solution.		

Definition	2.	The	function	fሺxሻ	with	domain	Rn	and	codomain	R	is	denoted	by		

fሺxሻ:	Rn →R.	

Definition	3.	If	fሺx*ሻ	൑	fሺxሻ	for	all	x	є	X,	x*	is	called	an	optimal	solution.		

Definition	4.	If	fሺx*ሻ	൑	fሺxሻ	for	all	x	є	X’  X,	x*	is	called	a	local	optimal	solution.		

Definition	5. The	set	TR  X	is	called	a	trust	region	where	an	approximation	of	fሺxሻ	is	

‘trusted.’	

Definition	6. The	set	X	is	convex,	if	and	only	if	αx1	൅	ሺ1	‐	αሻx2	є	X	for	all	α	є	ሾ0,1ሿ	and	

all	x1,	x2	є	X.	

Definition	7. The	function	fሺxሻ	is	convex	in	domain	X	if	and	only	if	fሺαx1	൅	ሺ1	‐	αሻx2ሻ	

൑	αfሺx1ሻ	൅	ሺ1	‐	αሻfሺx2ሻ	for	any	α	є	ሾ0,1ሿ	and	all	x1,	x2	є	X.	

Definition	8. If	d	є	Rn	and	x	൅	sd	є	X	for	some	s	൐	0,	d	is	called	a	feasible	direction.	
Page 5

Definition	9.	The	function	Fሺxkሻ	ൌ	xk൅1	is	the	result	of	iteration	k	and	the	approach	is	

a	descent	algorithm	if	fሺxk൅1ሻ	൏	fሺxkሻ	ሺthe	reduction	conditionሻ.	

Defintion	10.	The	L2‐norm	is	|	|.	

Definition	11.	If	fሺxሻ	is	once	differentiable,	the	vector	of	first	derivatives	of	fሺxሻ	with	

respect	to	x	is	fሺxሻ;	fሺxሻ	is	also	called	the	Jacobian	or	gradient	vector.	

Definition	12.	If	fሺxሻ	is	twice	differentiable,	the	matrix	of	second	derivatives	of	fሺxሻ	

with	respect	to	x	is	Hሺxሻ	ൌ	2fሺxሻ;	Hሺxሻ	is	called	the	Hessian.	

Definition	13.	If	Q	is	a	symmetric	matrix	and	d1

TQd2	ൌ	0,	d1	and	d2	are	Q‐orthogonal,	

also	known	as	conjugate	directions.		

Definition	14.	If	ሺA	‐	eIሻv	ൌ	0,	e	is	an	eigenvalue	and	v	is	an	eigenvector	of	A.	

Definition	15.	For	unconstrained	optimization	we	solve	inf	ሼfሺxሻ|	x	є	Xሽ	where	the	

objective	function	fሺxሻ:	Rn →	R

Definition	16.	For	equality	constrained	optimization	we	solve	inf	ሼfሺxሻ|	gሺxሻ	ൌ	0,	x	є	

Xሽ	where	the	constraints	gሺxሻ:	Rn →	Rmሺm	൑	nሻ.	

Definition	17.	For	inequality	constrained	optimization	we	solve	inf	ሼfሺxሻ|	gሺxሻ	൑	0,	x

є	Xሽ	where	the	constraints	gሺxሻ:	Rn →	Rm ሺm	may	be	larger	than	nሻ.	

Definition	18.	The	index	set	AS	is	the	active	set	of	constraints	where	ASሺxሻ	ൌ	ሼj	|	

gjሺxሻ	ൌ	0ሽ.	

Definition	19.	The	Lagrangian	function	is	Lሺx,	λሻ	ൌ	fሺxሻ	൅	λTgሺxሻ.				

Definition	20.	The	Lagrange	multiplier	λ	exists	such	that	ሺx*,	λ*ሻ	is	a	stationary	point	

of	the	Lagrangian	function	where	λLሺx*,	λ*ሻ	ൌ	0,	which	implies	gሺx*ሻ	ൌ	0.			

3.	Unconstrained	Nonlinear	Optimization		

	

We	start	with	the	‘unconstrained’	nonlinear	optimization	problem.	We	want	to	

minimize	the	nonlinear	function	fሺxሻ:		

f*ൌ	inf	ሼfሺxሻ|	x	є	Xሽ	 	 	

where	x*	is	a	solution	to	the	above	problem	ሺPሻ.	Methods	for	solving	this	problem	

date	back	to	Newton	and	Gauss.	In	unconstrained	optimization,	the	set	X	is	Rn.	If	X	is	

a	subset	of	Rn,	X	could	be	the	set	of	integers	Zn	or	denote	a	constraint	set	of	

equations.	If	P	is	infeasible,	we	define	f*ൌ	∞;	if	P	is	unbounded,	we	define	f*ൌ	‐∞.	

The	function	f	is	typically	called	the	objective	function.	Equivalently,	the	problem	

can	be	stated	as	a	maximization	of	‐fሺxሻ.	We	present	the	following	theorems	without	

proof.	For	more	details	see	Mangasarian	ሺ1969ሻ	and	Luenberger	ሺ1984ሻ.		
Page 6

Theorem	1.	If	fሺxሻ	is	differentiable,	x*	is	a	relative	or	local	minimum	if	for	any	

feasible	direction	d,	fሺx*ሻd	൒	0.	

Theorem	2.	If	fሺxሻ	is	twice	differentiable,	x*	is	a	relative	or	local	minimum	if	for	any	

direction	feasible	d,	fሺx*ሻd	൒	0;	and	if	fሺx*ሻd	ൌ	0,	dTHሺx*ሻd	൒	0.	

Theorem	3.	If	X	is	a	convex	set	and	fሺxሻ	is	convex	function	on	X,	a	relative	or	local	

minimum	is	a	global	minimum.		

Theorem	4.	If	X	is	a	convex	set	and	fሺxሻ	is	convex	function	on	X	and	twice	

differentiable,	Hሺx*ሻ	is	positive	semidefinite	on	X.		

Theorem	5.	If	A	is	symmetric,	all	eigenvalues	are	real.		

Theorem	6.	If	all	eigenvalues	are	positive	ሺnegativeሻ,	A	is	positive	ሺnegativeሻ	

definite.	

Theorem	7.	If	one	or	more	eigenvalues	are	0,	A	is	singular.	

Theorem	8.	If	A	is	symmetric,	V	is	the	eigenvector	matrix	and	E	is	diagሺeሻ	where	AV

ൌ	VE	and	V	can	be	chosen	orthonormal	ሺV‐1	ൌ	VTሻ,	thus	A	ൌ	VEVT.	

Iterative	Methods.	The	nonlinear	function	fሺxሻ	cannot	be	solved	directly	through	

factorization	methods.	The	algorithms	for	solving	unconstrained	nonlinear	

optimization	problems	can	be	broadly	defined	as	derivative‐free	methods,	methods	

using	first	derivatives	and	methods	using	second‐order	derivatives.	Since	the	

standard	ACOPF	is	continuous	and	differentiable,	we	focus	on	derivative‐based	

methods.		

	 Most	of	the	methods	described	are	iterative	methods	which	generate	a	

sequence	of	XK	ൌ	ሼx1,	x2,	...	xk,	...	xK	ሽ	for	K	iterations.	We	call	any	specific	process	a	

solver	or	an	algorithm.	Generally,	all	iterative	search	methods	have	a	five	step	

process:	

Step	1.	Choose	a	function	to	optimize;	set	k	ൌ	0.	Choose	an	initial	point:	x0.	

Step	2.	Choose	a	search	direction:	dk.	

Step	3.	Choose	step	size	sk		where	sk		is	a	positive	scalar	and	calculate	a	new	point:	

xk൅1	ൌ	xk	൅	sk	dk.	

Step	4.	Test	for	stopping:	If	xk൅1	satisfies	the	convergence	criteria	or	exceeds	the	

time	allotted,	then	set	K	ൌ	k൅1	and	stop.		Otherwise,	if	xk൅1	does	not	satisfy	the	

convergence	criteria,	then	set	k	ൌ	k	൅	1	and	go	to	Step	2.	

There	are	numerous	approaches	for	each	step.	Whereas	line	search	methods	

compute	a	search	direction	dk	and	then	a	step	size	sk	,	trust	region	methods	define	a	
Page 7

region	TR  X,	typically	an	elipse	around	the	current	iterate	xk,		to	choose	a	step	size	

sk		є	TR	and	direction	dk	simultaneously.	The	choice	of	the	step	size	and	then	

direction	in	line	search	methods,	or	in	the	case	of	trust	region	methods	the	choice	of	

the	trust	region	and	then	determining	sk	dk,	is	important	in	promoting	convergence	

from	remote	starting	points.	

	 In	Step	1,	the	user	can	supply	the	starting	point.	If	the	initial	point	is	

infeasible,	finding	a	feasible	solution	can	be	an	optimization	problem	in	itself.		For	

parallel	algorithms,	the	process	may	start	with	many	initial	points.	The	simplest	

parallel	algorithm	is	the	horse	race	where	the	approaches	are	started	in	parallel	and	

the	first	to	converge	terminates	the	algorithm.	More	complex	parallel	algorithms	

interact.	

	 Steps	2	and	3	are	often	defined	as	a	single	step	xk൅1	ൌ	Fሺxkሻ	where	F	is	the	

descent	algorithm.	Whether	a	linear	search	or	trust	region	method	is	applied,	the	

result	is	where		

fሺxk	൅	skdkሻ	൑	fሺxkሻ		

and	hopefully		

fሺxk	൅	skdkሻ	൏൏	fሺxkሻ.		

Typically	dk	is	a	descent	direction	where	

dk Tfሺxkሻ	൏	0.	

For	instance	we	can	determine	the	search	direction	dk	by	Newton‐Raphson’s	

method:	

2fሺxkሻdk	ൌ	‐fሺxkሻ.	

We	will	discuss	Netwon‐Raphson	in	further	detail	below.	Note	that	we	may	replace	

the	above	Hessian	Hሺxkሻ	ൌ	2fሺxkሻ	with	any	symmetric	and	nonsingular	matrix	Bk.	If	

Bk ൌ	I,	the	descent	algorithm	is	called	steepest	descent.	If	Bkൎ	Hሺxkሻ,	the	method	is	

called	quasi‐Newton.	For	example,	the	restricted	step	varies	from	a	Newton	step	to	a	

small	steepest	decent	step.	The	Newton	step	approach	has	quadratic	convergence,	

but	can	fail	to	converge.	The	steepest	decent	approach	has	slower	linear	

convergence,	but	fails	less	frequently.	

Step	4	tests	for	convergence.	Numerous	strategies	assume	Lipschitz	

continuity	of	the	gradient	ሺNocedal,	1999ሻ.	In	theory,	to	prove	convergence	an	

infinite	sequence	is	shown	to	converge	to	a	stationary	point.	A	subsequence	has	a	

limiting	property	that	satisfies	first‐order	necessary	or	second‐order	necessary	

conditions.	In	practice,	a	solver	stops	in	a	finite	number	of	iterations	when	the	

solution	is	‘close	enough’	to	the	optimal	solution.	For	example,	for	some	user	

defined	convergence	criteria	δ	൐	0,	stop	if	

|ሺxk	–	xk൅1ሻ|/|xk|	൑	δ	or	fk	–	fk൅1	൑	δ	or	|fሺxkሻ|	൏	δ.	
Page 8

Furthermore,	the	globally	convergent	algorithms	are	those	which	converge	from	

any	starting	point	to	a	local	optima,	and	therefore	have	the	property	that	the	

gradient	norms	converge	to	zero:	

	 limk∞	||fk||	ൌ	0	

Rate	of	Convergence.	The	error	term	is	defined	as		

εk	ൌ	xk	–	x*		

If	0	൑	β	ൌ	lim	sup	k є	∞	|εk൅1|/|εk|p	൏	∞,	β	is	called	the	convergence	ratio	and	p	is	the	

order	of	convergence.	If	p	ൌ	1	and	β	൐	0,	the	convergence	rate	is	called	linear.	If	p	ൌ	

1	and	β	ൌ	0,	the	convergence	rate	is	called	superlinear.	If	p	ൌ	2	and	β	൐	0,	the	

convergence	rate	is	called	quadratic.		

Feasible	Direction	Methods.	Feasible	direction	methods	require	that	for	a	given	

iterate	xk	є	X,	we	can	find	a	descent	direction	dk	which	is	also	a	feasible	direction	at	

xk.	Therefore	at	each	iteration	there	must	exist	a	new	feasible	point	of	the	form	xk	൅	

skdk	that	meets	the	reduction	condition.		

	 Feasible	direction	methods	include	both	derivative	free	and	derivative	based	

methods.	Line	search	without	derivatives,	such	as	the	Fibonacci	search,	only	

requires	function	evaluations.	The	Fibonacci	sequence	is	the	basis	for	this	approach	

in	which	sequential	points	are	chosen	such	that	the	discrepancy fሺxk൅1ሻ	െ	fሺxkሻ	is	

minimized.	However	if	f	is	differentiable,	we	solve	sfሺxk	൅	sdkሻ	ൌ	0.	This	is	a	

standard	algorithm	in	the	initial	stage	of	a	line	search.	

	 First‐derivative	gradient	approaches	use	the	first	derivative	as	the	feasible	

direction,	dk	ൌ	‐fሺxkሻ.	This	is	the	direction	at	which	f	decreases	most	quickly.	Below	

we	summarize	the	first‐derivative	methods	Gauss‐Seidel,	steepest	descent,	

conjugate	gradient,	and	Quasi‐Newton,	and	then	also	the	second‐derivative	based	

Newton’s	method.		

Gauss‐Seidel	Method.	The	Gauss‐Seidel	approach	chooses	the	unit	vector	as	a	

coordinate	direction	dk	ൌ	ሺ0,...	1k,	...0ሻ	and	a	line	search	in	one	variable	is	performed	

for	some	s	൐	0:	

	 mins fሺxk	൅	sdkሻ		

This	problem	is	sequentially	solved	in	each	coordinate	direction	and	repeated	until	

the	convergence	criteria	is	met.	This	approach	has	linear	convergence	rate.	

Steepest	Descent	Method.	The	steepest	descent	method,	which	is	fundamentally	

xk൅1	ൌ	Fሺxkሻ	ൌ	xk	–fሺxkሻ	
Page 9

for	unit	step	size	and	more	specifically,		

xk൅1	ൌ	Fሺxkሻ	ൌ	xk	൅	skdk

for	sk,	dk	൐	0	suffers	from	slow	convergence	due	to	zigzagging,	as	shown	in	

Figure	1.	Numerous	methods	have	been	proposed	to	avoid	zigzagging.		

Ideally	in	choosing	the	step	size	s,	we	find	the	local	minimizer	to	the	above	

objective	function	in	order	to	obtain	a	substantial	reduction	in	f.	However	this	may	

require	too	many	evaluations	of	f,	and	in	the	case	if	f		is	differentiable,	then	sf.	

Therefore	strategies	that	perform	inexact	line	search	identify	an	adequate	reduction	

in	f	by	trying	candidate	values	for	s.	As	we	shall	see,	exact	line	searches	are	not	

needed	and	may	not	be	the	best	approach	because	algorithms	that	achieve	rapid	

convergence	can	sometimes	conflict	with	the	global	convergence	requirements,	and	

vice	versa.	

Figure	1.	Zigzagging	with	steepest	descent	and	lack	of	zigzag	with	conjugate	

directions.	

Conjugate	Gradient	Method.	The	introduction	of	the	conjugate	gradient	method	by	

Fletcher	and	Reeves	in	1964	was	the	inception	of	large‐scale	nonlinear	optimization.	

The	conjugate	gradient	method	determines	conjugate	directions	to	the	Hessian	

(x0, y0)

Steepest Descent

Conjugate Gradient
Page

10

through	evaluating	f	and	f,	but	without	directly	evaluating	the	Hessian	2f.	The	

procedure	is	as	follows:	

Step	1.	Let	g0ൌ	fሺx0ሻ	and	d0	ൌ	‐g0;	set	k	ൌ	0.	

Step	2.	Solve	for	xk൅1	ൌ	xk	൅	skdk	where sk	ൌ	argmins fሺxk ൅	sdkሻ.	

Step	3.	Determine	direction	dk൅1ൌ	gk൅1൅	rkdk	where		

rkൌ	gk൅1Tgk൅1/gk

Tgk,,	and	

	 gk൅1ൌ	fሺxk൅1ሻ.	

Step	4.	Go	to	Step	2	and	repeat	for	k	ൌ	1	...,	nെ1.	

	

Quasi‐Newton	Method.	

	 In	the	late	1950’s,	Davidon	introduced	quasi‐Newton	methods,	based	on	

Newton’s	method	and	also	known	as	secant	methods,	which	use	a	sequence	of	

positive	definite	matrices	to	approximate	the	Hessian	ሺor	inverse	Hessianሻ.	These	

methods	are	inexact	line	search	methods	and	have	a	superlinear	convergence	rate.	

Fletcher	and	Powell	demonstrated	that	Davidon’s	approach	is	equivalent	to	

conjugate	gradient	method	when	applied	with	exact	line	searches	to	convex	

quadratic	functions.	Quasi‐Newton	methods	only	require	first	derivatives,	

approximate	the	Hessian	and	adhere	to	the	descent	property	due	to	the	positive	

definiteness.	

	 For	example,	if	fሺxሻ	is	twice	differentiable,	the	second‐order	Taylor’s	series	

expansion	of	fሺxሻ	at	xk	is:	

fሺxሻ	ൌ	fሺxkሻ	൅	fሺxkሻTሺx	െ	xkሻ	൅	ሺx	െ	xkሻTHሺxkሻሺx	െ	xkሻ/2	൅	Oሺ|x	െ	xkሻ|2ሻ		

Dropping	Oሺ|x	െ	xkሻ|2ሻ,	fሺxሻ	is	a	quadratic	approximation	of	fሺxሻ	

fሺxሻ	ൌ	fሺxkሻ	൅	fሺxkሻTሺx	െ	xkሻ	൅	ሺx	െ	xkሻTHሺxkሻሺx	െ	xkሻ/2		

If	fሺxሻ	is	quadratic,	then	

fሺxሻ	ൌ	bTx	൅	xTQx/2,		

fሺxሻ	ൌ	bT	൅	xTQ,		

Hሺxሻ	ൌ	Q,	and		

	 Oሺ|x	െ	xkሻ|2ሻ	ൌ	0.		

The	above	expansion	is	exact	and	x*	is	a	solution	to		

Qx*	ൌ	െb.		

If	Q	is	positive	semidefinite	at	x*,	we	have	that	

fሺxሻ	ൌ	0,		

x*TQx*	൒	0,	and		
Page

11

fሺxሻ	൒	fሺx*ሻ.		

Therefore	in	quasi‐Newton	methods	the	Hessian	approximation	Hk	is	chosen	to	

satisfy	

fሺxk	൅	Δxሻ	ൌ	fሺxkሻ	൅	HkΔx

which	is	called	the	secant	equation	and	is	the	Taylor	series	of	the	gradient	itself.	

Steps	in	the	Quasi‐Newton	Method:	

Step	1.	An	approximate	initial	value	of	H0	ൌ	I	is	frequently	sufficient	for	

convergence.	Set	k	ൌ	0.	

Step	2.	Hk	is	positive	definite.	Let	dk	ൌ	െ	Hk

‐1fሺxkሻ	and	compute	sk	ൌ	argmins fሺxk ൅	

sdkሻ.	Set	xk൅1	ൌ	xk	൅	skdk,	where	sk satisfies	Wolfe	conditions:	

fሺxk ൅	skdkሻ	൑	fሺxkሻ	൅	α1sk	dk Tfሺxkሻ	

	 dk Tfሺxk ൅	skdkሻ൒	α2skdk Tfሺxkሻ	

where	0	൏	α2	൏	α2	൏	1.	The	above	first	inequality	is	sometimes	referred	to	as	the

Armijo	rule,	which	ensures	that	the	step	length	sk	results	in	a	sufficient	decrease	in	f,	

and	the	second	is	the	curvature	condition,	which	ensures	that	the	slope	has	been	

reduced	sufficiently	ሺNocedal,	1999ሻ.	Another	way	of	stating	the	curvature	

condition	is	as	follows:		

vk

	Tሺfሺxk ൅	skdkሻെ	fሺxkሻሻ		ൌ	vk	THkvk൐	0	

where	vk	ൌ	skdk.	The	Wolfe	conditions	can	be	applied	in	most	line	search	methods	

and	are	important	in	implementing	quasi‐Newton	approaches.	

Step	3.	Compute	Hk൅1	by	a	rank‐two‐update,	where	the	inverse	of	the	Hessian	Bk	ൌ	

Hk

‐1	is	updated	by	the	sum	of	two	symmetric	rank	1	matrices:			

Bk൅1	ൌ	Bk	൅	βkvkvk

Tെ	γk	HkykykTHk

where	vk	ൌ	skdk,	yk	ൌ	fሺxk ൅	skdkሻെ	fሺxkሻ,	βkൌ	1/vk

Tyk,	and	γk	ൌ	1/ykTHkyk.	The	

Hessian	approximation	is	chosen	to	satisfy	the	quasi‐Newton	condition:	

Hk൅1ሺfሺxk൅1ሻെ	fሺxkሻሻ	ൌ	xk൅1	െ	xk.	

Set	k	ൌ	k	൅	1;	go	to	Step	2.	

The	above	algorithm	is	known	as	the	Davidon‐Fletcher‐Powell	ሺDFPሻ;	the	Broyden‐

Fletcher‐Goldfarb‐Shanno	ሺBFGSሻ	method	is	the	dual	of	DFP	because	whereas	the	

DFP	converges	to	the	inverse	of	the	Hessian,	the	BFGS	method	converges	to	the	

Hessian	in	itself	and	therefore	is	a	more	direct	approach.	 	
Page

12

	 In	1970,		Broyden,	Fletcher,	Goldfarb,	and	Shanno	independently	introduced	

a	secant	approach	now	known	as	the	Broyden‐Fletcher‐Goldfarb‐Shanno	ሺBFGSሻ	

method.		In	this	direct	approach,	

Hk൅1	ൌ	Hk	൅	βkykyk

Tെ	γk	HkvkvkTHk

where	instead	γk	ൌ	1/vk

THkvk.	

	 The	Levenberg‐Marquardt	method	uses	a	trust‐region	to	solve	a	system.	The	

method	determines	a	value	ν	൒	0	such	that	Hሺxkሻ	൅	νI	is	positive	definite	and	then	

solves	

sk∆xk	ൌ	െ	ሺHሺxkሻ	൅	νIሻ‐1fሺxkሻ.	

An	iteration	radius	rk	may	impose	a	limit	on	the	step	size	sk,	or	alternatively	ν

controls	the	length.	However	by	not	controlling	the	length	rk,	degeneracy	ሺi.e.	non‐

invertible	ሺHሺxkሻ	൅	νIሻሻ	and	ill‐conditioning	issues	can	result	in	a	step	size	and	

direction	that	are	badly	determined	by	ν.	

Newton’s	Method.	Newton’s	method	has	been	applied	with	both	line	search	and	

trust	region	approaches.	Non‐monotone	methods	in	which	the	function	values	are	

allowed	to	increase	at	some	iterations	have	been	useful	in	solving	highly	nonconvex	

problems.		

The	basic	Newton’s	algorithm	is:		

Step	1.	Guess	x0,	and	set	k	ൌ	0.	

Step	2.	At	xk,	evaluate	fሺxkሻ.	

Step	3.	Evaluate	Bk	ൌ	Hሺxkሻ	or	an	approximation	Bk	ൌ	Hሺxkሻ.	For	example,	when	

Hሺxkሻ	is	not	positive	definite,	then	modifications	could	be	applied	to	determine	

Hሺxkሻ	and	guarantee	a	decrease	in	fሺxkሻ.	Such	methods	may	modify	Hሺxkሻ	so	that		

fሺxkሻTBkfሺxkሻ	൐	0,		

or	such	methods	may	compute	a	negative	search	direction	which	satisfies	

		 ∆xk

TBk	∆xk	൐	0	and		

xfሺxkሻT∆xk	൑	0.		

Step	4.	Solve:	Bkdk	ൌ	െfሺxkሻ.	

Step	5.	If	the	convergence	error	is	less	than	the	tolerance,	ሺe.g.,	||fሺxkሻ	||	൑	ε	and	

||dk||	൑	εሻ	then	stop;	else	continue.	

Step	6.	Find	sk	so	that	0	൏	sk	൑	1	and	fሺxk	൅	skdkሻ	൏	fሺxkሻ.		

Step	7.	xk൅1	ൌ	xk	൅	skdk.	Set	k	ൌ	k	൅	1	and	go	to	Step	2.	
Page

13

Note	that	recalculating	fሺxkሻ,	Bk,	and	skdkൌ	xk൅1	െ	xk	at	each	iteration	may	be	very	

costly	in	terms	of	computational	time.	Let	dk	ൌ	fሺxkሻ,	xk൅1	ൌ	xk	൅sdk.	If	Hሺxkሻ	is	

positive	semidefinite,		

u	ൌ	dk

THሺxkሻdk	൒	0	and	v	ൌ	dkTdk൒	0	

fሺxk൅1ሻ	െ	fሺxkሻ	ൌ	dk

Tsdk	൅	sdkTHሺxkሻsdkൌ	vs	൅	us2

for	some	s	൏	0,	where	we	have	fሺxk൅1ሻ	൏	fሺxkሻ.		

For	globally	convergent	Newton	and	quasi‐Newton	methods,	Bk	must	be	positive	

definite	and	also	have	a	bounded	condition	number	because	a	high	condition	

number	can	result	in	ill‐conditioned	problems	ሺNocedal,	1999ሻ.	Intuitively,	the	

search	direction	dk	is	not	near	orthogonal	to	the	gradient	in	order	to	promote	

convergence.	However	for	problems	that	are	ill‐conditioned,	it	may	be	necessary	to	

search	along	directions	that	are	nearly	orthogonal	to	the	gradient.	

Numerical	Efficiency	and	Stability.	Applying	the	appropriate	techniques	results	in	

computationally	efficient	and	numerically	stable	procedures.	For	example	the	core	

of	solving	a	system	of	linear	algebraic	equations,	Ax	ൌ	b,	is	decomposing	or	factoring	

the	coefficient	matrix.	Through	this	process,	the	solution	can	be	obtained	with	much	

less	computational	effort,	and	better	methods	are	solved	faster	with	less	numerical	

error.	Matrix	inversion	is	usually	unstable	and	results	in	dense	matrices	and	

unnecessary	floating‐point	calculations.	Sometimes	for	simplicity	and	sometimes	

out	of	ignorance,	an	algorithm	is	presented	as	inverting	a	matrix	at	each	iteration:	x

ൌ	A‐1b.	This	approach	is	not	only	computationally	expensive,	but	also	numerically	

unstable.			

	 For	large	problems,	A	is	usually	sparse,	but	the	inverse	is	usually	dense.	

Sparse	matrices	are	stored	efficiently	by	storing	only	nonzeroes	and	their	row	and	

column	index.	When	one	of	the	operands	in	a	floating‐point	operation	is	zero,	it	can	

be	skipped,	saving	time.	Problems	are	referred	to	as	ill‐conditioned	when	small	

perturbations	to	the	coefficients	result	in	an	unpredictably	large	change	of	the	

solution.	Certain	conditions	such	as	positive	definiteness	can	be	checked,	and	

corrections	to	an	indefinite	matrix	can	be	applied.			

	 Accumulating	floating	point	round‐off	and	also	cancellation	errors	cause	

search	vectors	to	lose	orthogonality	in	conjugate	search	methods.	Similarly,	steepest	

descent	methods	can	accumulate	floating	point	round‐off,	which	may	cause	xk	to	

converge	to	some	point	near	x.	This	effect	can	be	deterred	by	occasionally	

recomputing	the	correct	residual.	
Page

14

Scaling.	Numerical	stability	starts	with	scaling	the	problem.	Scaling	is	part	theory	

and	part	art	form.	Proper	consistent	scaling	means	primal,	dual	solution	values,	

constants	and	derivatives	of	nonlinear	terms	ሺJacobian	elementsሻ	in	absolute	value	

be	around	1,	e.g.	from	0.01	to	100.	Most	solvers	have	automatic	scaling	options.		

Factoring	Linear	Equations.	In	practice,	at	each	iteration	dk	ൌ	െ	Bkfሺxkሻ	

is	solved	by	sparse	matrix	factorization	techniques	often	a	matrix	is	represented	as	

the	product	or	sum	of	matrices;	such	a	representation	is	termed	as	a	decomposition	

or	factorization	of	the	original	matrix.	The	Gaussian	elimination	product	form	of	the	

inverse	of	A	stores	the	inverse	as	a	product	of	elementary	column	matrices,	where	

Aെ1	ൌ	Ek	Ekെ1	൉	൉	൉	E2E1.		

and	E	is	an	elementary	column	operation	where	all	but	one	column	is	the	identity	

matrix.	Therefore	we	have	that			

x	ൌ	Aെ1 b	ൌ	Ek	Ekെ1	൉	൉	൉	E2E1b.	

For	iterative	procedures,	as	k	become	large,	numerical	errors	build‐up	and	

excessive	floating‐point	operations	occur.	Therefore	the	matrix	needs	to	be	

refactored	into	k	ൌ	n	factors.	

Another	technique	is	LU	factorization,	where	L	is	a	lower	triangular	matrix	

and	U	is	an	upper	triangular	matrix.	For	any	matrix	A,		

Ax	ൌ	LUx	ൌ	LሺUxሻ	ൌ	b		

We	can	solve	Ly	ൌ	b	by	substitution	and	then	Ux	ൌ	y	by	substitution.	For	additional	

details,	see	Markowitz	LU	decomposition	ሺReid,	1976;	1982ሻ	and	Bartels‐Golub	

update	ሺBartels	and	Golub,	1969;	1971ሻ.	

Factoring	Linear	Equations	with	Symmetric	Matrices.	Symmetric	matrices	have	more	

efficient	factorizations.	The	Cholesky	decomposition	of	a	Hermitian,	positive‐

definite	matrix,	A,	is	the	product	of	a	lower	triangular	matrix	and	its	conjugate

transpose.	A	ൌ	LL*	where	L*	is	the	conjugate	transpose.	Cholesky	decomposition	is	

roughly	twice	as	efficient	as	the	LU	decomposition	for	solving	systems	of	linear	

equations.	With	QR	factorization,	A	results	in	an	orthogonal	matrix	Q	that	is	a	basis	

for	the	column	space	of	A	and	an	upper	triangular	matrix	R	where	AൌQR.	These	

approaches	lose	certain	invariance	properties	associated	with	Newton’s	method	and	

Newton’s	method	with	line	search,	and	therefore	may	not	transform	correctly,	when	

the	Hessian	is	ill‐conditioned.		

4.	Constrained	Nonlinear	Optimization	

	

The	constrained	nonlinear	problem	or	program	in	canonical	primal	ሺPሻ	form	is:	
Page

15

f*ൌ	inf	ሼfሺxሻ	|	gሺxሻ	൑	0,	hሺxሻ	ൌ0,	xXሽ			ሺPሻ	

This	formulation	explicitly	introduces	an	additional	set	of	inequality,	gሺxሻ	൑	0,	and	

equality	constraints,	hሺxሻ	ൌ	0,	that	are	linear	or	nonlinear.	If	xX,	hሺxሻ	ൌ	0	and	gሺxሻ	

൑	0,	x	is	said	to	be	a	feasible	solution.	Nonlinear	equality	constraints	can	be	

represented	as	gሺxሻ	൑	0	and	‐gሺxሻ	൑	0.	We	define	cሺxሻ	ൌ	0	as	ሺgሺxሻ൅sൌ0,	hሺxሻൌ0ሻ.	

Lagrangian	Function.	The	algorithms	for	solving	constrained	nonlinear	optimization	

problems	typically	start	with	a	transformation	of	the	Lagrangian.	The	Lagrangian	

dual	ሺDሻ	formulation	is		

L*	ൌ	supλ	൒	0ሼinf	Lሺx,	λሻ|	x	X,	hሺxሻ	ൌ	0ሽ	

where	Lሺx,	λሻ	ൌ	fሺxሻ	൅	λTgሺxሻ	is	the	Lagrangian	function	and λ	is	the	Lagrangian	

multiplier	vector.	The	variable	λ	is	sometimes	called	the	dual	or	slack	variable.	Let	

ሺx*,	λ*ሻ	be	the	solution	to	ሺDሻ.		If	Lሺx,	λሻ	is	unbounded	for	x	X,	then	L* ൌ	∞.	

Furthermore,	if	the	feasible	region	to	ሺDሻ	is	an	empty	set,	then	L* ൌ	‐∞.	

	 Iterative	methods	such	as	penalty	and	augmented	Lagrangian	methods,	

barrier	or	interior	point	methods,	variable	metric	methods,	sequential	linear	

programming	methods	and	sequential	quadratic	programming	methods	often	solve	

a	sequence	of	generalized	Lagrangian	functions.	The	generalized	Lagrangian	

function	is:	

	 Lሺx,	λሻ	ൌ	fሺxሻ	൅	Λሺgሺxሻ,	λ,	μሻ	

that	includes	a	positive	penalty	parameter	μ	൐	0	with	a	penalty	function	Λ.	

Karush‐Kuhn‐Tucker	ሺKKTሻ	Conditions.	If	fሺxሻ	and	gሺxሻ	are	differentiable	the	KKT	

necessary	ሺfirst‐orderሻ	conditions	are						

fሺxሻ	൅	λTgሺxሻ	ൌ	0,			

	 	λ	൒	0,				

λTgሺxሻ	ൌ	0,		and	

hሺxሻ	ൌ	0.	

Under	certain	constraint	qualifications,	for	example,	if	there	is	gሺxሻ	൏	0	and	hሺxሻ	ൌ	

0,	or	ሾg ሺxሻ,	hሺxሻሿ	has	full	rank	where	the	active	set	AS	denotes	an	active	

constraint,	the	KKT	are	necessary	for	a	local	optimum.	If	x	is	a	local	optimal	solution	

to	P	and	also	satisfies	the	KKT	sufficiency	conditions	with	strict	complementarity	to	

λ,	then	ሺx,	λሻ	 X	ൈ	Rm൅	is	a	saddle	point	of	the	Lagrangian:	

Lሺx,	λሻ	൑	Lሺx,	λሻ	൑	Lሺx,	λሻ	

for	all	ሺx,	λሻ	 X	ൈ	Rm൅.	Therefore,	the	x	and	λ	are	optimal	solutions	to	the	primal	

and	dual	problems,	respectively	and	
Page

16

infሼfሺxሻ	|	gሺxሻ	൑	0,	hሺxሻ	ൌ	0,	xXሽ	൒	supλ	൒	0	inf	ሼLሺx,	λሻ	|	hሺxሻ	ൌ	0,	xXሽ.				

This	property	is	known	as	weak	duality	when	x	is	feasible	in	P	and	λ	is	feasible	in	D.	

f*	‐	L*	൒	0	is	called	the	duality	ሺcomplementarityሻ	gap.	Strong	duality	results	in	a	

zero	duality	gap:		

f*	‐	L*	ൌ	0.	

	 Furthermore,	if	x	is	feasible	and	there	exists	λ	such	that	strong	duality	holds	

under	certain	constraint	qualifications	where	x	is	a	local	optima.	Certain	stability	

properties	lead	to	necessary	and	sufficient	conditions	which	hold	with	equality	for	

the	existence	of	a	global	saddle	point	of	the	Lagrangian	function	with	respect	to	X	ൈ	

Rm൅.	If	gሺxሻ	and	fሺxሻ	are	convex	functions,	then	the	KKT	are	sufficient	for	global	

optimality.	If	fሺxሻ	and	gሺxሻ	are	twice	differentiable,	and	the	Hessian	

Hሺx,	λሻ	ൌ	2Lሺx,	λሻ	ൌ	2fሺxሻ	൅		λT2gሺxሻ		

is	positive	definite,		then	we	have	a	local	minima,	or	if	Hሺx,λሻ	is	negative	definite,	

then	we	have	a	local	maxima.	The	result	is	analogous	to	strong	duality.		

	 Often	the	KKT	may	not	hold	at	the	optimal	solution	if	the	problem	is	

nonconvex,	and	therefore	such	KKT	necessary	conditions	are	not	sufficient	to	prove	

global	optimality.	When	the	primal	problem	is	nonconvex,	there	may	be	local	optima	

that	are	not	globally	optimal.		

	 The	following	subsections	detail	the	dominant	approaches,	which	are	based	

directly	on	solving	the	KKT	conditions	for	constrained	optimization.	

Augmented	Lagrangian	and	Penalty	Methods.	In	1969	in	separate	papers,	Hestenes	

and	Powell	introduced	the	augmented	Lagrangian.	In	the	1970s,	it	gained	a	strong	

following	and	is	still	applied	today.	The	augmented	Lagrangian	is	the	Lagrangian	

function	with	a	quadratic	penalty	term:		

Λሺgሺxሻ,	λ,	μሻ	ൌ	λTgሺxሻ	൅	μTgሺxሻ2					

The	basic	algorithm	for	the	augmented	Lagrangian	and	general	penalty	function	

methods:		

Step	1.	Choose	λ0,	μ0	in	iteration	k	ൌ	0.	

Step	2.	Find	xk൅1	ൌ	argminx	ሼfሺxkሻ	൅	λk

Tgሺxkሻ	൅	μkTgሺxkሻሽ.	

Step	3.	Update	μk൅1൐	μk		and	λk൅1	ൌ	λk	൅	μk

Tgሺxkሻ.	

Step	4.	If	convergence	criteria	is	met,	then	stop;	else	go	to	Step	2.	

Barrier/Interior‐Point	Methods.	Barrier	methods	require	an	interior	point	x0,	where	

gሺx0ሻ	൏	0,	and	construct	a	sequence	of	unconstrained	problems	with	μ	൐	0	
Page

17

controlling	the	weight	of	the	barriers.	Barrier	terms	cannot	be	defined	for	equality	

constraints	because	μ	increases	tends	to	zero	as	the	algorithm	approaches	the	

boundary	to	the	feasible	region.	For	example,	the	most	popular	barrier	functions	are	

of	the	form		

Λሺgሺxሻ,	λ,	μሻ	ൌ	μTlogሺെgሺxሻሻ	and	Λሺgሺxሻ,	λ,	μሻ	ൌ	െ1/μTgሺxሻ.	

	 Originally	proposed	by	Frisch	ሺ1955ሻ	and	later	developed	by	Fiacco	and	

McCormick	ሺ1968ሻ,	the	logarithmic	barrier	method	was	applied	extensively	in	the	

1960s	and	1970s	and	has	become	foundational	for	nonlinear	interior‐point	

methods.	In	1984,	interior	point	methods	gained	momentum	with	Karmarkar’s	

polynomial‐time	LP.	in	the	log‐barrier	approach,	steps	are	generated	by	iteratively	

solving	an	unconstrained	optimization	problem	similar	to	the	form:	

	 min		 fሺxሻ	െ	μk

Tlogሺെgሺxሻሻ			

where		

μk	൐	0,	and	

	 ∂ሺμk

Tlogሺെgiሺxkሻሻ	ሻ/∂xk	ൌ	μkTgiሺxkሻ/giሺxkሻ.	

As	k	→	∞,		we	have	μk	→	0	and	μk/giሺxkሻ	→	λi		for	all	i.	

	 Since	the	logarithm	of	a	small	number	is	negative,	the	െμklogሺെgiሺxkሻሻ	term	

is	positive	for	any	variable	near	the	bound;	this	effectively	creates	barriers	pushing	

the	variables	into	the	interior	of	the	feasible	region.	Barrier	methods	suffer	from	ill‐

conditioning	in	the	Hessian	as	the	optimum	is	approached.	More	generally,	interior	

point	methods	are	often	sensitive	to	the	initialization	and	reduction	of	the	barrier	

ሺor	centeringሻ	parameter	μ.	Such	barrier	methods	are	also	sensitive	to	the	step	size	

when	a	global	convergence	strategy	is	employed.	However	applying	factorization	

and	global	convergence	strategies	specifically	designed	for	the	given	barrier	

function	can	address	some	of	the	above	complications.		

In	1992,	Polyak	introduced	a	modified	barrier	method	that	scales	the	barrier	

function	so	that	initial	infeasibility	is	not	a	problem.	Also	in	1992,	Mehrorta	

introduced	a	predictor‐corrector	mechanism	to	promote	convergence	by	reducing	

the	number	of	matrix	factorizations	when	determining	search	directions.	The	

method	uses	the	Cholesky	decomposition	to	find	two	different	directions:	a	

predictor	and	a	corrector.		The	predictor	computes	an	optimizing	search	direction	

based	on	a	first	order	term.	The	corrector	step	uses	the	same	Cholesky	

decomposition	found	during	the	predictor	step.	The	search	direction	is	the	sum	of	

the	predictor	direction	and	the	corrector	direction.		By	determining	the	centering	

parameter	adaptively	ሺrather	than	prior	to	the	affine‐scaling	directionሻ	and	enabling	

corrections	in	the	predicted	direction,	Mehrotra’s	approach	introduces	extra	

computation	per	iteration	but	significantly	reduces	the	number	of	iterations	
Page

18

required.	In	1996,	Gondzio	introduced	multiple	centrality	corrections	strategy	that	

maximizes	the	possible	step	size	increase	in	Mehrotra’s	approach	by	computing	up	

to	several	corrections,	which	leads	to	an	even	larger	decrease	in	the	residuals.	Most	

interior	point	approaches	have	leveraged	the	predictor‐corrector	method.	Although	

there	is	no	theoretical	complexity	bound,	the	method	is	widely	used	in	practice.	It	

appears	to	be	computationally	efficient	and	to	converge	very	fast	when	close	to	the	

optimum.		

Primal‐Dual	Methods.	Primal‐dual	approaches	apply	ሺquasi‐ሻNewton’s	method,	and	

similar	to	other	iterative	methods	we’ve	covered	thus	far,	includes	a	search	step	and	

direction,	and	a	convergence	criteria.		

Consider	the	primal	program	

	 minx,y fሺxሻ		

	 subject	to		

gሺxሻ	൅	y	ൌ	0	

y	൒	0		

with	the	Lagrangian		

Lሺx,	s,	λ,	μሻ	ൌ	fሺxሻ	൅	λTሺgሺxሻ	൅	yሻ	െ	μTy.		

The	first	order	conditions	are	

fሺxሻ	൅	λTgሺxሻ	ൌ	0	,	

μ	െ	ν		ൌ	0,	

gሺxሻ	൅	y		ൌ	0	,	

y	൒	0	,		

μ	൒	0,	

νiyi		–	1/t	ൌ	0,	and		

νy	ൌ	0		

where		is	element‐by‐element	multiplication	and	t	→	0.		

The	optimization	problem	is	

	 min	νy		

	 subject	to		

fሺxሻ	െ	λTgሺxሻ	ൌ	0		

μ	െ	ν	ൌ	0		

gሺxሻ	െ	y	ൌ	0		

y	൒	0			

μ	൒	0.		
Page

19

Therefore,	x,	s,	λ,	μ	and	t	are	updated	after	each	Newton	step.	Primal‐dual	interior	

point	methods	allow	for	primal	and	dual	iterates	to	be	infeasible,	e.g.	gሺxሻ	൅	s	്	0.	

Sequential	Linear	Programming	ሺSLPሻ	Methods.	The	SLP	method	determines	the	

solution	to	a	nonlinear	problem	through	solving	a	sequence	of	linear	

approximations,	mostly	using	first‐order	Taylor	series	expansions.	For	each	

iteration	k,	we	solve	the	linear	program	

	 xk൅1	ൌ	argmin	ሼfሺxkሻTd	| gሺxkሻTd	൅	gሺxkሻ	൒	0,	dl	൑	d	൑	duሽ	

where	fሺxkሻ can	be	replaced	with	linearized	generalized	Lagrangians.	Since	the	

linearizations	are	not	necessarily	bounded,	trust	regions	where xk൅1	є	TRk	can	be	

used	to	ensure	convergence.	Furthermore,	this	approach	may	be	used	in	

conjunction	with	a	penalty	or	merit	function	and	step	restrictions.	If	the	problem	is	

convex,	the	linear	constraints	are	always	outside	the	feasible	region	to	the	nonlinear	

problem,	but	are	guaranteed	to	achieve	convergence	as	k	→	0.	

	 The	SLP	can	handle	very	large	problems,	benefits	from	using	LP	solvers,	but	

may	converge	slowly	and	may	violate	nonlinear	constraints.	SLP	methods	are	

suitable	for	solving	large‐scale	nonlinear	programming	problems,	since	large	linear	

programs	can	be	solved	efficiently.	Because	of	the	linearity	of	the	approximation,	

special	scarcity	patterns	of	the	Jacobian	matrix	gሺxkሻ	are	passed	to	the	constraint	

matrix	of	the	linear	program	directly.	A	particular	advantage	of	sequential	linear	

programming	methods	is	that	although	second	order	information	is	not	used,	

convergence	is	linear	even	in	the	case	of	highly	nonconvex	problems.	

Sequential	Quadratic	Programming	ሺSQPሻ	Methods.This	method	is	also	referred	to	as	

a	projected	Lagrangian	or	Lagrange‐Newton	approach	because	the	accompanying	

subproblem	minimizes	a	quadratic	approximation	to	the	Lagrange	function	subject	

to	a	linearized	constraint	set.	

	 SQP	methods	solve	a	sequence	of	quadratic	subproblems	to	determine	the	

active‐set	and	a	search	direction	d	as	the	solution	for	each	iteration.	The	key	idea	is	

to	approximate	second	order	information	to	get	a	fast	final	convergence speed.	Thus,	

we	define	a	quadratic	approximation	of	the	generalized	Lagrangian	function	Lሺx,	λ,	

μሻ	and	an approximation	of	the	Hessian	matrix	by	a	quasi‐Newton	Newton	where	Bk		

ൌ	2xLሺxk,	λk,	μkሻ.	Then	we	obtain	the	subproblem:	

	 min	dk

TBkdk/2	൅	fሺxkሻTdk

	 subject	to	

gሺxkሻTdk	൅	gሺxkሻ	ൌ	0			

dk	TRk
Page

20

	 The	KKT	conditions	are	applied	to	the	objective	of	the	subproblem	to	solve	

for	dk.	A	merit	function	is	typically	the	sum	of	the	objective	function	and	the	amount	

of	infeasibility	of	the	constraints.		

	 A	SQP	usually	requires	fewer	functions	and	gradient	evaluations	than	a	SLP,	

but	may	violate	nonlinear	constraints	and	is	harder	to	solve	than	a	SLP.	Another	SQP	

based	approach,	the	sequential	linear‐quadratic	program	ሺSLQPሻ,	decouples	the	

active‐set	identification	and	the	step	computation	in	the	QP	subproblem.	SLQP	omits	

the	quadratic	term	in	the	objective	function	and	solves	the	subproblem	for	a	trust	

region.	With	only	the	active‐set	of	constraints,	a	subsequent	linear	program	is	

solved	to	determine	dk.	

Generalized	Reduced	Gradient	Methods.	In	1961,	Rosen	introduced	a	gradient	

projection	method.	In	1969,	Abadie	and	Carpentier	introduced	the	generalized	

reduced	gradient	ሺGRGሻ	for	nonlinear	constraints	combining	quasi‐Newton	methods	

and	Wolfe’s	reduced	gradient	method.	The	actual	implementation	has	many	

modifications	to	make	it	efficient	for	large	models	see	in	Drud	ሺ1985	and	1992ሻ.	

There	are	several	GRG	variations	including	choices	of	Lagrangians	and	line	searches.	

We	start	with	the	problem:	

min	fሺxሻ	

subject	to	

gሺxሻ	ൌ	0	

0	൑	x	൑	xu		

At	major	iteration	k,	we	solve	the	following	nonlinear	program:	

GRGk:	min	Fሺx,	λk,	μkሻ	 	 	

subject	to	

gk	ൌ	gሺxkሻ	൅	gሺxkሻTሺx	‐	xkሻ	ൌ	0		 	 ሺdual	variable	λk൅1ሻ

where		

Fሺx,	λk,	μkሻ	ൌ	fሺxሻ	–	λk

Tሺgሺxሻ	–	gkሻ	൅	μkሺgሺxሻ	–	gkሻTሺgሺxሻ	–	gkሻ,	and		

μk	is	a	scalar	penalty	parameter.		

For	the	set	of	basic	ሺdependentሻ	variables	B,	the	set	of	superbasic	ሺindependentሻ	

variables	S,	and	the	set	of	nonbasic		ሺdependent,	fixed	at	a	boundሻ	variables	N	we	

can	partition	gሺxkሻ	as		

gሺxkሻ	ൌ	ሾB,	S,	Nሿ	

where	B	and	S	are	nonsingular,	and	we	can	partition	xk	as		

xB	ൌ	ሺxkሻB,	
Page

21

xS	ൌ	ሺxkሻS,	and	

		 xN	ൌ	ሺxkሻN

where		

BTxB	൅	STxS	൅	NTxN	ൌ	0.	

At	a	solution,	the	basic	and	superbasic	variables	are	between	their	bounds,	0	൑	xB	൑	

xu B	and	0	൑	xS	൑	xu S	while	nonbasic	variables	will	normally	be	equal	to	their	lower	or	

upper	bound.	At	a	solution,	S	will	be	no	more	than	the	number	of	nonlinear	

variables	and	xs	is	regarded	as	a	set	of	independent	variables	that	are	allowed	to	

move	in	a	direction	that	will	decrease	the	objective	function	value.	The	basic	

variables	are	easily	adjusted	to	satisfy	the	linear	constraints.	

If	no	improvement	can	be	made	with	the	current	definition	of	B,	S	and	N,	

some	of	the	nonbasic	variables	are	selected	to	be	added	to	S,	and	the	process	is	

repeated	with	an	increased	value	of	S.	If	a	basic	or	superbasic	variable	encounters	

one	of	its	bounds,	the	variable	is	made	nonbasic	and	the	value	of	S	is	reduced	by	one.	

A	step	of	the	reduced‐gradient	method	is	called	a	minor	iteration.	For	linear	

problems,	simplex	method	is	the	reduced‐gradient	method	with	the	number	of	

superbasic	variable	oscillating	between	0	and	1.	

The	steps	in	a	GRG	algorithm	are:	

Step	1.	Set	k	ൌ	0	and	start	with	a	feasible	solution,	x0.	

Step	2.	Compute	gሺxkሻ	ൌ	J	ൌ	ሾB,	S,	Nሿ,	fሺxkሻ,	and	the	reduced	gradient,	gr.	Select	

the	set	of	superbasic	variables,	xs,	as	a	subset	of	the	nonbasic	variables	that	can	be	

gainfully	changed.	For	the	superbasic	variables		

grሺxSሻ	ൌ	fሺxkሻTെJ	Tλk		

and	for	the	basic	variables		

grሺxBሻ	ൌ	0	ሺBTλ	ൌ	∂f/∂xb ሻ.	

Step	3.	Find	a	search	direction,	ds,	for	the	superbasic	variables	that	is	based	on	gr

and	possibly	some	second	order	information.	The	reduced	Hessian	of	fሺxሻ	may	be	

approximated	by	solving	a	system	of	the	form		

RTRq	൅	Zfሺxሻ	ൌ	0		

where		

Z	ൌ	ሾ‐B‐1S,	I,	0ሿ,			

Zfሺxሻ	is	the	reduced	gradient,	and		

RTR	ൎ	ZTHZ.	
Page

22

The	dense	upper	triangular	matrix	R	is	updated	in	various	approaches	in	order	to	

approximate	2fሺxሻ.	Therefore,	ds	ൌ	Zq	where	gr	is	applied	to	solve	the	line	search	

min0൏α൏β Fሺx	൅αds,	λk,	μkሻ.	

Step	4.	Solve	the	subproblem	GRGk;	let	ሺxk൅1,	λk൅1 μk൅1ሻ	ൌ	argmin	ሼGRGkሽ.		

Step	5.	If	gr	projected	on	the	bounds	is	smaller	than	the	convergence	criteria,	then	

stop;	if	not,	increment	k	ൌ	k	൅1	go	to	Step	2.	

	 Generalized	reduced	gradient	methods	can	be	extended	easily	to	very	large	

problems	and	problems	with	special	structure.	GRG	is	probably	more	robust	than	

SLP	and	SQP	methods,	and	once	a	feasible	solution	is	found	in	GRG	methods,	it	

remains	feasible.	Moreover,	GRG	approaches	are	related	to	SQP	methods	and	

therefore	there	exists	combinations	of	both	approaches.	

Conic	and	Semidefinite	Programming	Methods.	These	methods	are	forms	of	

quadratic	programming,	and	are	similar	to	linear	programming	because	they	are	

convex	programs.	For	linear	program	

minx cTx

subject	to		

												 Ax	=	b	

	 x	≥	0	

the	equivalent	conic	program	is	

minx cTx

subject	to		

												 ||Aix	+	bi||2	≤	ci

Tx	+	di,	i	ൌ	1,	...,	m

and	the	equivalent	semidefinite	program	is	of	the	form	

minx cTx

subject	to		

	 	

Geometrically,	linear,	conic,	and	semidefinite	programs	are	of	the	form	

minx cTx

subject	to		

												 Ax	+	b K

ሺci

Tx	൅	diሻI Aix	൅	bi

≻ 0,				i ൌ	1,	..., m

ሺAix	൅	biሻT		 ሺci

Tx	൅	diሻI
Page

23

where	K	is	a	pointed1	convex	cone	with	a	nonempty	interior.	Typically	these	

methods	are	applied	to	avoid	local	minima	by	relaxing	the	objective	function	and	

domain	to	be	convex.

	

5.		Decomposition	Techniques	

	

	 Decomposition	methods	are	applied	for	two	purposes,	to	fit	the	decomposed	

problem	into	high‐speed	memory	and	to	decompose	the	problem	into	easier	to	

solve	subproblems.	Decomposition	methods	usually	employ	either	column	

generation	or	cutting	planes,	where	one	approach	is	the	dual	of	the	other.	

	 Bender’s	decomposition	is	a	cutting	plane	method	where	the	primal	

nonlinear	problem	is	

f*	ൌ	minx,y		fሺxሻ ൅	cTy		

subject	to	 	 	 	 	 dual	variable	

gሺxሻ ൅	Ay	൑	b		 	 			 λ		

and	its	dual		

L*	ൌ	maxλminx,y		fሺxሻ	൅	cTy	൅	λሺgሺxሻ	൅	Ay	െ	bሻ.	 	

The	problem	can	be	partitioned	where,	given	x,	we	can	solve	for	variable	y	using	the	

following	inner	minimization	of	the	primal	

miny	cTy			

subject	to		

Ay	൑	b	െ	gሺxሻ		

and	maximization	of	the	dual

maxλ	ሺb‐	gሺxሻሻTλ	

subject	to	

λA	൒	cT.		

The	above	inner	minimization	is	a	linear	program	such	that		

	 λk	ൌ	argmaxλ	ሼሺb	‐	gሺxkሻሻTλ	|	λA	൑	b,	λ	൒	0ሽ,

f*	ൌ	minimize	x		fሺxሻ	൅	y0,	and

y0		൒	ሺb	‐	gሺxሻሻTλk		k	ൌ	0,	1,	...	

for	k	ൌ	0,	1,	...	K.	

Bender’s	Algorithm:		

Step	1.	Set	k	ൌ0;	choose	a	feasible	solution	x0.	

																																																								 1

 K is pointed if it does not contain any subspace except the origin. 
Page

24

Step	2.	Solve	λk൅1		ൌ argminλ	ሼሺb‐	gሺxkሻሻTλ	|	λA	൒	cTሽ.	

Step	3.	Solve	xk൅1	ൌ	argminxሼfሺxkሻ	൅	y0	|	 y0൒	ሺb	‐	gሺxkሻሻTλkሽ	where	the	inner	

minimization	solves	for	y0.	

Step	4.	If	xk൅1	satisfies	the	convergence	criteria,	stop;	otherwise,	set	k	ൌ	k	൅	1	and	go	

to	Step	2.	

Dantzig‐Wolfe	solves	a	simpler	or	smaller	subproblem	that	generates	a	column	for	

the	master	problem:			

min	x,y		fሺxሻ	൅	cTy	

subject	to	 	 	 	 	 dual	variable	

gሺxሻ	൅	Ay	൑	b	.	 	 			 λ

y	є	Y

Therefore	the	Dantzig‐Wolfe	algorithm	is	as	follows:	

Step	1.	Set	k	ൌ	0;	choose	a	feasible	solution	ሺx0,	y0ሻ	where	gሺx0ሻ	൅	Ay0	൑	b.	

Step	2.	Generate	a	column	and	solve	the	modified	master	problem	for	y	and	uk:	

minuk,y		∑k	fሺxkሻuk		൅	cTy	 	

subject	to	 	 	 	 	 dual	variable	

∑k	gሺxkሻuk	൅	Ay		൑	λk

∑k	uk	ൌ	1		 	 	 	 	 	 βk 		

uk		൒	0	

y	є	Y	

Step	3.	Solve	the	subproblem	xk൅1	ൌ	argminx ሼfሺxሻ ൅	λkgሺxሻሽ.		

Step	4.	If	xk		ൌ	∑k	xkuk	satisfies	the	convergence	criteria,	stop;	else	k	ൌ	k	൅	1	and	go	

to	Step	2.		

	 Equivalent	to	both	steps	2	and	3,	the	dual	can	be	solved	by	generating	a	

cutting	plane:	

	 max	λk,y		λkb	൅	βk		

subject	to	

	 λkA	ൌ	cT		

	 λkgሺxkሻ	൅	βk		൒		fሺxkሻ		

	 λk൒	0.	
Page

25

	 Then	in	Step	4,	If	fሺxሻ and	gሺxሻ are	convex, ∑k fሺxkሻuk ൒	∑k	fሺxkሻ	and	∑k

gሺxkሻuk	൒	∑k	gሺxkሻ.	Furthermore,	if	u0*	,...	uk*	and	yk*		is	a	feasible	solution	then	xk	is	

a	feasible	solution	and	as	k	→	∞,	we	have	that	xk and	yk*	is	an	optimal	to	the	NLP.

Approximations	of	the	ACOPF.	The	ACOPF	can	be	formulated	in	several	ways	ሺsee	

Cain	et	al.	2012ሻ.	The	canonical	ACOPF	in	polar	coordinates	is:		

min	∑n cnሺpnሻ	൅	cnሺqnሻ	

subject	to

pn	െ	pd nൌ	∑m	vnvmሺgnmcosθnm	൅	bnmsinθnmሻ	

qn	െ	qd nൌ	∑m	vnvmሺgnmsinθnm	െ	bnmcosθnmሻ

pn

min	൑	pn	൑	pnmax 		 	

qn

min	൑	qn	൑	qnmax 	 	

vn

min	൑	vn	൑	vnmax

θmin nm	൑	θnm	ൌ	θn	െ	θm	൑	θmax nm	 	

Please	refer	to	Cain	et	al.	ሺ2012ሻ	for	notation	definitions	and	the	equivalent	ACOPF	

formulation	in	rectangular	coordinates.	For	the	real	and	reactive	power	flow	

equations,	the	first	partial	derivatives	on	the	power	injections	pn	and	qn	are:	

				∂pn/∂vmൌ	vnሺgnmcosθnm	൅	bnmsinθnmሻ	 		n	്	m

			∂pn/∂θmൌ	vnvmሺgnmsinθnm െ	bnmcosθnmሻ			n	്	m

			∂pn/∂vnൌ	∑mሺ്nሻ	vmሺgnmcosθnm	൅	bnmsinθnmሻ	൅	2vngnn

			∂pn/∂θnൌ	∑mሺ്nሻ	vnvmሺെgnmsinθnm ൅	bnmcosθnmሻ	െ	vn

2bnn

			∂qn/∂vm	ൌ	vnሺgnmsinθnm	െ	bnmcosθnmሻ			n	്	m

			∂qn/∂θm		ൌ	െvnvmሺgnmcosθnm	൅	bnmsinθnmሻ			n	്	m

			∂qn/∂vn	ൌ	∑mሺ്nሻ vmሺgnmsinθnm	െ	bnmcosθnmሻ	െ	2vnbnn

		∂qn/∂θn		ൌ	∑mሺ്nሻ	vnvmሺgnmcosθnm	൅	bnmsinθnmሻ	െ	vn

2gnn

The	first	partial	derivatives	are	used	in	approximation	methods.	Note	that	in	these	

derivatives	the	magnitude	of	the	cosine	and	sine	terms	are	generally	scaled	by	the	

multiple	of	the	conductance	or	susceptance.	Since	reactance	of	a	transmission	line	is	

much	larger	than	resistance,	the	conductance	ሺreal	part	of	admittance	matrix,	Bሻ	is	

much	smaller	than	the	susceptance	ሺimaginary	part	of	admittance	matrix,	Gሻ.	In	

most	cases	the	sine	terms	tend	to	be	small	whereas	the	cosine	term	is	near	unity.	

Upon	inspection,	the	result	of	these	physical	properties	results	in	real	power	that	is	

highly	sensitive	to	changes	in	voltage	angle	and	reactive	power	that	is	highly	

sensitive	to	changes	in	voltage	magnitude.	
Page

26

Decoupled	Power	Flow	Model.	Initially	proposed	by	Stott	and	Alsac	in	1974,	the	

Decoupled	Power	Flow	ሺDPFሻ	is	based	on	the	principle	that	when	the	Jacobian	of	the	

power	flow	equations	

∂p/∂θ ∂p/∂v

∂q/∂θ ∂q/∂v

is	evaluated	numerically,	the	off‐diagonal	submatrices	are	much	smaller	in	

magnitude	than	the	diagonal	submatrices:	

	 ∂p/∂θ		൐൐	∂q/∂θ		and	

	 ∂q/∂v		൐൐	∂p/∂v.	

We	therefore	set	the	off‐diagonal	entries	of	the	Jacobian	to	zero	

∂q/∂θ	ൎ	0		and

	 ∂p/∂v	ൎ	0,	

and	then	decompose	the	problem	into	a	pair	of	subproblems	where	the	p‐θ	real	

power	model	minimizes	the	system	costs	and	the	q‐v	reactive	power	model	

minimizes	the	real	power	transmission	losses.	In	this	approach	we	assume	that	θnm

ൎ	0,	sinθnm	ൎ	θnm,	cosθnm	ൎ	1, and	gnm	൏൏	bnm.	Let	

b’nm	ൌ	ሾ∂pn	/∂θm	ሿ	ൌ	െvnvmbnm					n	്	m,

b’nn	ൌ	0,

	 b”nm	ൌ	ሾ∂qn	/∂vm	ሿ	ൌ	െ	vnbnm 					n	്	m,	and		

b”nn	ൌ	െ	2vnbnn	

where	b’	is	an	approximation	of	the	matrix	of	partial	derivatives	of	the	real	power	

flow	equations	with	respect	to	the	bus	voltage	phase	angles	and b’’	is	an	

approximation	of	the	matrix	of	partial	derivatives	of	the	reactive	power	flow	

equations	with	respect	to	the	bus	voltage	magnitudes.		

The	Decoupled	Power	Flow	model	is:	

	 min	∑n cnሺpnሻ	൅	cnሺqnሻ

	 subject	to	 	 	 	 	 	

	 pn	‐	pd n	ൌ∑m	pnm	൅	sp n		 	

	 qn	‐	qd n	ൌ	∑m	qnm൅	sq n	 	 	 	 		

	 pmin nm	൑	pnm	ൌ	b’nmθmn	൑	pmax nm	

	 qmin nm	൑	qnm	ൌ	b”nmvm	൑	qmax nm	

	 pmin n	൑	pn	൑	pmax n		 	

	 qmin n	൑	qn	൑	qmax n		 	

	 vmin n	൑	vn	൑	vmax n

	 θmin nm	൑	θnm	ൌ	θn	‐	θm	൑	θmax nm	
Page

27

where	

	 pn	ൌ	∑m	vnvmbnmθnm		

	 qn	ൌ	∑m	‐vnvmbnm.

The	Decoupled	Power	Flow	algorithm:	

Step	1.	Set	k	ൌ	0,	and	initialize	vnk	ൌ	1.		

Step	2.	Solve	DPF	for	pnk,	θnk	,	qnk,	and	vnk.	

Step	3.	If	pnk,	θnk	,	qnk,	and	vnk satisfy	the	convergence	criteria,	e.g.,	ac	feasibility,	then	

stop.	If	not,	set	b’nm ൌ	vnkvmkbnm	and	b”nm ൌ	vnkbnm	and	k	ൌ	k	൅1,	and go	to	Step	2.	

Since	the	underlying	problem	is	nonconvex,	there	is	no	guarantee	of	convergence	in		

the	original	ACOPF	problem.	

Bθ	Model.	A	further	simplification	drops	reactive	power	completely.	If	gij	ൌ	0,	sinθ	ൌ	

θ,	cosθ	ൌ	1,	and	vm	ൌ	1,	the	Btheta	models	solves:	

min	∑n cnሺpnሻ		

subject	to	 	 	 	 	 	

pn െ	pd n	ൌ	∑mpnm												 	 	 	

pmin nm	൑	pnm	ൌ	bnmθmn	൑	pmax nm	

pmin n	൑	pn	൑	pmax n		 	

θmin nm	൑	θnm	ൌ	θn െ	θm	൑	θmax nm

Distribution	Factor	Model.	A	further	simplification	is	the	distribution	factor	ሺDFሻ	

model	where	all	transactions	are	decomposed	into	a	‘sale’	to	a	reference	node	and	

‘purchase’	from	a	reference	node.		

	 min	∑n	cnpn	 	 	 	

subject	to	 	 	 	 	 	

	 ∑n	pn െ		∑npd n ൌ	0	 	 	 	 	 	

	 pmin n	൑	pn	൑	pmax n		 	 	 	

	 	pnmk	ൌ	∑n	dfknሺpnെpd nሻ		൑	pmax nmk		 	 	

Since	we	have	assumed	losses	are	zero,	superpositioning	makes	DF	roughly	

equivalent	to	the	Bθ formulation	if	θmax nm	is	equivalent	to pmax nmk.	The	difference	is	

the	size	of	the	problem;	whereas	Bθ	has	N	൅	2K	constraints,	the	distribution	factor	

model	has	1	൅	KM	constraints	where	M	is	the	subset	of	lines	monitored	for	binding	

thermal	limits.	
Page

28

6.	Commercial	Solvers	

	

Commercial	optimization	solvers	for	nonlinear	problems	vary	in	algorithmic	

techniques	and	implementation.	We	survey	five	well‐known	nonlinear	solvers:	

MINOS,	IPOPT,	SNOPT,	KNITRO,	and	CONOPT.	Please	refer	to	ሺCastillo,	2013ሻ	for	a	

computational	study	on	the	performance	of	these	solvers	applied	to	the	ACOPF.	

MINOS.	MINOS	is	a	GRG	method	and	is	designed	to	solve	large‐scale	optimization	

problems	ሺMurtagh	and	Saunders,	1982	and	2003ሻ.	MINOS	partitions	variables	into	

linear	and	nonlinear	elements	and	then	iteratively	solves	the	subproblems	with	

linearized	constraints	and	an	augmented	Lagrangian	objective	function	ሺRobinson,	

1972ሻ.	Instead	of	a	nonlinear	conjugate‐gradient	approach,	a	quasi‐Newton	

approach	is	applied	to	certain	subspaces.	Sparse	LU	basis	factors	are	maintained	by	

LUSOL	ሺGill	et	at,	1987ሻ.	The	nonlinear	constraints	may	be	satisfied	only	in	the	limit,	

so	that	feasibility	and	optimality	may	occur	simultaneously.	An	important	feature	is	

a	stable	implementation	of	a	quasi‐Newton	algorithm	for	optimizing	the	superbasic	

variables.		

IPOPT	ሺInterior	Point	OPTimizerሻ.	IPOPT	converts	the	problem	to	a	barrier	problem.	

Ideally	the	functions	fሺxሻ	and	gሺxሻ	are	twice	continuously	differentiable.	IPOPT	uses	

line	filtered	searches	and	includes	a	feasibility	restoration	phase.	Filter	methods	

promote	global	convergence	through	measuring	decreases	in	the	objective	function	

and	infeasibility	as	two	separate	criteria	that	are	controlled	simultaneously.	IPOPT	

has	options	in	line	search	strategies	for	globalization,	including	an	exact	penalty	

merit	function,	augmented	Lagrangian	merit	function,	filter	method	ሺR.	Fletcher,	S.	

Leyffer,	and	P.	Tointሻ,	Hessian	and	several	Hessian	approximation	methods.	

SNOPT	ሺSparse	Nonlinear	OPTimizerሻ.	SNOPT	implements	a	sparse	active‐set	

sequential	quadratic	programming	ሺSQPሻ	method	that	employs	quasi‐Newton	

approximations	to	determine	the	Hessian	in	the	quadratic	programming	

subproblem;	then	an	augmented	Lagrangian	merit	function	guides	the	line	search	

direction.	Sparse	basis	factors	are	maintained	by	LUSOL.	If	only	the	objective	is	

nonlinear,	the	problem	is	linearly	constrained	ሺLCሻ	and	tends	to	solve	more	easily	

than	the	general	case	with	nonlinear	constraints	ሺNCሻ.	SNOPT	uses	limited‐memory	

quasi‐Newton	approximations	to	the	Hessian	of	the	Lagrangian.	The	merit	function	

for	step‐length	control	is	an	augmented	Lagrangian,	as	in	the	dense	SQP	solver	

NPSOL	ሺGill	et	at,	1992ሻ.	In	general,	SNOPT	requires	less	matrix	computation	than	

NPSOL	and	fewer	evaluations	of	the	functions	than	the	nonlinear	algorithms	in	

MINOS.	It	is	most	efficient	if	only	some	of	the	variables	enter	nonlinearly,	or	there	

are	relatively	few	degrees	of	freedom	at	a	solution	ሺi.e.,	many	constraints	are	activeሻ.	

	
Page

29

KNITRO.	KNITRO	implements	both	interior‐point	and	active‐set	methods	for	solving	

nonlinear	optimization	problems.	The	variables	can	be	continuous,	binary,	or	

integer.	In	the	barrier/interior	method,	KNITRO	solves	a	series	of	barrier	sub‐

problems	controlled	by	a	barrier	parameter.	The	algorithm	uses	trust	regions	and	a	

merit	function	to	promote	convergence.	The	algorithm	performs	one	or	more	

minimization	steps	on	each	barrier	problem,	then	decreases	the	barrier	parameter,	

and	repeats	the	process	until	the	original	problem	has	been	solved	to	the	desired	

accuracy.	KNITRO	provides	two	procedures	for	computing	the	steps.	One	version	

where	each	step	is	computed	using	a	projected	conjugate	gradient	iteration.	This	

approach	factors	a	projection	matrix	to	approximately	minimize	a	quadratic	model	

of	the	barrier	problem.	The	other	procedure	always	attempts	to	compute	a	new	

iterate	by	solving	the	primal‐dual	KKT	matrix	using	direct	linear	algebra,	but	if	the	

step	quality	cannot	be	guaranteed	or	if	negative	curvature	is	detected,	then	the	new	

iterate	is	computed	by	the	first	procedure.	KNITRO	also	implements	an	active‐set	

sequential	linear‐quadratic	programming	ሺSLQPሻ	algorithm	that	uses	linear	

programming	sub‐problems	to	estimate	the	active‐set	at	each	iteration.	KNITRO	

includes	an	active‐set	SLQP	with	trust	region	to	promote	convergence.	

	

CONOPT	ሺCONstrained	OPTimizationሻ.	CONOPT	is	a	generalized	reduced	gradient	

ሺGRGሻ	method	that	searches	along	the	steepest	descent	direction	in	the	superbasic	

variables.	The	generalized	projection	method	projects	the	search	direction	into	the	

subspace	tangent	to	the	active‐set	of	constraints.	The	active‐set	is	the	subset	of	

equality	and	inequality	constraints	at	a	point	that	satisfy	with	equality;	for	example,	

equality	constraints	are	active	at	all	feasible	points.	Therefore,	the	search	direction	

is	projected	into	the	null	space	of	the	gradients	for	the	equality	and	binding	

inequality	constraints.	The	projected	gradient	could	be	infeasible,	which	then	

requires	a	correction	step.		

7.	History	of	ACOPF	Solution	Techniques	

	

	 In	1962,	Carpentier	introduced	the	alternating	current	optimal	power	flow	

ሺACOPFሻ	for	economic	dispatch	based	upon	Karush‐Kuhn‐Tucker	ሺKKTሻ	conditions.		

Carpentier	employed	the	Gauss‐Seidel	method	representing	the	load	flows	as	power	

injections	in	the	voltage	polar	form.	Carpentier	included	operational	constraints	on	

real	power	control,	generator	bus	voltage	magnitude	limits,	reactive	power	control	

of	switchable	VAR	sources,	and	transformer	tap	setting.	This	formulation	has	

nonconvex	equality	constraints	with	quadratic	and	trigonmetric	functions	ሺsee	Cain	

et	al,	2012ሻ.	The	complex	variables	can	be	expressed	in	polar	or	rectangular	

coordinates	and	results	in	different	types	of	nonconvex	constraints.		Since	the	

ACOPF	problem	is	a	nonconvex	mathematical	program,	the	KKT	conditions	yield	
Page

30

only	a	local	optimal	solution;	even	feasibility	cannot	be	guaranteed	by	most	

nonlinear	programming	ሺNLPሻ	algorithms.		

	 In	1969,	Carpentier	and	Abadie	published	the	generalized	reduced	gradient	

method,	which	is	a	generalization	of	Carpentier’s	‘differential	injections’	method	

that	was	originally	conceived	in	1964	to	solve	the	optimal	power	flow	problem	

ሺCarpentier,	1979ሻ.	Since	Carpentier’s	contributions,	there	has	been	a	wealth	of	

research	done	on	algorithmic	methods	to	solve	the	ACOPF.	In	1968,	Dommel	and	

Tinney	presented	an	approach	based	on	power	flow	solution	by	Newton’s	method	

with	a	gradient	adjustment	factor	for	the	penalty	function	to	account	for	dependent	

constraints.	In	1969,	Sasson	used	the	Fletcher‐Powell	method;	then	in	1973,	Sasson	

directly	computed	the	Hessian	and	utilized	its	sparsity.		In	1973,	Alsac	and	Stott	

incorporated	exact	outage‐contingency	constraints	and	used	an	augmented	

Lagrangian.		In	1977,	Biggs	and	Laughton	used	a	recursive	equality	quadratic	

program.		In	the	above	studies,	test	problems	with	less	than	30‐buses	were	used.	

During	this	time	the	advancements	of	gradient	techniques	with	Newton’s	method	

were	widely	adopted	in	the	early	ACOPF	algorithms,	see	Happ	ሺ1977ሻ.		In	1979,	Wu	

et	al.	used	a	modified	reduced	gradient	with	penalty	function	methods	and	included	

the	largest	test	case	published	to‐date,	a	1410‐bus	system.		

	 In	the	following	table	we	summarize	the	published	optimal	power	flow	

studies	by	the	test	problems,	convergence	criteria,	initialization,	and	performance	

metrics	reported.	The	studies	vary	in	all	the	metrics	listed	in	Table	1,	and	very	few	

of	the	studies	included	a	comparison	against	commercial	solvers	or	previously	

published	data.	For	example	out	of	the	studies	listed,	Burchett	et	al.	ሺ1982ሻ,	Aoki	

and	Kanezashi	ሺ1985ሻ,	Huneault	and	Galiana	ሺ1990ሻ,	Ponrajah	and	Galiana	ሺ1990ሻ,	

Momoh	and	Zhu	ሺ1999ሻ,	and	Wang	et	al.	employed	MINOS	in	their	algorithmic	

implementation	or	compared	computational	results	to	that	of	MINOS.		

	 The	abbreviations	in	Table	1	are	as	follows:	N/A	ሺNot	Applicableሻ,	NR	ሺNot	

Reportedሻ,	FS	ሺFlat	Startሻ,	LFS	ሺLoad	Flow	Startሻ,	RS	ሺRandom	Startሻ,	HS	ሺHot	Startሻ,	

WSሺWarm	Startሻ,	US	ሺUser	Specified	Startሻ,	C.G.	ሺComplementarity	Gapሻ,	Opt.	

ሺOptimality	Toleranceሻ,	and	Feas.	ሺFeasibility	Toleranceሻ.	Note	for	the	SDP	methods	

in	Bai	ሺ2008ሻ	and	Lavaei	ሺ2012ሻ,	initialization	is	not	required.	

	 The	studies	with	US	ሺUser	Specifiedሻ	listed	as	the	initialization	method	

employed	a	non‐random,	selective	approach.	Numerous	studies	reported	additional	

metrics	to	those	listed	in	Table	1.	However,	we	identify	this	handful	of	metrics	as	

basic	information	needed	for	fair	testing	and	comparison	of	proposed	algorithms	by	

other	researchers.	

	 The	following	Table	2	summarizes	the	test	problems	published	on	in	ACOPF	

studies	to‐date.	Numerous	of	the	test	problems	have	a	network	origin	that	is	not	

described	well.	Although	many	of	the	test	problems	are	used	across	multiple	

studies,	as	denoted	in	the	“Source	Data”	column	a	small	percentage	of	the	test	
Page

31

problems	are	readily	available	in	the	public	domain	to	use;	these	test	problems	are	

posted	online	by	the	University	of	Washington	Electrical	Engineering	or	by	

MATPOWER.	The	raw	data	for	test	problems	23	and	39	are	available	in	Biggs	ሺ1977ሻ	

and	Pai	ሺ1989ሻ,	respectively.	We	list	the	referencing	study	if	the	reader	is	interested	

in	further	information	on	the	test	problem	or	to	contact	the	author	for	the	source	

data.		
Page

32

Table	1.	ACOPF	Studies	that	include	Numerical	Analysis	

Year	 Author	 Test	Problems Convergence	 Initial‐ Work‐ 

ሺNo. of	Nodesሻ Criteria	 ization station 

1969	 Sasson	 30	 10^‐5	Opt.	 NR x 1973	 Alsac	 30	 10^‐3	Opt. US 1974	 Mukherjee	 25	 NR US x 1977	 Biggs	 23	 NR US 1979	 Wu	 11,	136,	333,	1410 NR LFS x 1982	 Burchett	 118,	597	 NR NR x 1982	 Divi	 9,	10,	11	 10^‐5	Feas. FS 1982	 Shoults	 5,	30,	962	 NR US x 1984	 Burchett	 350,	1100,	1600,	1900 NR RS x 1984	 Sun	 912	 NR US x 1985	 Aoki	 14,	135	 NR NR x 1988	 Santos	 118,	129	 NR NR x1989	 Nanda	 14,	30,	89	 10^‐4	Feas.	 LFS 1989	 Ponrajah	 6,	10,	30,	118 NR FS x 1990	 Alsac	 1330,	1200,	700 NR HS x 1990	 Huneault	 30,	118	 NR US 1990	 Salgado	 14,	30,	39,	57,	89,	118 10^‐3	Feas. LFS x 1993	 Almeida	 14,	30,	34	 NR FS,LFS 1994	 Momoh	 9,	14,	30,	118 NR US x 1994	 Wu	 9,	30,	39,	118,	244 10^‐6	C.G. FS x 1995	 Chebbo	 706	 50^‐5 Opt.,	50^‐4	Feas.	 US x 1997	 Lai	 30	 NR US x 1998	 Torres	 30,	57,	118,	300 10^‐4	Opt. FS,LFS 1998	 Wei	 14,	30,	57,	118,	344,	703,	1047 10^‐6	Opt.	 US x 
Page

33

Year Author Test	Problems Convergence Initial‐ Work‐ 

ሺNo.	of	Nodesሻ Criteria ization station 

1999	 Momoh	 14	 10^‐6	Opt.	 US 1999	 Yan	 118,	1062	 10^‐8	Opt.,	Feas.		 US x 2000	 Nejdawi	 30,	57,	118,	300 NR FS 2000	 Torres	 30,	57,	118,	300 10^‐4	Opt.	 US x 2001	 Castronuovo	 118,	352,	750,	1704 10^‐6	Opt.,	Feas.		 NR x 2001	 Lima	 810,	2256	 10^‐3	Feas. US 2001	 Torres	 118,	256,	300,	555,	2098 10^‐4	Opt. US x 2001	 Xie	 30,	57,	118	 10^‐5	C.G.,	Feas.		 WS 2002	 Jabr	 14,	24,	30,	57,	118,	175,	300 10^‐4	Opt.	 US x 2002	 Lin	 118,	244	 10^‐3	Opt.	 US 2002	 Torres	 30,	57,	118,	256,	30,	555,	2098 10^‐4	Opt.	 US x 2003	 Oliveira	 30,	118,	1564,	1732,	1993 Sqrt	of	machine	eps.	 US x 2004	 Lin	 118,	244	 10^‐3	Opt.	 US 2005	 Min	 14,	30,	57,	118,	254,	300, 662 10^‐3	Opt.	 LFS x 2005	 Tate	 2,	118,	10274 10^‐3	Feas. FS 2007	 Capitanescu	 60,	118,	300 10^‐6	C.G. RS x 2007	 Lin	 9,	14,	30,	57,	118 10^‐4	C.G.,	Feas.		 FS x 2007	 Sousa	 30,	57,	118,	300 10^‐4	Opt.	 US 2007	 Wang	 30,	57,	118,	300,	2383,	2935 NR FS x 2008	 Bai	 4,	14,	57,	118,	300 10^‐5	Opt. N/A x 2008	 Jabr	 9,	14,	30,	39,	57,	118,	300,	2383 10^‐8	Opt. FS,LFS x 2008	 Lin	 9,	14,	30,	57,	118,	300 10^‐4	C.G.,	Feas.		 FS x 2009	 Bedrinana	 2,	14,	57	 NR NR 2009	 Chiang	 678,	2052,	2383 10^‐6	C.G.,	Feas.		 FS x 
Page

34

Year Author Test	Problems Convergence Initial‐ Work‐ 

ሺNo.	of	Nodesሻ Criteria ization station 

2009	 Jiang	 14,	30,	39,	57,	118,	300,	701,	2052,	

2383,	2736,	2746

10^‐5	C.G.,	Feas.		 FS x 

2009	 Sousa	 30,	118,	300,	2256 10^‐3	Feas.	 NR x 2009	 Yang	 30,	118,	300 NR NR 2010	 Jiang	 5,	6,	9,	14,	30,	39,	57,	118,	2383 NR FS x 2010	 Jiang	et	al.	 118,	300,	678,	2052,	2383,	2746 NR NR x 2010	 Xie	 57,	118,	300,	2052,	2790 10^‐6	C.G. US x 2011	 Chung	 14,	118	 10^‐4	Feas.	 FS x 2011	 Phan	 6,	9,	14,	30,	39,	57,	118,	300,	2746 10^‐5	Opt.	 US x 2011	 Sousa	 30,	57,	118,	300,	1211 10^‐4	Opt.	 FS,LFS,RS 2011	 Zimmerman	 9,	30,	36,	118,	300,	2383,	2736,	3120,	

2935,	21000,	42000

NR NR x 

2012	 Lavaei	 14,	30,	57,	118,	300 NR N/A 
Page

35

Table	2.	Test	Problems	used	in	ACOPF	Studies	

Test	Problem	 Network	Origin Source	Data Referencing	

14	 Midwestern	US	System UW	Electrical	Eng Multiple	stu23	 N/A Biggs,	1977 Biggs,	1924	 IEEE	Reliability	Test System UW	Electrical	Eng Multiple	stu25	 N/A N/A	 Mukherjee,	30	 Midwestern	US	System UW	Electrical	Eng Multiple	stu34	 N/A N/A	 Almeida,	139	 New	England US System Pai,	1989 Phan,	2057	 Midwestern	US System UW	Electrical	Eng Multiple	stu60	 Nordic32	System N/A	 Capitanescu89	 N/A N/A	 Salgado,	1

118	 Midwestern	US System UW	Electrical	Eng Multiple	stu129	 Companhia	Energetica	de Sao	Paulo N/A	 Santos,	19135	 Chugoku	Electric	Power	Co.,	Japan N/A	 Aoki,	19136	 N/A N/A	 Wu,	197175	 N/A N/A	 Jabr,	200244	 N/A N/A	 Wu,	1994	ሺIniti254	 N/A N/A	 Min,	200256	 N/A N/A	 Torres,	2001	a300	 Midwestern	US System UW	Electrical	Eng Multiple	stu333	 N/A N/A	 Wu,	197344	 Japan	System N/A	 Wei,	199350	 Northeastern	US Utility N/A	 Burchett,	1352	 South‐Southeastern	Brazil	System N/A	 Castronuovo555	 N/A N/A	 Torres,	2001	a597	 Interconnection	of	Several	Utilities N/A	 Burchett,	1662	 N/A N/A	 Min,	200678	 N/A N/A	 Chiang,	20
Page

36

Test	Problem Network	Origin Source	Data Referencing	

700	 Interconnection	of	Several	Utilities N/A	 Alsac,	19701	 Undisclosed	Real	Power	System N/A	 Jiang,	20703	 China	System N/A	 Wei,	199706	 N/A N/A	 Chebbo,	1750	 South‐Southeastern	Brazil	System N/A	 Castronuovo810	 South‐Southeastern	Brazil	System N/A	 Lima,	20912	 Northeastern	US System N/A	 Sun,	198962	 16	Interconnected	Areas N/A	 Shoults,	1

1047	 Simulation	Power	System N/A	 Wei,	1991062	 N/A N/A	 Yan,	1991100	 Eastern	United	States	Pool N/A	 Burchett,	11200	 Utility	Company N/A	 Alsac,	191211	 Undisclosed	Real	Power	System N/A	 Sousa,	201330	 Interconnection	of	Several Utilities N/A	 Alsac,	191410	 N/A N/A	 Wu,	1971564	 South‐Southeastern	Brazil	System N/A	 Oliveira,	21600	 Western	United	States	Utility N/A	 Burchett,	11704	 South‐Southeastern	Brazil	System N/A	 Castronuovo1732	 South‐Southeastern	Brazil	System N/A	 Oliveira,	21900	 Northeastern	US Utility N/A	 Burchett,	11993	 South‐Southeastern	Brazil	System N/A	 Oliveira,	22052	 N/A N/A	 Multiple	stu2098	 Modified	Brazil	System N/A	 Torres,	2001	a2256	 South‐Southeastern	Brazil	System N/A	 Lima,	202383	 Polish	System Zimmerman,	2011 Zimmerman2736	 Polish	System Zimmerman,	2011 Zimmerman2746	 Polish	System Zimmerman,	2011 Zimmerman
Page

37

Test	Problem Network	Origin Source	Data Referencing	

2790	 Undisclosed	Real	Power	System N/A	 Xie,	2012935	 Polish	System Zimmerman,	2011 Zimmerman3120	 Polish	System Zimmerman,	2011 Zimmerman

10274	 N/A N/A	 Tate,	2021000	 N/A N/A	 Zimmerman42000	 N/A N/A	 Zimmerman
Page

38

	 Sequential	quadratic	programming	ሺSQPሻ	methods	applied	to	the	OPF	were	

either	linear	programming	ሺLPሻ‐based	QP	methods	or	QP	methods	based	on	the	

concept	of	an	active‐set	of	linearly	independent	constraints.		The	LP	approaches,	

based	on	Wolfe’s	or	Beale’s	algorithm,	solve	a	linear	program	by	a	revised	simplex	

technique.		Exploiting	many	of	the	strengths	of	linear	programming,	these	

approaches	permit	the	use	of	artificial	variables	to	guide	feasibility,	parametric	

programming,	and	ease	of	equality	constraints	into	the	formulation.		Furthermore,	

some	of	the	successive	quadratic	programs	apply	Newton’s	method	to	the	

subproblem.		In	1984	by	comparison,	Burchett	et	al.	present	a	method	where	the	QP	

subproblems	are	equivalent	to	a	sequence	of	Newton	steps	to	the	optimal	solution.		

	 In	1982	Shoults	and	Sun	decomposed	the	problem	into	active	and	reactive	

subproblems.	In	1984,	Sun	applied	Newton’s	method	to	the	SQP	and	advanced	

sparsity	techniques	to	the	decomposition,	but	the	method	exhibited	problems	in	

initialization	and	ill‐conditioning.	In	1986,	Contaxis	et	al.	and	in	1989,	Nanda	et	al.	

applied	Fletcher’s	method	to	the	decoupled	subproblems.		

	 In	1981,	Giras	et	al.	applied	Powell’s	quasi‐Newton	method	to	the	ACOPF,	

which	performs	a	Broyden‐Fletcher‐Goldfarb‐Shanno	ሺBFGSሻ	update	and	shows	

convergence	from	infeasible	starting	points.	In	1982,	Divi	et	al.	applied	Fletcher’s	

quasi‐Newton	method	with	a	shifted	penalty	function,	and	the	Broyden‐Fletcher‐

Shanno	ሺBFSሻ	updating	formula	to	numerically	stabilize	the	Hessian	as	a	positive	

definite	matrix	and	promote	global	convergence	by	means	of	a	penalty	function	that	

ensures	sufficient	progress	along	the	line‐search	towards	the	active‐set	of	

constraints.	In	1982,	Burchett	et	al.	and	in	1985	Aoki	et	al.	solve	a	sequence	of	

sparse,	linearly	constrained	subproblems	based	on	MINOS.	In	1988,	Santos	

proposed	the	dual	augmented	Lagrangian	method	that	attempts	to	address	ill‐

conditioning	of	the	Hessian	by	applying	a	quasi‐Newton	method	to	the	dual	

function.	In	1992,	Monticelli	et	al.,	present	an	improvement	upon	Maria	et	al.	ሺ1987ሻ	

in	which	an	adaptive	penalty	strategy	is	used	to	ensure	positive	definiteness	of	the	

Hessian,	without	negatively	affecting	the	quadratic	convergence	characteristic	of	

Newton’s	method.	

	 These	methods	build	upon	Newton’s	method	for	unconstrained	optimization,	

Lagrange’s	method	for	optimization	with	equalities,	and	Fiacco	and	McCormick’s	

barrier	method.	Meliopoulos	and	Xia	ሺ1993ሻ,	Vargas	et	al.	ሺ1992ሻ,	Momoh	et	al.	

ሺ1994ሻ,	Momoh	and	Zhu	ሺ1999ሻ,	and	Nejdawi	et	al.	ሺ2000ሻ	applied	an	interior	point	

algorithm	ሺIPMሻ	to	a	LP	or	QP,	where	the	constraints	are	linearized. 	Moreover,	

primal‐dual	IPMs	have	been	successful	in	solving	the	ACOPF	by	introducing	a	

logarithmic	barrier	function	in	place	of	the	inequality	constraints.	Wu	et	al.	ሺ1994ሻ,	

Wei	et	al.	ሺ1998ሻ,	Torres	and	Quintana	ሺ1998ሻ,	Yan	and	Quintana	ሺ1999ሻ,	

Castronuovo	et	al.	ሺ2001ሻ,	Xie	and	Song	ሺ2001ሻ,	Wang	et	al.	ሺ2007ሻ,	Xie	and	Chiang	
Page

39

ሺ2010	and	2011ሻ,	Sousa	et	al.	ሺ2011ሻ,	and	Chung	et	al.	ሺ2011ሻ	present	a	primal‐dual	

interior	point	method	to	solving	the	ACOPF.			

	 In	1994,	Granville	solves	the	primal‐dual	without	Mehrotra’s	predictor‐

corrector	mechanism,	and	notes	that	the	contribution	of	barrier	terms	into	the	

diagonal	of	the	Hessian	matrix	are	very	effective	in	bringing	positive	definiteness	to	

the	problem,	therefore	making	auxiliary	penalty	functions	unnecessary.		In	1998,	

Torres	et	al.	exploit	the	rectangular	formulation	of	the	power	flow	constraints,	

which	is	quadratic,	using	an	interior	point	method.		This	results	in	a	constant	

Hessian	and	a	Taylor	expansion	terminating	at	the	second‐order	term	that	reduces	

the	computational	burden	and	iterations	to	convergence.		Furthermore,	Torres	et	al.	

ሺ1998ሻ	perturb	the	boundary	to	deal	with	numerical	ill‐conditioning	that	may	occur	

with	binding	constraints.		In	1998,	Wei	et	al.	use	an	interior‐point	method	based	on	

applying	Newton’s	method	to	the	nonlinear	system	of	perturbed	KKT	conditions,	

which	is	similar	to	the	approach	by	Torres	and	Quintana	ሺ2001ሻ;	this	method	

promotes	global	convergence.	In	1998,	Wei	et	al.	present	a	data	structure	that	

further	reduces	the	nodal	block	fill‐in	elements.		In	comparison	to	storing	the	

augmented	system	in	compact	blocks,	as	noted	by	Sun	et	al.	ሺ1984ሻ	and	Wei	et	al.	

ሺ1998ሻ,	Castronuovo	et	al.	ሺ2001ሻ	propose	a	vectorization	technique	that	only	

considers	nonzero	terms	in	order	to	decrease	computational	cost	per	operation.	In	

1999,	Yan	and	Quintana	ሺ1999ሻ	present	a	dynamic	adjustment	of	the	step	size	and	

tolerance	to	improve	convergence	speed.	

	 Wu	et	al.	ሺ1994ሻ,	Torres	and	Quintana	ሺ1998ሻ,	Yan	and	Quintana	ሺ1999ሻ,	

Castronuovo	et	al.	ሺ2001ሻ,	Lin	et	al.	ሺ2007ሻ,	Wang	et	al.	ሺ2007ሻ,	Lin	et	al.	ሺ2008ሻ,	Xie	

and	Chiang	ሺ2010	and	2011ሻ,	Sousa	et	al.	ሺ2011ሻ	and	Chung	et	al.	ሺ2011ሻ	use	

Mehrotra’s	predictor‐corrector.		In	2009,	Sousa	et	al.	present	a	predictor‐corrector	

modified	barrier	approach,	based	on	Polyak’s	modified	barrier	method,	to	address	

ill‐conditioning.	However,	the	corrector	step	can	lead	to	very	slow	convergence	or	

failure.	In	2001,	Torres	and	Quintana	apply	Gondzio’s	multiple	centrality	corrections	

strategy.	Furthermore,	Xie	and	Song	ሺ2001ሻ	and	Chung	et	al.	ሺ2011ሻ	leverage	a	

nodal	block	data	structure	within	the	interior	point	method	iterates	in	order	to	

reduce	the	correction	equation.	

	 	Due	to	the	nonlinearities	and	ill‐conditioning	of	the	ACOPF,	recent	research	

has	focused	on	applying	methods	with	particularly	robust	global	convergence	

properties.	Jabr	et	al.	ሺ2002ሻ	present	a	primal‐dual	interior	point	method	that	

replaces	the	Hessian	with	a	“2‐norm	positive	approximant”	with	a	“watch‐dog”	

strategy.	Chiang	et	al.	ሺ2009ሻ	present	a	two‐stage	solution	algorithm	where	an	

active‐set	quotient	gradient	method	is	applied	in	stage	one	to	induce	global	

convergence,	and	an	interior‐point	method	is	applied	in	stage	two	to	obtain	a	local	

optimal	solution.	Min	and	Shengsong	ሺ2005ሻ,	Pajic	and	Clements	ሺ2005ሻ,	Zhou	et	al.	

ሺ2005ሻ,	Sousa	and	Torres	ሺ2007ሻ,	Wang	et	al.	ሺ2007ሻ,	Chiang	et	al.	ሺ2009ሻ,	and	
Page

40

Sousa	et	al.	ሺ2011ሻ	apply	a	trust	region	method.	Similar	to	Torres	and	Quintana	

ሺ2001ሻ,	Min	and	Shengsong	ሺ2005ሻ	also	apply	Gondzio’s	multiple	centrality	

corrections	strategy	to	the	primal‐dual	IPM,	but	instead	incorporate	the	trust‐region	

method	in	solving	the	subproblems.		The	infinity	norm,	which	forms	a	closed,	

compact,	convex	hypercube	in	n‐dimensional	space,	is	used	to	define	the	trust	

region,	and	a	merit	function	is	used	to	determine	whether	to	accept	or	reject	a	trial	

step	to	the	trust	region	subproblem.		The	authors	also	introduce	a	feasibility	

restoration	variable,	which	controls	the	feasibility	requirements	of	the	trust	region	

subproblem.			

	 Also	known	as	parametric	continuation	methods,	homotopy	methods	are	

path‐following	approaches	which	explicitly	progress	towards	a	solution	to	the	

original	nonlinear	problem.	In	doing	so	these	methods	construct	and	solve	a	new	

simpler	problem	compared	to	the	original	one	and	then	gradually	reconstruct	the	

original	system	of	equations	in	order	to	solve	for	the	unknown	solution.	Ponrajah	

and	Galiana	ሺ1989ሻ,	Huneault	and	Galiana	ሺ1990ሻ,	Almeida	et	al.	ሺ1993ሻ,	Lima	et	al.	

ሺ2001ሻ,	and	Jiang	et	al.	ሺ2010ሻ	apply	homotopy	methods	to	solve	the	ACOPF.	

	 Most	recently,	convex	optimization	techniques	have	been	applied	to	the	

ACOPF.	In	2007,	Jabr	reformulated	the	ACOPF	as	a	second‐order	conic	program	and	

applied	an	interior	point	method	in	MOSEK	that	is	specific	to	conic	quadratic	

optimization	that	has	polynomial	convergence.	In	2012,	Lavaei	and	Low	apply	

semidefinite	programming	optimization	to	solve	the	ACOPF,	and	prove	global	

optimality	under	a	sufficient	zero‐duality	gap	condition.	

	 Another	class	of	optimization,	derivative	free	optimization,	is	typically	

applied	when	first‐	and	second‐derivatives	are	not	available	or	are	expensive	to	

compute.	Lai	et	al.	ሺ1997ሻ,	Iba	ሺ1994ሻ,	Abido	ሺ2002ሻ	and	many	others	apply	

derivative	free	approaches	to	solve	the	OPF.	

8.	Summary	and	Conclusions	

	

	 For	the	last	fifty	years,	the	latest	developments	in	nonlinear	optimization	

have	been	applied	to	the	ACOPF	in	hopes	of	better	solution	techniques	for	large‐

scale,	practical	network	operations	and	planning.	Until	recently,	the	dominant	

formulation	has	been	in	polar	coordinates.	Although	the	research	to‐date	presents	a	

relatively	positive	picture,	clearly	the	published	experimental	results	have	been	

limited.	The	lack	of	reported	metrics	and	available	test	problems	makes	it	difficult	to	

perform	a	comparative	assessment	of	proposed	solution	techniques	to‐date.	

Furthermore,	there	is	a	significant	lack	of	independent	studies	that	compare	the	

numerous	approaches	in	a	systematic	manner;	for	example,	refer	to	the	

benchmarking	study	completed	by	Mittelman	ሺ2012ሻ	where	he	reports	on	solver	

performance	for	optimizing	general	nonlinear	problems.	In	a	companion	study	
Page

41

Castillo	and	O’Neill	ሺ2013ሻ	report	numerical	results	from	testing	CONOPT,	IPOPT,	

KNITRO,	MINOS,	and	SNOPT	on	various	sized	test	problems	in	which	we	apply	

various	mathematically	equivalent	AC‐OPF	formulations	with	numerous	starting	

points.		

	 Although	an	approximation	of	the	ACOPF	is	solved	in	practice,	current	

market	operations	and	planning	do	not	co‐optimize	real	power	dispatch	with	

voltage	and	reactive	power	management.	A	recent	white	paper	by	Stott	and	Alsac	

ሺ2012ሻ	argue	that	the	ACOPF	is	still	very	much	a	‘work	in	progress.’	In	fact,	the	

authors	contend	that	most	of	the	work	completed	on	ACOPF	solution	techniques	to‐

date	is	not	adaptable	to	real	world	applications.	The	ACOPF	would	often	be	

embedded	within	bigger	calculations	that	make	convergence	and	computational	

effort	even	more	difficult.	Therefore	in	developing	and	testing	future	solution	

techniques,	it	is	important	to	understand	how	market	operations	and	planning	

would	leverage	the	ACOPF	in	real	world	applications.	

References	

J.	Abadie,	J.	Carpentier.		“Generalization	of	the	Wolfe	Reduced	Gradient	Method	to	

the	Case	of	Nonlinear	Constraints”	in	Optimization	by	R.	Fletcher,	Pp.	37‐49,	

Academic	Press,	London,	1969.	

M.A.	Abido.	“Optimal	Power	Flow	Using	Particle	Swarm	Optimization.”	International	

Journal	of	Electrical	Power	&	Energy	Systems,	Vol.	24,	No.	7,	Pp.	563‐571,	2002.	

K.C.	Almeida,	F.D.	Galiana,	S.	Soares.	“A	General	Parametric	Optimal	Power	Flow.”	

Power	Industry	Computer	Application	Conference	Proceedings,	1993.	

O.	Alsac,	B.	Stott.	“Fast	Decoupled	Load	Flow.”	IEEE	Transactions	on	Power	

Apparatus	and	Systems,	PAS‐93,	Pp.	859‐869,	May	1974.	

O.	Alsac,	B.	Stott.	“Optimal	Load	Flow	with	Steady‐State	Security.”	IEEE	PES	Summer	

Meeting,	July	15‐20,	1973.		

O.	Alsac,	J.	Bright,	M.	Prais,	B.	Stott.	“Further	Developments	in	LP‐Based	Optimal	

Power	Flow.”		IEEE	Transactions	on	Power	Systems,	Vol.	5,	No.	3,	August	1990.	

K.	Aoki,	M.	Kanezashi.		“A	Modified	Newton	Method	for	Optimal	Power	Flow	Using	

Quadratic	Approximated	Power	Flow.”		IEEE	Transactions	on	Power	Apparatus	and	

Systems,	Vol.	PAS‐104,	No.	8,	August	1985.	

X.	Bai,	H.	Wei,	K.	Fujisawa,	Y.	Wang.	“Semidefinite	Programming	for	Optimal	Power	

Flow	Problems.”	Electrical	Power	and	Energy	Systems,	Vol.	30,	Pp.	383‐392,	2008.	

R.H.	Bartels.	“A	Stabilization	of	the	Simplex	Method.”	Numerische	Mathematik,	Vol.	

16,	Pp.	441‐434,	1971.	

R.H.	Bartels	and	G.H.	Golub.	“The	Simplex	Method	of	Linear	Programming	Using	the	

LU	Decomposition.”	Communications	of	the	ACM,	Vol.	12,	Pp.	266‐268.	
Page

42

E.M.L.	Beale.		“An	Introduction	to	Beale’s	Method	of	Quadratic	Programming”	in	

Nonlinear	Programming,	J.	Abadie,	Amsterdam:	North‐Holland,	1967.	

M.	F.	Bedrinana,	M.J.	Rider,	C.A.	Castro.	“Ill‐Conditioned	Optimal	Power	Flow	

Solutions	and	Performance	of	Non‐Linear	Programming	Solvers.”	IEEE	Bucharest	

Power	Tech	Conference,	June	28th‐July	2nd,	2009.	

L.	T.	Biegler,	“Nonlinear	Programming:	Concepts	and	Algorithms	for	Process	

Optimization.”	Chemical	Engineering	Department,	Carnegie	Mellon	University	

M.C.	Biggs,	M.A.	Laughton.		“Optimal	Electric	Power	Scheduling:	A	Large	Nonlinear	

Programming	Test	Problem	Solved	by	Recursive	Quadratic	Programming.”		

Mathematical	Programming,	Vol.	13,	Pp.	167‐182,	1977.	

R.C.	Burchett,	H.H.	Happ,	K.A.	Wirgau.		“Large	Scale	Optimal	Power	Flow.”		IEEE	

Transactions	on	Power	Apparatus	and	Systems,	Vol.	PAS‐101,	No.	10,	October	1982.	

R.C.	Burchett,	H.H.	Happ,	D.R.	Vierath.		“Quadratically	Convergent	Optimal	Power	

Flow.”		IEEE	Transaction	on	Power	Apparatus	and	Systems,	Vol.	PAS‐103,	No.	11,	

November	1984.	

R.	Byrd,	J.C.	Gilbert,	J.	Nocedal.	"A	Trust	Region	Method	Based	on	Interior	Point	

Techniques	for	Nonlinear	Programming."	Math	Prog	A,	Vol.	89,	Pp.	149‐185,	2000.	

R.	H	Byrd,	R.	B.	Schnabel,	and	G.	A.	Schultz.	"A	Trust	Region	Algorithm	for	

Nonlinearly	Constrained	Optimization",	SIAM	J.	Numer.	Anal.,	Vol.	24,	Pp.	1152–

1170,	1987.	

M.B.	Cain,	R.P.	O’Neill,	A.	Castillo.	“Optimal	Power	Flow	Papers:	Paper	1.	History	of	

Optimal	Power	Flow	and	Formulations.”	Federal	Energy	Regulatory	Commission	

ሺFERCሻ	website,	2012.	

J.	Carpentier.		“Contribution	a	l’Etude	du	Dispatching	Economique.”	Bull.	Soc.	

Francaise	Elect.,	Vol.	3,	Pp.	431‐447,	1962.	

J.	Carpentier.		“Optimal	Power	Flows.”		International	Journal	of	Electrical	Power	&	

Energy	Systems,	Vol.	1,	Iss.	1,	Pp.	3‐15,	April	1979.	

A.	Castillo,	R.P.	O’Neill.	“Optimal	Power	Flow	Papers:	Paper	5.	Computational	

Performance	of	Solution	Techniques	Applied	to	the	ACOPF.”	Federal	Energy	

Regulatory	Commission	ሺFERCሻ	website,	2013.	

E.D.	Castronuovo,	J.M.	Campagnolo,	R.	Salgado.	“On	the	Application	of	High	

Performance	Computation	Techniques	to	Nonlinear	Interior	Point	Methods.”	IEEE	

Transactions	on	Power	Systems,	Vol.	16,	No.	3,	August	2001.		

M.	Celis,	J.	E.	Dennis,	and	R.	A.	Tapia.	"A	Trust	Region	Strategy	for	Nonlinear	Equality	

Constrained	Optimization"	in	Numerical	Optimization	ሺP.	Boggs,	R.	Byrd	and	R.	

Schnabel,	edsሻ,	Philadelphia:	SIAM,	Pp.	71–82,	1984.	

R.M.	Chamberlain,	M.J.D.	Powell,	C.	Lemarechal,	H.C.	Pedersen.		“The	Watchdog	

Technique	for	Forcing	Convergence	in	Algorithms	for	Constrained	Optimization.”	

Math.	Prog.	Study	16,	Pp.	1‐17,	1982.	
Page

43

A.M.	Chebbo,	M.R.	Irving,	N.H.	Dandachi.	“Combined	Active	Reactive	Dispatch.”	IEEE	

Proc.	Gener.	Transm.	Distrib.,	Vol.	142,	No.	4,	July	1995.		

H.‐D.	Chiang,	B.	Wang,	Q.‐Y.	Jiang.		“Applications	of	TRUST‐TECH	Methodology	in	

Optimal	Power	Flow	of	Power	Systems”	in	Energy	Systems:	Optimization	in	the	

Energy	Industry,	Chapter	13,	2009.	

C.Y.	Chung,	W.	Yan,	F.	Liu.	“Decomposed	Predictor‐Corrector	Interior	Point	Method	

for	Dynamic	Optimal	Power	Flow.”	IEEE	Transactions	on	Power	Systems,	Vol.	26,	

No.	3,	August	2011.	

A.	R.	Conn,	Nicholas	I.	M.	Gould,	Philippe	L.	Toint.	"Trust‐Region."	MPS‐SIAM	Series	

on	Optimization,	1987.	

N.	Demirovic,	S.	Tesnjak,	A.	Tokic.		“Hot	Start	and	Warm	Start	in	LP	Based	Interior	

Point	Method	and	Its	Application	to	Multiperiod	Optimal	Power	Flows.”		Power	

Systems	Conference	and	Exposition	ሺPSCEሻ,	2006.	

R.	Divi,	H.K.	Kesavan.		“A	Shifted	Penalty	Function	Approach	for	Optimal	Load‐Flow.”	

IEEE	Transactions	on	Power	Apparatus	and	Systems,	Vol.	PAS‐101,	No.	9,	September	

1982.	

H.W.	Dommel,	W.F.	Tinney.		”Optimal	Power	Flow	Solutions.”			IEEE	Transactions	on	

Power	Apparatus	and	Systems,	Vol.	PAS‐87,	No.	10,	October	1968.

A.V.	Fiacco,	G.P.	McCormick.	Nonlinear	Programming:	Sequential	Unconstrained	

Minimization	Techniques.	John	Wiley	&	Sons,	1968.	

R.	Fletcher,	M.J.D.	Powell.	“A	Rapidly	Convergent	Descent	Method	for	Minimization.”	

Computer	J.,	Vol.	6,	Pp.	163‐168,	1963.	

R.	Fletcher,	C.M.	Reeves.	“Function	Minimization	by	Conjugate	Gradients.”	Computer	

J,	Vol.	7,	Pp	149‐154,	1964.	

R.	Fletcher.	“Fortran	Subroutines	for	Minimization	by	Quasi‐Newton	Methods.”	

AERE	Report,	Harwell,	England,	1972.	

R.	Fletcher.	Practical	Methods	to	Optimization,	Volume	1:	Unconstrained	

Optimization.	John	Wiley	&	Sons,	1980.	

R.	Fletcher.	Practical	Methods	to	Optimization,	Volume	2:	Constrained	Optimization.	

John	Wiley	&	Sons,	1981.	

R.	Fletcher,	S.	Leyffer,	and	P.	Toint,	“A	Brief	History	of	Filter	Methods.”	ANL/MCS‐

P1372‐0906,	September	26,	2006	revised	October	9,	2006	

K.R.	Frisch.	“The	Logarithmic	Potential	Method	for	Convex	Programming.”	

Unpublished	manuscript.	Institute	of	Economics,	University	of	Oslo,	Oslo,	Norway,	

1955.	

P.	E.	Gill,	W.	Murray,	M.	A.	Saunders,	and	M.	H.	Wright.	“Some	Theoretical	Properties	

of	an	Augmented	Lagrangian	Merit	Function.”	Advances	in	Optimization	and	Parallel	

Computing	ሺP.	M.	Pardalos,	edሻ,	Pp.	101‐128,	North	Holland,	1992.	
Page

44

P.	E.	Gill,	W.	Murray,	M.	A.	Saunders,	and	M.	H.	Wright.	“Maintaining	LU	Factors	of	a	

General	Sparse	Matrix,	Linear	Algebra	and	its	Applications.”	Vol.	88/89,		Pp.	239‐

270,	1987.	

T.C.	Giras	and	S.N.	Talukdar.		“Quasi‐Newton	Method	for	Optimal	Power	Flows.”		

International	Journal	of	Electrical	Power	and	Energy	Systems,	Vol.	3,	No.	2,	Pp.	59‐

64,	April	1981.	

J.	Gondzio.	“Multiple	Centrality	Corrections	in	a	Primal‐Dual	Method	for	Linear	

Programming.”	Computational	Optimization	and	Applications,	Vol.	6,	Pp.	137‐156,	

1996.	

S.	Granville.	“Optimal	Reactive	Dispatch	Through	Interior	Point	Methods.”		IEEE	

Transactions	on	Power	Systems,	Vol.	9,	No.	1,	February	1994.	

H.H.	Happ.		“Optimal	Power	Dispatch	–	A	Comprehensive	Survey.”	IEEE	Transactions	

on	Power	Apparatus	and	Systems,	Vol.	PAS‐96,	No.	3,	May/June	1977.	

M.	R.	Hestenes.	“Multiplier	and	Gradient	Methods.”	J.	Optim.	Theory	Appl.,	Vol.	4,	Pp.	

303–320,	1969.	

N.J.	Higham.	“Computing	a	Nearest	Symmetric	Positive	Semidefinite	Matrix.”	Lin.	

Alg.	Appl.,	Vol.	103,	Pp.	103‐118,	1988.		

M.	Huneault,	F.D.	Galiana.	“An	Investigation	of	the	Solution	to	the	Optimal	Power	

Flow	Problem	Incorporating	Continuation	Methods.”	IEEE	Transactions	on	Power	

Systems,	Vol.	5,	No.	1,	February	1990.	

K.	Iba.	“Reactive	Power	Optimization	by	Genetic	Algorithm.”	IEEE	Transactions	on	

Power	Systems,	Vol.	9,	No.	2,	Pp.	685‐692,	1994.	

R.	A.	Jabr,	A.	H.	Coonick,	B.	J.	Cory.	“A	Primal‐Dual	Interior	Point	Method	for	Optial	

Power	Flow	Dispatching.”	IEEE	Transactions	on	Power	Systems,	Vol.	17,	No.	3,	

August	2002.	

R.A.	Jabr.	“A	Conic	Quadratic	Format	for	the	Load	Flow	Equations	of	Meshed	

Networks.”	IEEE	Transactions	on	Power	Systems,	Vol.	22,	No.	4,	November	2007.	

R.A.	Jabr.	“Optimal	Power	Flow	Using	an	Extended	Conic	Quadratic	Formulation.”	

IEEE	Transactions	on	Power	Systems,	Vol.	23,	No.	3,	August	2008.	

Q.Y.	Jiang,	H.‐D.	Chiang,	C.X.	Guo,	Y.J.	Cao.	“Power‐Current	Hybrid	Rectangular	

Formulation	for	Interior‐Point	Optimal	Power	Flow.”	IET	Gener.	Transm.	Distrib.,	

Vol.	3,	Iss.	8,	Pp.	748‐756,	2009.	

Q.Y.	Jiang,	H.‐D.	Chiang.	“Sequential	Feasible	Optimal	Power	Flow:	Theoretical	Basis	

and	Numerical	Implementation.”	European	Transactions	on	Electrical	Power,	Vol.	

20,	Pp.	695‐709,	2010.	

Q.	Jiang,	G.	Geng,	C.	Guo,	Y.	Cao.	“An	Efficient	Implementation	of	Automatic	

Differentiation	in	Interior	Point	Optimal	Power	Flow.”	IEEE	Transactions	on	Power	

Systems,	Vol.	25,	No.	1,	February	2010.	

D.S.	Johnson.	“A	Theoreticians	Guide	for	Experimental	Analysis	of	Algorithms.”	In	

Proceedings	of	the	5th	and	6th	DIMACS	Implementation	Challenges,	2002.	
Page

45

N.	Karmarker.		“A	New	Polynomial‐Time	Algorithm	for	Linear	Programming.”	

Combinatorica,	Vol.	4,	Pp.	373‐395,	1984.	

L.	L.	Lai,	J.	T.	Ma.	“Improved	Genetic	Algorithms	for	Optimal	Power	Flow	under	both	

Normal	and	Contingent	Operation	States.”	Electrical	Power	&	Energy	Systems,	Vol.	

19,	No.	5,	Pp.	287‐292,	1997.	

J.	Lavaei,	S.	H.	Low.	“Zero	Duality	Gap	in	Optimal	Power	Flow	Problem.”	IEEE	

Transactions	on	Power	Systems,	Vol.	27,	Iss.	1,	Pp.	92‐107,	2012.	

S.	Leyffer,	A.	Mahajan.		Foundations	of	Constrained	Optimization.		Argonne	National	

Laboratory,	2010.	

F.	G.	M.	Lima,	S.	Soares,	A.	Santos	Jr.,	K.	C.	Almeida,	F.	D.	Galiana.	“Numerical	

Experiments	with	an	Optimal	Power	Flow	Algorithm	Based	on	Parametric	

Techniques.”	IEEE	Transactions	on	Power	Systems,	Vol.	16,	No.	3,	August	2001.	

C.‐H.	Lin,	S.‐Y.	Lin,	S.‐S.	Lin.	“Improvements	on	the	Daulity	Based	Method	Used	in	

Solving	Optimal	Power	Flow	Problems.”	IEEE	Transactions	on	Power	Systems,	Vol.	

17,	No.	2,	May	2002.	

S.‐Y.	Lin,	Y.‐C.	Ho,	C.‐H.	Lin.	“An	Ordinal	Optimization	Theory‐Based	Algorithm	for	

Solving	the	Optimal	Power	Flow	Problem	With	Discrete	Control	Variables.”	IEEE	

Transactions	on	Power	Systems,	Vol.	19,	No.	1,	February	2004.	

W.‐M.	Lin,	C.‐H.	Huang,	T.‐S.	Zhan.	“Predictor‐Corrector	Interior	Point	Algorithm	for	

Current‐Based	Optimal	Power	Flow.”	IEEE	Power	Engineering	Society	General	

Meeting,	Pp.	1‐6,	2007.	

W.‐M.	Lin,	C.‐H.	Huang,	T.‐S.	Zhan.	“A	Hybrid	Current‐Power	Optimal	Power	Flow	

Technique.”	IEEE	Transactions	on	Power	Systems,	Vol.	23,	No.	1,	February	2008.	

David	G.	Luenberger,	Linear	and	Nonlinear	Programming,	Addison‐Wesley,	1984.	

Olvi	L.	Mangasarian,	Nonlinear	Programming,	McGraw‐Hill,	1969.	

A.G.	Maria,	J.A.	Findlay.			“A	Newton	Optimal	Power	Flow	Program	for	Ontario	Hydro	

EMS.”		IEEE	Transactions	on	Power	Systems,	Vol.	2,	Iss.	3,	August	1987.	

S.	Mehrotra.	“On	the	Implementation	of	a	Primal‐Dual	Interior	Point	Method.”	SIAM	

Journal	on	Optimization,	Vol.	2,	Iss.	4,	Pp.	575‐601,	1992.	

A.P.S.	Meliopoulos,	F.	Xia.		“Optimal	Power	Flow	Application	to	Composite	Power	

System	Reliability	Analysis,”	IEEE/NTUA	Athens	Power	Tech	Conference,	Athens,	

1993.	

W.	Min,	L.	Shengsong.	“A	Trust	Region	Interior	Point	Algorithm	for	Optimal	Power	

Flow	Problems.”	Electrical	Power	and	Energy	Systems,	Vol.	27,	Pp.	293‐300,	2005.	

H.	Mittelman.	“AMPL‐NLP	Benchmark,	IPOPT,	KNITRO,	LOQO,	PENNLP,	SNOPT,	&	

CONOPT.”	Benchmarks	for	Optimization	Software.	Posted:	

http://plato.asu.edu/bench.html,	Last	Accessed:	January,	2013.	

J.A.	Momoh,	S.X.	Guo,	E.C.	Ogbuobiri,	R.	Adapa.		“The	Quadratic	Interior	Point	method	

Solving	Power	System	Optimization	Problems.”		IEEE	Transactions	on	Power	

Systems,	Vol.	9,	No.	3,	August	1994.	
Page

46

J.A.	Momoh,	J.Z.	Zhu.		“Improved	Interior	Point	Method	for	OPF	Problems.”		IEEE	

Transactions	of	Power	Systems,	Vol.	14,	No.	3,	August	1999.	

A.	Monticelli,	W‐H	E.	Liu.	“Adaptive	Movement	Penalty	Method	for	the	Newton	

Optimal	Power	Flow.”		IEEE	Transactions	on	Power	Systems,	Vol.	7,	No.	1,	February	

1992.	

P.	K.	Mukherjee,	R.	N.	Dhar.	“Optimal	Load‐Flow	Solution	by	Reduced‐Gradient	

Method.”	Proc.	IEEE,	Vol.	121,	No.	6,	June	1974.	

B.A.	Murtagh,	M.A.	Saunders.		“A	Projected	Lagrangian	Algorithm	and	its	

Implementation	for	Sparse	Nonlinear	Constraints.”		Mathematical	Programming	

Study,	Vol.	16,	Pp.	84‐117,	1982.	

B.A.	Murtagh,	M.A.	Sunders.	“MINOS	5.51	User’s	Guide.”	Stanford	University,	2003.	

J.	Nanda,	D.P.	Kothari,	S.C.	Srivastava.		“New	Optimal	Power‐Dispatch	Algorithm	

Using	Fletcher’s	Quadratic	Programming	Method.”	IEE	Proceedings	C:	Generation,	

Transmission	&	Distribution,	Vol.	136,	Iss.	3,	1989.			

I.M.	Nejdawi,	K.A.	Clements,	P.W.	Davis.		“An	Efficient	Interior	Point	Method	for	

Sequential	Quadratic	Programming	Based	Optimal	Power	Flow.”	IEEE	Transactions	

on	Power	Systems,	Vol.	15,	No.	4,	November	2000.	

J.		Nocedal,	S.	Wright.		Numerical	Optimization.		Springer‐Verlag,	2nd	Edition,	1999.	

J.	Nocedal.	“Conjugate	Gradient	Methods	and	Nonlinear	Optimization.”	in	Linear	and	

Nonlinear	Conjugate	Gradient‐Related	Methods	ሺL.M.	Adams,	J.L.	Nazareth,	edsሻ.	

Proceedings	of	the	AMS‐IMS‐SIAM	Summer	Research	Conference	held	at	the	

University	of	Washington,	July	1995.	

A.	R.	L.	Oliveira,	S.	Soares,	L.	Nepomuceno.	“Optimal	Active	Power	Dispatch	

Combining	Network	Flow	and	Interior	Point	Approaches.”	IEEE	Transactions	on	

Power	Systems,	Vol.	18,	No.	4,	November	2003.	

M.	Pai.	Energy	Function	Analysis	for	Power	System	Stability.	Kluwer	Academic	

Publishers,	Boston,	MA,	1989.	

S.	Pajic,	K.A.	Clements.		“Power	System	State	Estimation	via	Globally	Convergent	

Methods.”	IEEE	Transactions	on	Power	Systems,	Vol.	20,	No.	4,	November	2005.	

P.	Pardalos,	D.‐Z.	Du.	Optimization	Theory	and	Methods:	Nonlinear	Programming.	

Springer,	2006.		

D.T.	Phan.		“Lagrangian	Duality‐Based	Branch	and	Bound	Algorithms	for	Optimal	

Power	Flow.”		IBM	Research	Report.	April	2011.	

R.	Polyak.	“Modified	Barrier	Functions.”	Math.	Prog.,	Vol.	54,	No.	2,	Pp.	177‐222,	

1992.	

B.T.	Polyak.	“Newton’s	Method	and	Its	Use	in	Optimization.”	European	Journal	of	

Operational	Research.	Vol.	181,	Pp.	1086‐1096,	2007.	

R.A.	Ponrajah,	F.D.	Galiana.		“The	Minimum	Cost	Optimal	Power	Flow	Problem	

Solved	Via	the	Restart	Homotopy	Continuation	Method.”	IEEE	Transactions	on	

Power	Systems,	Vol.	4,	No.	1,	February	1989.	
Page

47

M.	J.	D.	Powell.	“A	Method	for	Nonlinear	Constraints	in	Minimization	Problems”	in	

Optimization	ሺR.	Fletcherሻ,	Academic	Press,	Pp.	283–298,	London	and	New	York	,	

1969.	

J.K.	Reid.	“Fortran	Subroutines	for	Handling	Sparse	Linear	Programming	Bases.”	

Report	R8269.	Atomic	Energy	Research	Establishment.	Harwell,	England,	1976.	

J.K.	Reid.	“A	Sparsity‐Exploiting	Variant	of	the	Bartels‐Golub	Decomposition	for	

Linear	Programming	Bases.”	Mathematical	Programming	Vol.	24,	Pp	55‐69,	1982.	

S.	M.	Robinson.	“A	Quadratically	Convergent	Algorithm	for	General	Nonlinear	

Programming	Problems.”	Mathematical	Programming	3,	Pp.	145‐156,	1972.	

J.B.	Rosen.	“The	Gradient	Projection	Method	for	Nonlinear	Programming.	Part	II.	

Nonlinear	Constraints.”	Journal	of	the	Society	for	Industrial	and	Applied	

Mathematics,	Vol.	9,	No.	4,	Pp.	514‐532,	Dec.	1961.	

R.	Salgado,	A.	Brameller,	P.	Aitchison.	“Optimal	Power	Flow	Solutions	Using	the	

Gradient	Projection	Method.	Part	2:	Modelling	of	the	Power	System	Equations.”	IEE	

Proceedings,	Vol.	137,	Pt.	C,	No.	6,	November	1990.	

A.	Santos,	Jr.,	S.	Deckmann,	S.	Soares.		“A	Dual	Augmented	Lagrangian	Approach	for	

Optimal	Power	Flow.”		IEEE	Transactions	on	Power	Systems,	Vol.	3,	No.	3,	August	

1988.	

A.M.	Sasson.	“Combined	Use	of	the	Powell	and	Fletcher‐Powell	Nonlinear	

Programming	Methods	for	Optimal	Load	Flows.”		IEEE	Transactions	on	Power	

Apparatus	and	Systems,	Vol.	PAS‐88,	No.	10,	October	1969.	

A.M.	Sasson,	F.	Viloria,	F.	Aboytes.		“Optimal	Load	Flow	Solution	Using	the	Hessian	

Matrix.”	IEEE	Transactions	on	Power	Apparatus	and	Systems.		Vol.	PAS‐92,	Iss.	1,	

1973.	

R.	R.	Shoults,	D.	T.	Sun.	“Optimal	Power	Flow	Based	Upon	P‐Q	Decomposition.”	IEEE	

Transactions	on	Power	Apparatus	and	Systems,	Vol.	PAS‐101,	No.	2,	February	1982.	

A.	A.	Sousa,	G.	L.	Torres.	“Globally	Convergent	Optimal	Power	Flow	by	Trust‐Region	

Interior‐Point	Methods.”	IEEE	Lausanne	Power	Tech,	2007.	

V.A.	de	Sousa,	E.C.	Baptista,	G.R.M.	da	Costa.		“Loss	Minimization	by	the	Predictor‐

Corrector	Modified	Barrier	Approach.”		Electric	Power	Systems	Research,	Vol.	79,	

Pages	803‐808,	2009.	

A.A.	Sousa,	G.L.	Torres,	C.A.	Cañizares.		“Robust	Optimal	Power	Flow	Solution	Using	

Trust	Region	and	Interior‐Point	Methods.”		IEEE	Transactions	on	Power	Systems,	

Vol.	26,	No.	2,	May	2011.	

B.	Stott,	O.	Alsac.	“Optimal	Power	Flow—Basic	Requirements	for	Real‐Life	Problems	

and	their	Solutions.”	Posted:	http://www.ieee.hr/_download/repository/Stott‐

Alsac‐OPF‐White‐Paper.pdf,	July	1,	2012.	Last	Accessed:	January,	2013.		

D.I.	Sun,	B.	Ashley,	B.	Brewer,	A.	Hughes,	W.F.	Tinney.	“Optimal	Power	Flow	by	

Newton	Approach.”	IEEE	Transaction	on	Power	Apparatus	and	Systems,	Vol.	PAS‐

103,	No.	10,	October	1984.	
Page

48

J.	E.	Tate,	T.	J.	Overbye.	“A	Comparison	of	the	Optimal	Multiplier	in	Polar	and	

Rectangular	Coordinates.”	IEEE	Transactions	on	Power	Systems,	Vol.	20,	No.	4,	

November	2005.		

G.L.Torres,	V.H.	Quintana.		“Nonlinear	Optimal	Power	Flow	by	a	Noninterior‐Point	

Method	Based	on	Chen‐Harker‐Kanzow	NCP‐functions.”		In	Proceeding	of	Electrical	

and	Computer	Engineering,	Vol.	2,	1998.	

G.	L.	Torres,	V.	H.	Quintana.	“Optimal	Power	Flow	by	a	Nonlinear	Complementarity	

Method.”	IEEE	Transactions	on	Power	Systems,	Vol.	15,	No.	3,	August	2000.	

G.L.	Torres,	V.H.	Quintana.	“On	a	Nonlinear	Multiple‐Centrality‐Correction	Interior‐

Point	Method	for	Optimal	Power	Flow.”	IEEE	Transactions	on	Power	Systems,	Vol.	

16,	No.	2,	May	2001.	

G.	L.	Torres,	V.	H.	Quintana.	“A	Jacobian	Smoothing	Nonlinear	Complementarity	

Method	for	Solving	Nonlinear	Optimal	Power	Flows.”	14th	PSCC,	Sevilla,	24‐28	June	

2002.	

G.L.	Torres,	M.A.	De	Carvalho.	”On	Efficient	Implementation	of	Interior‐Point	Based	

Optimal	Power	Flows	in	Rectangular	Coordinates.”	IEEE	Power	Systems	Conference	

and	Exposition,	Pp.	1747‐1752,	2006.	

University	of	Washington	ሺUWሻ	Electrical	Engineering.	Power	Systems	Test	Case	

Archive.	Posted:	http://www.ee.washington.edu/research/pstca/.	Last	accessed:	

January,	2013.	

L.S.	Vargas,	V.H.	Quintana,	A.	Vannelli.		“A	Tutorial	Description	of	an	Interior	Point	

Method	and	Its	Application	to	Security‐Constrained	Economic	Dispatch.”	IEEE	

Summer	Meeting,	Paper	SM	92	416‐8	PWRS,	1992.	

H.	Wang,	C.E.	Murillo‐Sánchez,	R.D.	Zimmerman,	R.J.	Thomas.	“On	Computational	

Issues	of	Market‐Based	Optimal	Power	Flow.”	IEEE	Transactions	on	Power	Systems,	

Vol.	22,	No.	3,	August	2007.	

A.	Wächter	and	L.	Biegler,	“Line	Search	Filter	Methods	for	Nonlinear	Programming:	

Motivation	and	Global	Convergence.”	SIAM	Journal	on	Optimization,	Vol.	16,	No.	1,	

Pp.	1‐31,	2005.	

H.	Wei,	H.	Sasaki,	J.	Kubokawa,	R.	Yokoyama.	“An	Interior	Point	Nonlinear	

Programming	for	Optimal	Power	Flow	Problems	with	A	Novel	Data	Structure.”		IEEE	

Transactions	on	Power	Systems,	Vol.	13,	No.	3,	August	1998.	

P.	Wolfe.		“The	Simplex	Method	for	Quadratic	Programming,”	Econometrica,	Vol.	27,	

Pp.	382‐398,	1959.	

P.	Wolfe.	“Methods	of	nonlinear	programming.”	Nonlinear	Programming,	J.	Abadie,	

Ed.,	North	Holland,	Amsterdam,	Pp.	97‐131,	1967.	

Y.‐C.	Wu,	A.S.	Debs,	R.E.	Marsten.		“A	Direct	Nonlinear	Predictor‐Corrector	Primal‐

Dual	Interior	Point	Algorithm	for	Optimal	Power	Flows.”		IEEE	Transactions	on	

Power	Systems,	Vol.	9,	No.	2,	May	1994.	
Page

49

F.F.	Wu,	G.	Gross,	J.F.	Luini,	P.M.	Look.		“A	Two‐Stage	Approach	to	Solving	Large‐

Scale	Optimal	Power	Flows.”		IEEE	Conference	Proceedings	,	Power	Industry	

Computer	Applications	Conference,	PICA‐79,	Pp.	126‐136,	1979.	

L.	Xie,	H.‐D.	Chiang.		“A	Enhanced	Multiple	Predictor‐Corrector	Interior	Point	

Method	for	Optimal	Power	Flow.”		IEEE	Power	and	Energy	Society	General	Meeting,	

2010.	

L.	Xie,	H.‐D.	Chiang.	“Weighted	Multiple	Predictor‐Corrector	Interior	Point	Method	

for	Optimal	Power	Flow.”	Electric	Power	Components	and	Systems,	Vol.	39,	Pp.	99‐

112,	2011.	

K.	Xie,	Y.H.	Song.		“Dynamic	Optimal	Power	Flow	by	Interior	Point	Methods.”	IEE	

Proceedings:	Generation,	Transmission,	&	Distribution,	Vol.	148,	Iss.	1,	Pp.	76‐84,	

2001.		

X.	Yan,	V.H.	Quintana.		“Improving	An	Interior‐Point‐Based	OPF	by	Dynamic	

Adjustments	of	Step	Sizes	and	Tolerances.”		IEEE	Transactions	on	Power	Systems,	

Vol.	14,	No.	2,	May	1999.	

B.‐Y.	Yang,	J.‐C.	Peng.	“A	Novel	Optimal	Power	Flow	Model	and	Algorithm	Based	on	

the	Imputation	of	an	Impedance‐Branch	Dissipation	Power.”	Asian‐Pacific	Power	

and	Energy	Engineering	Conference,	2009.	

R.	Zhou,	Y.	Zhang,	H.	Yang.		“A	Trust‐Region	Algorithm	Based	on	Global	SQP	for	

Reactive	Power	Optimization.”		2nd	International	Conference	on	Electrical	and	

Electronics	Engineering	ሺICEEEሻ	and	XI	Conference	on	Electrical	Engineering	ሺCIE	

2005ሻ.		Mexico	City,	Mexico.		September	2005.	

R.D.	Zimmerman,	C.	E.	Murillo‐Sanchez,	R.J.	Thomas.	“MATPOWER:	Stead‐State	

Operations,	Planning,	and	Analysis	Tools	for	Power	Systems	Research	and	

Education.”	IEEE	Transactions	on	Power	Systems,	Vol.	26,	No.	1,	February	2011.