Transcript PowerPoint 簡報

Chapter 4
Sensitivity Analysis, Duality and
Interior Point Methods
Interior Point Methods
•Simplex methods reach the optimal solution by traversing a series
of basic solutions. The path followed by the algorithm moves
around the boundary of the feasible region.
•This can be cripplingly inefficient given that the number of
extreme points grows exponentially with n for m fixed.
•The running time of an algorithm as a function of the problem
size is known as its computational complexity.
• In practice, the simplex method works surprisingly well,
exhibiting linear complexity—i.e., its running time is proportional
to m + n.
•Researchers, however, have long tried to develop solution
algorithms whose worst-case running times are polynomial
functions of the problem size.
Interior Point Methods (IPMs)
•First Success: Russian Mathematician Khachian
proposed the ellipsoid method with complexity O(n6).
It is efficient only theoretically, no one ever realize an
implementation that matched the performance of current
simplex codes.
•Big Success: Karmarkar(1984)
IPMs have produced solutions to many industrial
problems that previously were intractable.
Three major types of IPMs:
(1) Potential Reduction algorithm: most closely
embodies the constructs of Karmarkar;
(2) Affine Scaling algorithm: the simplest to
implement;
(3) Path Following algorithms: combine excellent
behavior in theory and practice.
•
We discuss a member of the third category
known as the primal-dual path following
algorithm. The algorithm operates
simultaneously on the primal and dual linear
programming problems.
Example:
Primal model
Dual model
Max ZP = 2X1+ 3 X2
Min ZD = 8π1 + 6π2
s.t.
s.t.
2 X1 + X2 + X3 =8
X1 + 2X2
+ X4 = 6
Xj≧0, j=l,.,.,4
2π1 +π2 - Z1
π1 + 2π2
π1
=2
- Z2 =3
-Z3 =0
π2
-Z4 =0
Zj≧0, j=1,...,4
Note: Zj’s denote the slack variables of the dual model.
• The algorithm iteratively modifies the primal and
dual solutions until the optimality conditions of
primal feasibility, dual feasibility, and
complementary slackness are satisfied to a sufficient
degree. This event is signaled when the duality gap—
the difference between the primal and dual objective
functions—is sufficiently small.
• For the purposes of the discussion, "interior of the
feasible region" is taken to mean that the values of
the primal and dual variables are always strictly
greater than zero.
The Primal and Dual Problems
Basic Ideas:
• The application of the Lagrange multiplier method to the
transformation of an equality constrained optimization
problem into an unconstrained one.
• The transformation of an inequality constrained optimization
problem into a sequence of unconstrained problems by
incorporating the constraints in a logarithmic barrier function
that imposes a growing penalty as the boundary { Xj = 0, and
Zj = 0 for all j) is approached.
• The solution of a set of nonlinear equations using Newton's
method, thereby arriving at a solution of the unconstrained
optimization problem.
The Lagrangian Approach
Consider the general problem
Maximize f(x)
subject to gi (x) = 0, i= 1,...,m
where f(x) and gi (x) are scalar functions of the ndimensional vector x.
The Lagrangian for this problem is
m
L(x,π)=f(x)-   i g i ( x)
i 1
where the variables π = ( π1 , π2 ,..., πm ) are the
Lagrange multipliers
Necessary conditions for a stationary point (maximum
or minimum):
L
, j =1,...,n and
x j
L
 0,
