Page 229 - Numerical Methods for Chemical Engineering

218 5 Numerical optimization

In the conjugate gradient method, this zig-zag behavior is reduced by mixing-in a portion
of the previous search direction such that the trajectory tends not to double back upon itself,
p [k+1] =−γ [k+1] + β [k] [k] (5.13)
p
β [k+1] is commonly chosen by one of the following two formulas:
γ · γ
 [k+1] [k+1]

 , Fletcher–Reeves (CG-FR)
γ · γ
 [k] [k]



β [k] = (5.14)
 [k+1] [k+1] [k]
 γ · γ − γ


, Polak–Ribiere (CG-PR)


γ [k] · γ [k]
The CG-FR formula yields search directions with very favorable convergence properties for
quadratic cost functions, as is explained below. The CG-PR formula gives the same search
direction as CG-FR for quadratic cost functions, but its extra term biases the search direction
towards the direction of steepest descent when the cost function is far from quadratic. As
steepest descent is more robust in this case, CG-PR is preferred.
A gradient minimizer routine
gradient minimizer.m implements either the CG-PR method as the default or the steepest
descent method. It uses a weak line search starting from an initial step based on quadratic
approximation. Thus, when the cost function is (locally) quadratic, a strong line search is
conducted. The syntax is

[x, F, grad, iflag, x traj] = gradient minimizer( . . .
func name, x0, OptParam, ModelParam);

func name is the name of a routine that returns F(x) and γ(x),
function [F, grad, iOK]= func name(x, ModelParam);

Here, iOK is 1 if the cost function was evaluated correctly.
ModelParam is a structure that contains the fixed parameters of the cost function. x0 is
the initial guess, and OptParam is an optional structure that controls the behavior of the
minimizer. As output, x is the local minimum, F and grad are its cost function and gradient,
iflag is 1 if the minimizer converged. x traj is an optional array recording the progress of
the minimizer.
The steepest descent and gradient methods are demonstrated in Figure 5.6 and Figure 5.7
for the cost function
4
2
2
F(x) = (x 1 − 1) + 10(x 2 − 2) + c(x 1 − 1) + c(x 2 − 2) 4 (5.15)
The following routine returns the cost function and its gradient:
function [F, grad, iOK] = simple cost func(x, ModelParam);
iOK = 0; c = ModelParam.c;
dx1=x(1)-1;dx2=x(2)-2;