GREEN’S FUNCTIONS AND OPTIMAL SYSTEMS The Gradient Direction in Decision Space M. M. DENN’AND RUTHERFORD ARlS Department
of
Chemical Engineering, University
of Minnesota.
Minneapolis, M i n n .
The gradient direction in the space of decision vectors for multistage optimization problems is obtained under general constraints through the? solution of variational equations b y means of Green’s functions. This direction i s then used for a steepest ascent solution of the optimization problem.
N THE
determination of an optimal policy in a multistage
I-decision problem one of the most serious difficulties is the “curse of dimensionality,” by which we mean that the computation time has an exponential dependence on the number of state variables. Through the solution of the first-order variational equations of a system by means of Green’s functions we have shown ( 7 7) how successive “approximations to the problem” may be generated from the weak maximum principle (70, 7 7 , 27) and the dimensionality difficulties overcome. In this paper we continue our treatment of the Green’s function as a unified approach to optimization problems and show how the solution of the variational equations may be used to develop another computational method which overcomes the curse by a technique for which Dreyfus has coined the name “approximations to the solution” (72). In this scheme we determine the gradient direction in the space spanned by the set of decision vectors and use the method of steepest ascent to obtain the solution. Method of Steepest Ascent
The development of the method of steepest ascent is due to Cauchy (5) and perhaps independently to Sarrus (22). In its simplest form, if we seek to maximize a function P of a vector x then, given an initial value xo and estimate of the maximum P ( x o ) , we seek the set of direction cosines in x space which will result in the greatest rate of increase of P and then move in that direction to a new value xl--that is, if the x space is Euclidean and is parametrized by a distance variable, s, choose dxlds. subject to dxi dxj ai,, - - = 1 (1) ds ds in order to maximize dP/ds. (The summation convention, in which summation is carried out over all values of an index repeated once as a subscript and once as a superscript, is used throughout. ai, is Kronecker’s delta, which is equal to unity if i = j and zero otherwise.) The classical method (78) is to introduce a Lagrange multiplier, Xo, and to seek stationary points with respect to
1
Present address, University of Delaware, Newark, Del
Since
dP b P dxi ds a x i ds
(3)
it easily follows that the optimal choice of direction cosines is
which is simply the direction of the gradient vector. x is then changed by some amount in this direction and the process repeated until the maximum is approximately reached. The choice of a Euclidean metric is arbitrary, and for many physical problems the natural metric is non-Euclidean (23, 2 4 ) . Modifications upon the basic scheme are then necessary in order to obtain convergence within a reasonable amount of computing time. Wilde (24) has discussed such procedures, and perhaps the most powerful technique now available when the gradient direction is known analytically is the method described by Fletcher and Powell ( 7 3 ) . The calculation of the gradient direction in a function space for the solution of variational problems was first suggested by Hadamard ( 6 ) , and the method has been used with great success in the determination of optimal trajectories in continuous systems. Horn and Troltenier (76) have used this technique in the optimization of tubular chemical reactors. Kelley has written an excellent introduction to the method in a recent article (78) which contains a large number of references, and somewhat different points of view have been adopted by Dreyfus (721, Ho and Brentani ( 7 4 , and Denham and Bryson (3, 9 ) . In contrast to the extensive studies of gradient methods for continuous systems, the application of such techniques to multistage decision problems with discrete stages has been extremely limited, the significant exception being the work of Horn and Troltenier (77). Recently, Lee ( 7 9 ) has extended the work of Dreyfus on discretized continuous systems to unconstrained staged systems and studied the optimum solvent distribution in crosscurrent extraction. In this paper we show how the gradient direction in decision space for staged systems with a straight-chain topology may be obtained under general constraints on both the decision and state spaces. VOL. 4
NO. 2
MAY
1965
213
For completeness we also obtain an extension of the simplest result for continuous systems to systems with lag time. I n a later communication we will discuss the extension of these and previous results to mixed systems of general topology.
The state of the system produced by the (n - 1)st stage is denoted by the S-dimensional vector pn- 1. This is transformed by the action of stage n in a manner described by the transformation equation =
Tn(pn-1, q n ) ,
n = 1, 2, ' .
' >
N
bT n k ~
an,
bPn-1
(15)
(16) Equation 13 may then be written in an alternate form which is often convenient,
(5)
where the basis of the p space is Cartesian and qn is the vector of R, decisions which may be made at stage n. The optimal N stage policy is the set { q n ) which is admissible under the constraints of the problem and causes the prescribed objective function P ( P , ~ to ) take on a maximum value. If we specify a particular set of decision vectors, { qn),and a particular initial condition, PO,then application of Equation 5 defines a trajectory { &). The first-order variational equations about this trajectory as a result of small changes in ( q n )and P O are then
where 6qn is the change from the decision vector Zin, 6pc the change from Po, 6pn the corresponding variation in p n , and all partial derivatives are evaluated along the nominal trajectory {pni. Equation 6 is linear and may be solved by application of Green's identity ( 7 7) in the form
The Green's tensor,
=
where the Green's vector, a, is related to the Green's tensor by
Variational Equations
Pn
a,
Gradient Direction
Let us suppose that the metric in the space of decision vectors is such that distance is defined by h'
ds2 =
gntjdqnzdqn'
(18)
n=O
or
gnt3is the positive definite metric tensor at each stage, which serves both as a normalizing factor for decisions with different dimensions and as a weighting factor. For the special case of a Euclidean metric gnil = 6 t 3 . but in general it will be necessary to weight some decisions more heavily than others in order to prevent oscillations and speed convergence. We seek the maximum of Equation 13 subject to the constraint Equation 19. If we introduce the Lagrange multiplizr, 10,and seek the unconstrained stationary points of
anji,satisfies the difference equations then it follows at once that the gradient direction is
and the boundary condition
aN,i = 6ji
(9)
If p o is specified, then 6 p 0 vanishes. If we consider the feed p o as a decision vector which is to be chosen optimally, it is convenient to define a fictitious stage zero with transformation
gnzf is the associated metric tensor, whose components are formed from the inverse matrix of components of gn,,. The local direction of maximum rate of improvement in the profit is then given by Equation 21, or, alternatively,
and write Equation 7 as
If we then introduce aI distance measure, s, into the space spanned by (qn) we may divide Equation 11 by 6s and take the limit as 6s + 0 to obtain
anj =
~
b T 2 dqnk ~
dqnk
-
ds
(1 3)
Equation 8 is linear and homogeneous, so we may take its inner product with the vector bP/bpp,' to obtain the difference equation 214
l&EC FUNDAMENTALS
dqn' ds (The tensorial character of Equations 21 and 22 causes no difficulty, because we have identified the decision with the contravariant. rather than "physical" components of the vector.) Constraints on the decision space of the form
and the rate of change of the profit, P(p,) is
dP - bP dp.vN"- bP ds bPjvz ds bp,vz n
(21)
present no particular difficulty. If the restriction is not actually in force and inequality obtains, then the gradient direction is that defined by Equation 22. If we ultimately
reach a point in decision space where equality holds for one or more of the constraints (Equation 23) and the free gradient direction leads into a n inadmissible region, it is necessary to remain on the constraint surface for the next move. This is done by requiring that
dQnk ds
-
b Q n k 4,' bqnt ds
By substituting Equation 31 into Equation 30 we find that the X 1 must satisfy the linear equation
=o
Equation 24 may then be included as a constraint by means of Lagrange multipliers, At, and the projected gradient direction is
By slight modification of a proof by Greenspan in the appendix to (78) it follows that Equation 32 has a solution for 31 if the matrix of components
dqni ds
(33)
- - 3
is nonsingular and
(25) where it follows from Equation 24 that 3, satisfies the linear equation
has rank a t least L. The condition on Equation 33 is the Gram determinant criterion for independence of the constraints ( 8 ) . I t is more realistic to suppose, however, that rather than exactly satisfying the constraints (Equation 28), the nominal trajectory { p,] will be such that
In particular, for the important case that
where qnk is a constant, then qnk is simply maintained a t the value gnk. When more than one constraint surface is reached, all combinations must be tested to determine which constraints should be set a t equality.
where c 1 is some nonzero value. Instead of Equation 29, then, we seek to satisfy, to a linear approximation,
(36) where A, is the step length in the decision space. Proceeding as usual, we find that the gradient direction is given again by Equation 31, but the X 1 must satisfy the nonlinear equation
Constrained Terminal Points
Let us suppose that the effluent from the last stage must satisfy the L linearly independent constraints
where the gradient of Q l is nonvanishing when Equation 28 is satisfied. If we suppose initially that the nominal trajectory { F n } is such that the constraints are all satisfied, then we wish to maintain Q z a t a constant zero value while moving in the decision space. T h a t is, we require
2"aQL~PN' =o dPpNz ds
ds
or, using Equation 12,
bQi
N
I-
n =O
mn51
~PPN
bTn5 dqnk __ - -0 bqnk ds
We define L Lagrange multipliers, X 1 , and include the constraint (Equation 28) in the maximization of 13 by adding the inner product of 3L with Equation 30. I t then follows that the gradient direction projected onto the subspace spanned by Equation 28 is
/
Under the conditions placed on Expressions 33 and 34, Equation 31 will have a unique solution for tk = 0 and may be expected to have a solution in a neighborhood of ck = 0. In general, however, there will be values of ex for which Equation 37 cannot be solved for a, and one should not attempt a large correction in one step. If the constraints are inequalities, the methods of this section may be applied when the unconstrained gradients leads to an inadmissible region, but with more than one constraint all possibilities must be examined to determine which constraints are to be treated as equalities.
bP
VOL. 4
NO. 2
MAY 1965
215
Penalty Functions
.4n alternative method for incorporating end point constraints is by the use of penalty functions, as introduced by Courant (6, 7), in which we seek the unconstrained maximum of 1
J'(PN)
=
P(p.v) -
KijQ'(p,v)Q'(p.v)
(38)
where (K,,) must be positive definite and will usually be diagonal. For convenience we will assume that ( K i f ) is symmetric. It is clear fhat as the "penalty functions" K,, become arbitrarily large it is necessary for a finite solution that Qj --* 0, and hence a succession of problems is solved for increasing K i j until the constraints are satisfied within allowable limits. It can be shown that if the procedure converges, the result is the required constrained solution (7). The gradient direction for the problem defined by Equation 38 follows directly from Equation 21 as
When the maximum is reached, the derivatives of P in all directions must vanish, and in particular in the directions defined by the decision space grqdients of the constraints Q1. When we set the derivatives of P in those directions to zero, we obtain the linear equation
Equation 40 differs from Equation 32 only in sign. As K i j + m and the constraints are identically satisfied, the two equations will define identical functions a t the solution point. Thus X i = - lim (K,,Q')
(41)
In particular, when ( K t j ) is diagonal and K,, = K(,)aij, (no sum over i), then in the neighborhood of the solution to the constrained problem,
Kt A - X , / Q i
(42)
Kelley (78) has used this relation to suggest that successive values of K ( t ) be chosen by estimating an approximation x i from Equation 32 along the approximate solution determined by the form (Equation 38), and then choosing the new K ( i ) by
Ki
= ,it/;'
i0
(44)
IBEC FUNDAMENTALS
(46)
(47) By use of the transformation Equation 5 we may always express such constraints as
Qn'(pn-1, q n )
=
0
(48)
Suppose one such set of constraints is imposed at stage m (m < .\r). These may then be converted to end point con-
straints by defining the new state variables
and the theory of the two previous sections is applicable. Inequalities may be treated in the same way. Alternatively, problems with constraints of the form of Equation 48 may be solved by a method of projection which is equivalent to the method of Denham and Bryson for continuous systems ( 9 ) . \Vhen the constraint is satisfied along the nominal trajectory the variational equation is
If there are L constraints of the type of 48, then we may solve Equations 51 for the first L 6qnX. The variational Equation 6 then becomes
where Fka is the inverse matrix of the first L bQnk/&,j, summation over cy and k runs from 1 to L, and summation over p from L 1 to R,. The Green's tensor must then satisfy
+
and the gradient direction is described by Equation 21 or 22 for the R, - L free decisions, with the L additional relations (Equation'48). Inequalities may be treated in the same way when the free gradient direction extends into an inadmissible region. Continuous Systems with Lag Time
(45)
where So($?) is a Heaviside step function with the value zero for q f < 0 and unity for @f 2 0. The equality constraint de216
I t may frequently be necessary to consider constraints at stages other than the last in the form
Many physical systems may be represented by differencedifferential equations of the form
then we may follow Kelley (78) and define
Qj(p~= ) [@'(p.v)I' So(@')
Trajectory Constraints
(43)
where t i is the allowable error in Qi. Because we have not found it necessary to place a nonvanishing condition on the gradient of Q i ( p N ) in the penalty function treatment we may convert all inequality constraints on the endpoint to equality. If the trajectory must satisfy
@'(P'")
fined by Equation 45 is continuous, has continuous first derivatives, and is equivalent to the inequality constraint {Equation 44).
p
=
T [ p ( t ) ,p ( t - A ) , q ( t ) l ,
0
5t 5
0
(54)
where 8 is the total process duration and A is the "lag time"
in the system. tion”
I t is also necessary to specify an “initial func-
p(t) = p*(t),
5 t