i
i  1, 2,..., m
• For linear constraints (aix –bi = 0) , the conditions
are sufficient for a maximum if the function f(x) is
concave and sufficient for a minimum if f(x) is
convex.
The Barrier Approach
• Considering the non-negativity conditions, the idea of the
barrier approach, as developed by Fiacco and McCormick
[1968], is to start from a point in the strict interior of the
inequalities (xj > 0, and zj > 0 for all j) and construct a barrier
that prevents any variable from reaching a boundary (e.g.,
xj = 0).
• Adding log(xj) to the objective function of the primal
problem, for example, will cause the objective to decrease
without bound as Xj approaches 0.
• The difficulty is that if the constrained optimal solution is
on the boundary (that is, one or more Xj*= 0, which is
always the case in linear programming), then the barrier
will prevent us from reaching it. To get around this
difficulty, we use a barrier parameter μ that balances the
contribution of the true objective function with that of the
barrier term.
• The parameter μis required to be positive and controls the
magnitude of the barrier term. Because function log(x)
takes on very large negative values as x approaches zero
from above, as long as the problem variables x and z
remain positive the optimal solutions to the barrier
problems will be interior to the nonnegative orthants (Xj >
0 and Zj > 0 for all j).
• The barrier term is added to the objective function for a
maximization problem and subtracted for a minimization
problem.
• The new formulations have nonlinear objective functions
with linear equality constraints, and can be solved with the
Lagrangian technique for μ > 0 fixed. The solutions to
these problems will approach the solution to the original
problem as μ approaches zero.
Note:
• These conditions are also sufficient, because the primal
Lagrangian is concave and the dual Lagrangian is convex.
• The optimal conditions are nothing more than primal
feasibility, dual feasibility, and complementary slackness
satisfied to within a tolerance of μ.
To facilitate the presentation, we define two n × n diagonal
matrices containing the components of x and z,
respectively—i.e.,
X = diag{ x1, x2,.., xn }
Z = diag{ z1, z2,..., zn }
• Let e = (1,1,..., 1)T be a column vector of size n. The
necessary and sufficient conditions can be written as
Primal feasibility:
Ax - b = 0
(m linear equations)
Dual feasibility: A T πT - z - c T = 0
(n linear equations)
μ-Complementary slackness: XZe - μe = 0 (n nonlinear
equations)
We must now solve this set of nonlinear equations for the
variables (x, π, z).
Stationary Solutions Using Newton's Method
• Consider the single variable problem of finding y to satisfy
the nonlinear equation
f(y)=0
where f is once continuously differentiable. Let y* be the
unknown solution.
• Taking the first-order Taylor series expansion of f(y)
around yk yields
f ( y k 1 )  f ( y k )  f ' ( y k )
where △= yk+1– yk.
• Setting the approximation f(yk+1) equal to zero and solving
for △ yields
k
 f (y )
 ' k
f (y )
• The point yk+1 = yk + △ is an approximate solution
of the equation.
• It can be shown that if one starts at a point y°
sufficiently close to y*, then yk→y* as k→∞.
• Newton's method can be extended to multidimensional
functions. Consider the general problem of finding the
r-dimensional vector y that solves the set of r equations
fi (y)=0, i=l,...,r or
f(y)=0
• The Jacobian at yk is
 
J yK
 f1

 y1
 f 2
  y
 1
：
 f r
 y
 1
f1
y 2
f 2
y 2
：
f r
y 2




f1
y r
f 2
y r
：
f r
y r










All the partial derivatives are evaluated at yk.
• Taking the first-order Taylor series expansion around the
point yk and setting it equal to zero yields
f(yk)+ J(yk)d=0
where d = y k+1 - y k is an r-dimensional vector whose
components represent the change in position for the
(k+l)th iteration.
• Solving for d leads to
d= -J(yk) -1f(yk)
• The point yk+1 = yk + d is an approximation of the
solution of the set of equations.
Newton's Method for Solving Barrier Problems:
• We now use Newton's method to solve the optimality
conditions for the barrier problems for a fixed value of μ.
• For y = (x, π, z) and r = 2n + m, the corresponding
equations and Jacobian are
Ax  b  0

A


T T
T
A   z  c  0  J ( y )   0
Z
XZe  e  0 

0
AT
0
0 

I
X 
where e is an n-dimensional column vector of 1's.
• Assume that we have a starting point (x°, π°, z°) satisfying
x° > 0 and z° > 0, and denote by
δp=b-Ax°
δD= cT –AT(π°)T+ z°
• The optimality conditions can be written as
J(y)d = -f(y)
A

0
Z

0
AT
0
0  d x    p

  

 I  d     D

X  d z   e  XZe 
where the (2n + m)-dimensional vector d = (dX , dπ, dz)T is
called the Newton direction. We must now solve for the
individual components of d.
• In explicit form, the preceding system is
Adx   p
AT d  d z   D
Zd x  XdZ  e  XZe
• The first step is to find dπ. With some algebraic manipulation,
we get
( AZ 1 XAT )d  b  AZ 1e  AZ 1 X D
or
d  ( AZ 1 XAT )1 (b  AZ 1e  AZ 1 X D )
Note that Z-l = diag{ 1/z1, l/z2,.. , 1/zn} and is trivial to compute.
• Further multiplications and substitutions lead to
dz= -δD+A T dπ
and
dx = Z-l (μe-XZe-Xdz)
• From these results, we can see in part why it is necessary to
remain in the interior of the feasible region. In particular, if
either Z-l or X-l does not exist, the procedure breaks down.
• We move from the current point (x°, π°, z°) to a new point
(x1, π1, z1) by taking a step in the direction (dx, dπ , dz). The
step sizes αp and αD are chosen to preserve x > 0 and π > 0.
• This requires a ratio test similar to that performed in the
simplex algorithm. The simplest approach is to use
  x j

 p   min 
: (d x ) j  0 
j

 (d x ) j
  z j

 D   min 
: (d z ) j  0 

 (d z ) j
where γ is the step size factor that keeps us from actually
touching the boundary. Typically, γ= 0.995. The notation (dx) j
refers to the jth component of the vector dx
• The new point is
x1  x 0   p d x
 1   0   D d
z1  z 0   D d z
which completes one iteration.
• Ordinarily, we would now resolve using (x1, π1, z1) to find
a new Newton direction and hence a new point. Rather
than iterating in this manner until the system converges for
the current value of μ, it is much more efficient to reduce μ
at every iteration.
• It is straightforward to show that the Newton step reduces
the duality gap—the difference between the dual and
primal objective values at the current point. Assume that x°
is primal feasible and that (π0, z°) is dual feasible.
• Let "gap(0)" denote the current duality gap.
gap(0) =π0b - cx°
=π0(Ax0) - (π0A - z0) T x°
(primal and dual feasibility)
= (z0) T x°
• If we let α = min{ αp, αD}, then
gap( )  ( z 0  d z )T ( x0  d z )T
• It can be shown that gap(α) < gap(0) as long as
gap (0)

n
In our computations, we have used
k
k T k
gap
(

)
(
z
) x
k
 

2
n
n2
which indicates that the value of μk is proportional to the duality
gap.
• Termination occurs when the gap is sufficiently small—say, less
than 10-8.
Primal-Dual Interior Point Algorithm
Inputs:
1. The data of the problem (A, b, c), where the
m x n matrix A has full row rank
2. Initial primal and dual feasible solutions x°
> 0, z° > 0, π°
3. The optimality tolerance ε > 0 and the step
size parameter γ (0,1)

Step 1: (Initialization) Start with some feasible
point x° > 0, z° > 0, π° and set the iteration
counter k = 0.
Step 2: (Optimality Test) If (zk) Txk < ε, stop;
otherwise, go to Step 3.
Step 3: (Compute Newton Directions) Let

 k k
( z k )T x K
 k k
k
k
k
k
X  diag  x1 , x2 ,..., xn  , Z  diag  z1 , z2 ,..., zn  ,  
2
n



k
Solve the following linear system to get dxk , dk , and d zk
Ad  0
k
x
A d d  0
k
k
Zd x  Xd z  e  XZe
T
k

k
z
Note that δP=0 and δD=0 as a result of the feasibility
of the initial point.
Step 4: (Find Step Lengths) Let
k
k



z


  xj

j
k
k
 P   min k : (d x ) j  0 ,  D   min k : (d z ) j  0 
j
j
(
d
)




x
j

 (d z ) j
Step 5: (Update Solution) Take a step in the Newton direction to get
x

k 1
 x  Pd
k
k
x
K
    D d
k 1
k
k
z  z   D d
k 1
k
Put k← k + 1 and go to Step 2.
Example:
We have n = 4 and m = 2, with data
 2 1 1 0
8
, b   
c  (2,3,0,0), A  
1 2 0 1
 6
and variable definitions
 x1 
 z1 
 
 
 x2 
 z2 
x   , z   ,    1 ,  2 
x3
z3
 
 
x 
z 
 4
 4
For an initial interior solution, we take x° = (1,1, 5, 3), z° = (4, 3, 2,
2), and π° = (2, 2) .
•The primal residual vectorδP and the dual residual vector δD are
both zero, since x° and z° are interior and feasible to their respective
problems.
•The column labeled "Comp. slack" is the μ-complementary
slackness vector computed as
XZe-μe
Since this vector is not zero, the solution is not optimal.