- Email: [email protected]
Finite-element solution of nonlinear optimal control problems with a quadratic performance index
Comp~r. chm. Engag, Vol. 1I, No. I, pp. 77-81, 1987 Printed in Great Britain. All rights reserved
0098-1354/87
Copyright 0
1987
$3.00 +o.oll
Pergamon Journals Ltd
SHORT NOTES FINITE-ELEMENT SOLUTION OF NONLINEAR OPTIMAL CONTROL PROBLEMS WITH A QUADRATIC PERFORMANCE INDEX C. KIPARI~SDFSand A. GEORCIOU~ Department of Chemical Engineering, University of Thessaloniki, Thessaloniki, Greece (Received 18 January 1984;/inal revision received 31 Juiy 1986; received for publication 22 August 1986)
Abstract-A new iterative method is developed
for the solution of nonlinear optimal control problems with a quadratic performance index. The new method employs a local linearization technique together with the finite-element method to obtain the optimal state and control sequences. Preliminary results show that the proposed method is simple, accurate and convergence is achieved in a few iterations. The new technique is not very sensitive to initial guesses and it is shown to perform at least as well as the classical techniques (i.e. steepest descent, quasilinearization).
INTRODUCTION The numerical methods that have been developed to solve the optimal control problem can be broadly classified into direct and indirect methods, Furthermore, they can be of iterative or noniterative type. Nakamichi and Washizu [l] developed a direct noniterative approach for solving the linear optimal
Consider any controlled process operating over the fixed interval, t e[tO, t,], which can be described by a system of nonlinear differential equations ir = f(x, II, t),
(1)
where x is the state vector of dimension n and u is the control vector of dimension m. Given the initial condition vector for x,
control problem with a quadratic performance index. The main advantage of their method is that the necessary conditions for the augmented performance index to be stationary can be expressed in the form of a system of algebraic equations. As a result, the optimal state and control trajectories are obtained by solving a system of linear algebraic equations. However, application of the finite-element method (FEM) to nonlinear problems results in a system of nonlinear
x(4l) = %9
(2)
the solution of equation (1) is uniquely defined, once the control vector u is chosen. Consider the performance index J, If L(x, UJ) dt J = G(x(t,), 4) + (3) s @I
algebraic equations which are hard to solve. This seems to mitigate against the effectiveness and simplicity of the technique for nonlinear problems. To circumvent this problem, we propose a new iterative approach to the solution of nonlinear control problems with a quadratic performance index using the FEM. According to the new method, a local linearization technique is used to obtain a linear model, then the FEM is directly applied to the linearized model. Solution of the resulting algebraic equations gives some optima1 values for the state and control variables. Using the new values the linearization procedure is repeated and the linear optimal control problem is solved again. The iterative solution terminates when convergence has occurred. Computationally, the new method is simple and convergence is achieved in a few iterations. Moreover, the new method is not very sensitive to initial guesses. tcurrently
FORMULATIONAND SOLUTIONOF THE NONLINEARCONTROLPROBLEM
where to and tr denote the initial and final time, respectively, and G and L are scalar functions. It is required to find the optimal control u*(t) which causes the general nonlinear system (1) with initial conditions (2) to follow an optima1 trajectory x*(f) and minimizes the performance index J. In some very simple cases, linear system dynamics (1) and a quadratic performance index (3), it is possible to obtain the optimal control law by numerically integrating a nonlinear matrix differential equation of the Ricatti type [2]. For nonlinear systems with quadratic or nonquadratic performance indices the resulting TPBVP is nonlinear, which is usually hard to solve. The numerical techniques that have been developed to solve the nonlinear differential equations mostly belong to indirect iterative methods. Nakamichi and Washizu [l] proposed a direct noniterative technique using the FEM for solving the linear control problem with a quadratic performance
a Ph.D. student at Lenhigh University, Be-
thlehem, PA 18015, U.S.A. II
78
Short Notes
index. A general formulation of the method is presented in the Appendix. Application of the FEM to a nonlinear optimal control problem (1) results in a system of nonlinear algebraic equations which can be solved numerically using a Newton-Raphson procedure [3]. In the present work, in order to avoid the problems associated with the solution of nonlinear algebraic equations, we employ an apparent linearization technique [4], that changes the nonlinear system into a linearized one. The basis of this approach is the adoption of an instantaneously linearized control and state model of the general form k = A(x, u, t) . x + B(x, II, t) * II,
(4)
where A and B are general (n x n) and (n x m) system matrices, respectively, found by linearization at a certain time subinterval. Note the selection of the coefficient matrices A and Ii3may not be unique and the following example illustrates this method. Consider a nonlinear system of the form i = g(x, u, t) + h(x, u, t).
(5)
This can be written in a linearized from as
(vi) Evaluate the W and Ho matrices, equations (A.10) and (A.1 1) and compute the optimal sequences of the control and state variables using equation (A.12) and (A.13). If the state variables are specified at rf then equations (A. 12) and (A. 13) have to be properly modified [1,51. (vii) Compare the ith and (i - 1)st trajectories for convergence: Ixck)(i) - x$?(i - 1)l < 6 Y
v=l,2,...,N lul’)(i)-ut”)(i-_)I<~
lim lg/x 1< co, lim Ih/u) x-0
<
03.
u-0
Thus, by providing instantaneous state and control values, the dependence of (g/x) and (h/u) on x and u is eliminated. We will refer to this form of system rearrangement as an apparent linearization. Subsequently we describe the proposed iterative procedure for solving a nonlinear optimal control problem with a quadratic performance index using the FEM. It is assumed that the final time is fixed and no restrictions are imposed on the control and state variables. The steps required to determine the optimal policy o*(t) are summarized as follows: (i) Rearrange the nonlinear equations, so that, if x and u are known, the system appears to be linear, equation (4). (ii) Divide the time interval [to, t,] into N subintervals and select the interpolation functions for the state variables, control variables and Lagrange multipliers, equations (A.3) and (A.4) in the Appendix. (iii) Assume some initial trajectories for the state and control variables. (iv) Use the state and control trajectories of Step (iii) to evaluate the system matrices A(x, u, t) and B(x, II, t) for each element. (v) Calculate the element submatrices IE,, D,, 1 and D6, for each element. Form the overall matrices K, ILand MOand the matrices P,, Poe, P, and PO, from the row vectors of P = K’lL and PO= M-l Ho matrices.
p=1,2,...,m;
(7)
E is a preselected termination constant. If equation (7) is satisfied the iterative procedure is terminated. (viii) If the solution has not converged, return to Step (iv) and reevaluate the A and B matrices using the results of the ith iteration. Application of the new method to the solution of nonlinear optimal control problems with a quadratic criterion is illustrated by the following example. A CONTINUOUS
where
k = 1,2 ,...,n
STIRRED-TANK
REACTOR
Consider a nonisothermal CSTR where a first-order irreversible exothermic reaction takes place. The flow of coolant through the coil inserted in the reactor is to control the performance of the reactor. The state equations that describe the dynamics of the system are given in Ref. [6]: 1, = -2(x, + 0.25) + (x2 + 0S)exp - (x1 + 0.25)~ 1, = 0.5 -x2 - (xZ + 0.5)ex with the initial conditions x(0) = (0.05, O)T.The states of the plant respresent the deviations of temperature and concentration from their corresponding steadystate values: x,(t) = T(t) - T,;
x2(r) = C(t) - C,;
(9)
u(t) is a normalized control variable representing the effect of coolant flow on the chemical reaction. The objective is to find the control policy u(t) that minimizes the quadratic criterion 0.78
J=
[x;(t) + x;(t) + O.lu*(t)]dt. (10) s0 Following the general methodology developed for nonlinear problems, equations (8) are first linearized using the apparent linearization approach:
i, =
Short Notes
79
1
State O.OSr
6
’
0.01 t
-0.01
*”
1%
-0.02
-m
-0.03
-
-a04
-
-0.05 -0.06
-
’
’
’
’
’
’
’
’
’
0.1
0.2
0.3
0.4
0.5 r/t,
0.6
0.7
0.6
0.9
1.0
I
:
:
I 0.4
I 0.5
I 0.6
I 07
I
I 09
I
0.1
0.2
0:3
.a
State
0.6
1.0
2
.
‘t $
0.024
I
0
35 Number
of
elements,
N
Fig. 3. Effect of the number of elements, N, on the calculated minimum value of J.
%
la
reaches an almost constant value of 0.0265 as the number of elements increases. This value agrees well with the minimum value of J,,( = 0.0266) obtained Q -0.07 L by Kirk [6] using the classical methods of SDM, Fig. 1. Effect of the number of elements, N, on the calculated VEM and QLM. optimal state trajectories (*not all points are plotted). To test the sensitivity of the new method to the initial guesses of xi, the optimal control problem was solved for different initial values of x, . In all runs the number of finite elements was equal to 10. The initial guess for x, varied from 10s3 to 10. In particular, the following initial values were assumed for x, :O.OOl, (11) 0.01, 0.05, 0.1, 0.5, 1, 2, 3 and 10. In all cases convergence was achieved in less than six iterations Numerical calculations were carried out for which indicates that the proposed FEM is not, at least for the particular problem, very sensitive to different numbers of elements, N, and different initial initial guesses. guesses of the state variable x,(x$“; i = 1,2, . . ., N). On the other hand, the number of iterations reThe iterative procedure was terminated when the convergence criterion for quired for convergence of the classical methods (7) was satisfied (SDM, VEM and QLM) largely depends on the c = 5 x 10e4. For all runs convergence was accominitial guesses [6]. It should be noted also that, for the plished in less than six iterations although the state reactor problem, the number of unknown variables, and control variables did not vary significantly after which must be initially specified in order to start the the second iteration. To find the effect of the number iterative solution, is equal to 1, 1, 2 and 4 for the of elements on the calculated optimal state and control policies, the total time interval was succes- FEM, SDM, VEM and QLM, respectively. Moreover, the VEM and QLM require initial values for the sively divided into 5, 10, 15, 20, 25 and 30 elements. The results of these runs are plotted ‘in Figs 1 and 2. costate variables which are usually difficult to guess. The FEM is computationally simple and requires We can see that there is no significant change in the the solution of (n + m)N linear algebraic equations at calculated policies for different values of N. However, as the number of elements increases an almost con- each iteration. These equations are derived directly tinuous trajectory is obtained. The FEM results are for the linearized model from the stationary conditions of the discretized functional (A.6). Note the in excellent agreement with those obtained by Kirk difference between the FEM and the QLM. In the [6] using the classical methods of steepest descent QLM, a sequence of linearized two-point boundary(SDM), variation of extremals (VEM) and quasilinearization (QLM). In Fig. 3, the minimum value of value problems is solved. The iterative procedure the performance index Jmi, is plotted as a function of requires the selection of initial state and costate the parameter N. It is interesting to note that Jmin trajectories x(t), p(t), re[t,,, t,] to linearize the governing state and costate differential equations. It is needless to say that a poor selection of the initial values might not produce convergent results. Although the new method has not been extensively tested yet, the results of this study show its good convergent characteristics. Therefore, the proposed FEM could provide an alternative tool for solving nonlinear optimal control problems. 0
‘?a
0.1
0.2
0.3
‘h
0.4
0.5 t/t.
0.6
0.7
0.6
0.9
1.0
Fig. 2. Effect of the number of elements, N, on the calculated optimal control trajectory.
Acknowledgements-We would like to thank the reviewers for their valuable comments. We are also grateful to A. Karagiannis for his help with some computer simulations.
Short Notes
80
7. D. Tabak and B. C. Kuo, Optimal Control by Mathematical Programming. Prentice-Hall, Englewood Cliffs,
NOMENCLATURE
N.J. (1971). 8. B. A. Finlayson, Nonlinear Analysis in Chemical Engineering. McGraw-Hill, New York (1980). 9. 0. Strang and G. J. Fix, An Analysis of the Finite Element Method. Prentice-Hall, Englewood Cliffs, N.J. (1973). 10. J. W. Daniel, The Approximate Minimization of Functionals. Prentice-Hall, Englewood Cliffs, N.J. (1971).
State matrix (n x m) Control matrix (2n x 2x) Element matrix (mN x mN) System matrix, equation (A.10) (m x m) Element matrix (mN x mN) System matrix, equation (A.10) (mN x n) System matrix, equation (A. 11) Performance index (nN x nN) System matrix, equation (A.14) (n x 2n) Element matrix (nN x n) System matrix, equation (A.14) (nN x mN) System matrix, equation (A.14) (n x m) Element matrix Dimension of control vector (2n x n) Derivative element matrix Dimension of the state vector (2n x n) Interpolation element matrix (nN x mN) System matrix (2n x mN) Element matrix (n x mN) Element matrix (nN x n) System matrix (2n x n) Element matrix (n x n) Element matrix (n x n) Weighting matrix, equation (A.2) (m x m) Weighting matrix, equation (A.2) (n x n) Weighting matrix, equation (A.2) time Final time (m x 1) Control vector (mN x 1) Overall discrete control vector (m x 1) Element control vector (n x 1) State vector (nN x 1) Overall discrete state vector (n x 1) Element state vector (n x 1) Vector of final values (2n x 1) Element discrete vector (n x 1) Final element vector (n x 1) Vector of initial values
A = (n x n)
B= D, = E= E, = H= W, = J = K= DQ,= 1)6,= IL= IL,= m = M, = n = N, = P= P, = P, = P, = P, = P, = Q= W= S, = t= rr = u= U= u, = x= X= x, = xr= X, = x, = x, =
APPENDIX The Finite-element
i(t) = A(t)x(t) + B(t)u(t)
J = jxT(t,)§,x(t,) +;
t [xT(t)B(t)x(t) s ro + u(t)TR(t)u(t)] dr, (A.2) where Sr, Q(r) and R(t) are weighting matrices. Following the approach of Nakamichi and Washizu [I], the time interval [rO,rr] is divided into N time subintervals (elements), Fig. Al, and interpolation functions for all variables are defined in each element. Thus for an element “e” we obtain:
6.
I,(r) = L
(A.4)
u,(t) = u,;
XT(Z) x(*) e = [x(‘)(t), e e9(t)
. .,x’“‘(t)] c1
XT c = [xj”,, xp, xi”,, xp, . . . , xp,, xpy and
N:(t) = ti - t --.
REFERENCES
4. 5.
(A.3)
where
1 = (n x 1) Vector of Lagrange multipliers
3.
x,(t) = NT(t)X,, jr.(t) = &4:x, and
c = Preselected termination constant
2.
64.1)
with the initial conditions (2) and a quadratic performance index
Greek letters
1.
Approach of Nakamichi and Washizu [I]
Consider a linear differential system
t-t,-, At
At ’
Nakamichi and K. Washizu, .Z.Franklin Inst. 306,309 (1978). A. P. Sage, Optimum System Control. Prentice-Hall, Ennlewood Cliffs. NJ. (1968). J. hakamichi and K. ‘Washizu, Znt J. nurner. Meth. Engng 12, 1559 (1978). J. D. Pearson, J. Electron. Control 13, 453 (1962). A. Georgiou, Application of finite elements to the solution of optimal control problems. Diploma Thesis, Univ. of Thessaloniki, Thessaloniki, Greece (1982). D. E. Kirk, Optimal Control Theory. Prentice-Hall, Englewood Cliffs, N.J. (1970).
0
J.
1
0
2
I
‘0 lkl
*o
4
‘2
,
.
‘3
.
t--i-,
’
At’dl
ti - t
0
At ’
t--t,_, At
Each state variable is approximated by a linear function of time
= fxf’-lx!“+f’-xj~,, At ’ At
xqf) e
i
i-l
3
-ti - t
.
N-l
i+l
.k_,
.
5-r
f I
*(k) i-1
.
.
N
.,p,
Ir+,
tN-l
!v
*MI
*lkl
i
N
(P’
u/
$1
Fig. Al FEM discretization.
k:1,2,...,n p:l,2,...,m
(A.3
Short Notes where superscript (k) denotes the k th state variable (k: 1,2,..., n). Note that MT is the time derivative of NT(t). Then we define an augmented performance index Jo in terms of equations (A.1) and (A.2). After substitution of equations (A.3) and (A.4) into the new index Jo, we obtain the following discretized form of Jo:
81
equations of the following form: wu = WC&, w = lPT,s,p,+
(
w,= Jo = fx$§,x,
+ ; ;(x:D,x, r=,
+ n:@$x,
+ lu,),
(A.6)
where $
li E,=
R(t) dt, s &--I 4 L, = B(t) dt, I 11-I
D, =
N,(t)Q(t)Nf(t)
d&
(A.?
I 4-I K, =
’ [A(t)N:(t) - MT] dt (A.8) I h- I
and XN= [XC’,x!$, . . , , xl;l’]. TO establish the necessary conditions for a minimum, we require that the first variation in Jo vanishes for arbitrary variation in SU. This results in a system of linear algebraic
i P;D,P,+E c-1
-PT,§,P,-;
(
+ ufE,u,)
(A.9)
P,D,P, r=l
,
(A.10)
,
(A.ll)
> >
u = w-‘W,x,
(A.12)
X=(PW-‘w,+P,)x,
(A.13)
The matrices P,, Plk, PN and P, are formed from the row vectors of the P = K-IL and P, = K-‘06, matrices. Where K, ILand W, are the assembled matrices obtained by the summation of all element submatrices K, and 1: KX = LU + I&&). (A.14) For more details the reader should refer to the original work of Nakamichi and Washizu [l] or alternatively to Georgiou [5]. A general computer program has been developed capable of solving the linear optimal control problem with a quadratic criterion using the FEM. Computation of various element submatrices and overall matrices is straightforward. The dimensions bf the different matrices are given in the Nomenclature.
0098-1354/87
Copyright 0
1987
$3.00 +o.oll
Pergamon Journals Ltd
SHORT NOTES FINITE-ELEMENT SOLUTION OF NONLINEAR OPTIMAL CONTROL PROBLEMS WITH A QUADRATIC PERFORMANCE INDEX C. KIPARI~SDFSand A. GEORCIOU~ Department of Chemical Engineering, University of Thessaloniki, Thessaloniki, Greece (Received 18 January 1984;/inal revision received 31 Juiy 1986; received for publication 22 August 1986)
Abstract-A new iterative method is developed
for the solution of nonlinear optimal control problems with a quadratic performance index. The new method employs a local linearization technique together with the finite-element method to obtain the optimal state and control sequences. Preliminary results show that the proposed method is simple, accurate and convergence is achieved in a few iterations. The new technique is not very sensitive to initial guesses and it is shown to perform at least as well as the classical techniques (i.e. steepest descent, quasilinearization).
INTRODUCTION The numerical methods that have been developed to solve the optimal control problem can be broadly classified into direct and indirect methods, Furthermore, they can be of iterative or noniterative type. Nakamichi and Washizu [l] developed a direct noniterative approach for solving the linear optimal
Consider any controlled process operating over the fixed interval, t e[tO, t,], which can be described by a system of nonlinear differential equations ir = f(x, II, t),
(1)
where x is the state vector of dimension n and u is the control vector of dimension m. Given the initial condition vector for x,
control problem with a quadratic performance index. The main advantage of their method is that the necessary conditions for the augmented performance index to be stationary can be expressed in the form of a system of algebraic equations. As a result, the optimal state and control trajectories are obtained by solving a system of linear algebraic equations. However, application of the finite-element method (FEM) to nonlinear problems results in a system of nonlinear
x(4l) = %9
(2)
the solution of equation (1) is uniquely defined, once the control vector u is chosen. Consider the performance index J, If L(x, UJ) dt J = G(x(t,), 4) + (3) s @I
algebraic equations which are hard to solve. This seems to mitigate against the effectiveness and simplicity of the technique for nonlinear problems. To circumvent this problem, we propose a new iterative approach to the solution of nonlinear control problems with a quadratic performance index using the FEM. According to the new method, a local linearization technique is used to obtain a linear model, then the FEM is directly applied to the linearized model. Solution of the resulting algebraic equations gives some optima1 values for the state and control variables. Using the new values the linearization procedure is repeated and the linear optimal control problem is solved again. The iterative solution terminates when convergence has occurred. Computationally, the new method is simple and convergence is achieved in a few iterations. Moreover, the new method is not very sensitive to initial guesses. tcurrently
FORMULATIONAND SOLUTIONOF THE NONLINEARCONTROLPROBLEM
where to and tr denote the initial and final time, respectively, and G and L are scalar functions. It is required to find the optimal control u*(t) which causes the general nonlinear system (1) with initial conditions (2) to follow an optima1 trajectory x*(f) and minimizes the performance index J. In some very simple cases, linear system dynamics (1) and a quadratic performance index (3), it is possible to obtain the optimal control law by numerically integrating a nonlinear matrix differential equation of the Ricatti type [2]. For nonlinear systems with quadratic or nonquadratic performance indices the resulting TPBVP is nonlinear, which is usually hard to solve. The numerical techniques that have been developed to solve the nonlinear differential equations mostly belong to indirect iterative methods. Nakamichi and Washizu [l] proposed a direct noniterative technique using the FEM for solving the linear control problem with a quadratic performance
a Ph.D. student at Lenhigh University, Be-
thlehem, PA 18015, U.S.A. II
78
Short Notes
index. A general formulation of the method is presented in the Appendix. Application of the FEM to a nonlinear optimal control problem (1) results in a system of nonlinear algebraic equations which can be solved numerically using a Newton-Raphson procedure [3]. In the present work, in order to avoid the problems associated with the solution of nonlinear algebraic equations, we employ an apparent linearization technique [4], that changes the nonlinear system into a linearized one. The basis of this approach is the adoption of an instantaneously linearized control and state model of the general form k = A(x, u, t) . x + B(x, II, t) * II,
(4)
where A and B are general (n x n) and (n x m) system matrices, respectively, found by linearization at a certain time subinterval. Note the selection of the coefficient matrices A and Ii3may not be unique and the following example illustrates this method. Consider a nonlinear system of the form i = g(x, u, t) + h(x, u, t).
(5)
This can be written in a linearized from as
(vi) Evaluate the W and Ho matrices, equations (A.10) and (A.1 1) and compute the optimal sequences of the control and state variables using equation (A.12) and (A.13). If the state variables are specified at rf then equations (A. 12) and (A. 13) have to be properly modified [1,51. (vii) Compare the ith and (i - 1)st trajectories for convergence: Ixck)(i) - x$?(i - 1)l < 6 Y
v=l,2,...,N lul’)(i)-ut”)(i-_)I<~
lim lg/x 1< co, lim Ih/u) x-0
<
03.
u-0
Thus, by providing instantaneous state and control values, the dependence of (g/x) and (h/u) on x and u is eliminated. We will refer to this form of system rearrangement as an apparent linearization. Subsequently we describe the proposed iterative procedure for solving a nonlinear optimal control problem with a quadratic performance index using the FEM. It is assumed that the final time is fixed and no restrictions are imposed on the control and state variables. The steps required to determine the optimal policy o*(t) are summarized as follows: (i) Rearrange the nonlinear equations, so that, if x and u are known, the system appears to be linear, equation (4). (ii) Divide the time interval [to, t,] into N subintervals and select the interpolation functions for the state variables, control variables and Lagrange multipliers, equations (A.3) and (A.4) in the Appendix. (iii) Assume some initial trajectories for the state and control variables. (iv) Use the state and control trajectories of Step (iii) to evaluate the system matrices A(x, u, t) and B(x, II, t) for each element. (v) Calculate the element submatrices IE,, D,, 1 and D6, for each element. Form the overall matrices K, ILand MOand the matrices P,, Poe, P, and PO, from the row vectors of P = K’lL and PO= M-l Ho matrices.
p=1,2,...,m;
(7)
E is a preselected termination constant. If equation (7) is satisfied the iterative procedure is terminated. (viii) If the solution has not converged, return to Step (iv) and reevaluate the A and B matrices using the results of the ith iteration. Application of the new method to the solution of nonlinear optimal control problems with a quadratic criterion is illustrated by the following example. A CONTINUOUS
where
k = 1,2 ,...,n
STIRRED-TANK
REACTOR
Consider a nonisothermal CSTR where a first-order irreversible exothermic reaction takes place. The flow of coolant through the coil inserted in the reactor is to control the performance of the reactor. The state equations that describe the dynamics of the system are given in Ref. [6]: 1, = -2(x, + 0.25) + (x2 + 0S)exp - (x1 + 0.25)~ 1, = 0.5 -x2 - (xZ + 0.5)ex with the initial conditions x(0) = (0.05, O)T.The states of the plant respresent the deviations of temperature and concentration from their corresponding steadystate values: x,(t) = T(t) - T,;
x2(r) = C(t) - C,;
(9)
u(t) is a normalized control variable representing the effect of coolant flow on the chemical reaction. The objective is to find the control policy u(t) that minimizes the quadratic criterion 0.78
J=
[x;(t) + x;(t) + O.lu*(t)]dt. (10) s0 Following the general methodology developed for nonlinear problems, equations (8) are first linearized using the apparent linearization approach:
i, =
Short Notes
79
1
State O.OSr
6
’
0.01 t
-0.01
*”
1%
-0.02
-m
-0.03
-
-a04
-
-0.05 -0.06
-
’
’
’
’
’
’
’
’
’
0.1
0.2
0.3
0.4
0.5 r/t,
0.6
0.7
0.6
0.9
1.0
I
:
:
I 0.4
I 0.5
I 0.6
I 07
I
I 09
I
0.1
0.2
0:3
.a
State
0.6
1.0
2
.
‘t $
0.024
I
0
35 Number
of
elements,
N
Fig. 3. Effect of the number of elements, N, on the calculated minimum value of J.
%
la
reaches an almost constant value of 0.0265 as the number of elements increases. This value agrees well with the minimum value of J,,( = 0.0266) obtained Q -0.07 L by Kirk [6] using the classical methods of SDM, Fig. 1. Effect of the number of elements, N, on the calculated VEM and QLM. optimal state trajectories (*not all points are plotted). To test the sensitivity of the new method to the initial guesses of xi, the optimal control problem was solved for different initial values of x, . In all runs the number of finite elements was equal to 10. The initial guess for x, varied from 10s3 to 10. In particular, the following initial values were assumed for x, :O.OOl, (11) 0.01, 0.05, 0.1, 0.5, 1, 2, 3 and 10. In all cases convergence was achieved in less than six iterations Numerical calculations were carried out for which indicates that the proposed FEM is not, at least for the particular problem, very sensitive to different numbers of elements, N, and different initial initial guesses. guesses of the state variable x,(x$“; i = 1,2, . . ., N). On the other hand, the number of iterations reThe iterative procedure was terminated when the convergence criterion for quired for convergence of the classical methods (7) was satisfied (SDM, VEM and QLM) largely depends on the c = 5 x 10e4. For all runs convergence was accominitial guesses [6]. It should be noted also that, for the plished in less than six iterations although the state reactor problem, the number of unknown variables, and control variables did not vary significantly after which must be initially specified in order to start the the second iteration. To find the effect of the number iterative solution, is equal to 1, 1, 2 and 4 for the of elements on the calculated optimal state and control policies, the total time interval was succes- FEM, SDM, VEM and QLM, respectively. Moreover, the VEM and QLM require initial values for the sively divided into 5, 10, 15, 20, 25 and 30 elements. The results of these runs are plotted ‘in Figs 1 and 2. costate variables which are usually difficult to guess. The FEM is computationally simple and requires We can see that there is no significant change in the the solution of (n + m)N linear algebraic equations at calculated policies for different values of N. However, as the number of elements increases an almost con- each iteration. These equations are derived directly tinuous trajectory is obtained. The FEM results are for the linearized model from the stationary conditions of the discretized functional (A.6). Note the in excellent agreement with those obtained by Kirk difference between the FEM and the QLM. In the [6] using the classical methods of steepest descent QLM, a sequence of linearized two-point boundary(SDM), variation of extremals (VEM) and quasilinearization (QLM). In Fig. 3, the minimum value of value problems is solved. The iterative procedure the performance index Jmi, is plotted as a function of requires the selection of initial state and costate the parameter N. It is interesting to note that Jmin trajectories x(t), p(t), re[t,,, t,] to linearize the governing state and costate differential equations. It is needless to say that a poor selection of the initial values might not produce convergent results. Although the new method has not been extensively tested yet, the results of this study show its good convergent characteristics. Therefore, the proposed FEM could provide an alternative tool for solving nonlinear optimal control problems. 0
‘?a
0.1
0.2
0.3
‘h
0.4
0.5 t/t.
0.6
0.7
0.6
0.9
1.0
Fig. 2. Effect of the number of elements, N, on the calculated optimal control trajectory.
Acknowledgements-We would like to thank the reviewers for their valuable comments. We are also grateful to A. Karagiannis for his help with some computer simulations.
Short Notes
80
7. D. Tabak and B. C. Kuo, Optimal Control by Mathematical Programming. Prentice-Hall, Englewood Cliffs,
NOMENCLATURE
N.J. (1971). 8. B. A. Finlayson, Nonlinear Analysis in Chemical Engineering. McGraw-Hill, New York (1980). 9. 0. Strang and G. J. Fix, An Analysis of the Finite Element Method. Prentice-Hall, Englewood Cliffs, N.J. (1973). 10. J. W. Daniel, The Approximate Minimization of Functionals. Prentice-Hall, Englewood Cliffs, N.J. (1971).
State matrix (n x m) Control matrix (2n x 2x) Element matrix (mN x mN) System matrix, equation (A.10) (m x m) Element matrix (mN x mN) System matrix, equation (A.10) (mN x n) System matrix, equation (A. 11) Performance index (nN x nN) System matrix, equation (A.14) (n x 2n) Element matrix (nN x n) System matrix, equation (A.14) (nN x mN) System matrix, equation (A.14) (n x m) Element matrix Dimension of control vector (2n x n) Derivative element matrix Dimension of the state vector (2n x n) Interpolation element matrix (nN x mN) System matrix (2n x mN) Element matrix (n x mN) Element matrix (nN x n) System matrix (2n x n) Element matrix (n x n) Element matrix (n x n) Weighting matrix, equation (A.2) (m x m) Weighting matrix, equation (A.2) (n x n) Weighting matrix, equation (A.2) time Final time (m x 1) Control vector (mN x 1) Overall discrete control vector (m x 1) Element control vector (n x 1) State vector (nN x 1) Overall discrete state vector (n x 1) Element state vector (n x 1) Vector of final values (2n x 1) Element discrete vector (n x 1) Final element vector (n x 1) Vector of initial values
A = (n x n)
B= D, = E= E, = H= W, = J = K= DQ,= 1)6,= IL= IL,= m = M, = n = N, = P= P, = P, = P, = P, = P, = Q= W= S, = t= rr = u= U= u, = x= X= x, = xr= X, = x, = x, =
APPENDIX The Finite-element
i(t) = A(t)x(t) + B(t)u(t)
J = jxT(t,)§,x(t,) +;
t [xT(t)B(t)x(t) s ro + u(t)TR(t)u(t)] dr, (A.2) where Sr, Q(r) and R(t) are weighting matrices. Following the approach of Nakamichi and Washizu [I], the time interval [rO,rr] is divided into N time subintervals (elements), Fig. Al, and interpolation functions for all variables are defined in each element. Thus for an element “e” we obtain:
6.
I,(r) = L
(A.4)
u,(t) = u,;
XT(Z) x(*) e = [x(‘)(t), e e9(t)
. .,x’“‘(t)] c1
XT c = [xj”,, xp, xi”,, xp, . . . , xp,, xpy and
N:(t) = ti - t --.
REFERENCES
4. 5.
(A.3)
where
1 = (n x 1) Vector of Lagrange multipliers
3.
x,(t) = NT(t)X,, jr.(t) = &4:x, and
c = Preselected termination constant
2.
64.1)
with the initial conditions (2) and a quadratic performance index
Greek letters
1.
Approach of Nakamichi and Washizu [I]
Consider a linear differential system
t-t,-, At
At ’
Nakamichi and K. Washizu, .Z.Franklin Inst. 306,309 (1978). A. P. Sage, Optimum System Control. Prentice-Hall, Ennlewood Cliffs. NJ. (1968). J. hakamichi and K. ‘Washizu, Znt J. nurner. Meth. Engng 12, 1559 (1978). J. D. Pearson, J. Electron. Control 13, 453 (1962). A. Georgiou, Application of finite elements to the solution of optimal control problems. Diploma Thesis, Univ. of Thessaloniki, Thessaloniki, Greece (1982). D. E. Kirk, Optimal Control Theory. Prentice-Hall, Englewood Cliffs, N.J. (1970).
0
J.
1
0
2
I
‘0 lkl
*o
4
‘2
,
.
‘3
.
t--i-,
’
At’dl
ti - t
0
At ’
t--t,_, At
Each state variable is approximated by a linear function of time
= fxf’-lx!“+f’-xj~,, At ’ At
xqf) e
i
i-l
3
-ti - t
.
N-l
i+l
.k_,
.
5-r
f I
*(k) i-1
.
.
N
.,p,
Ir+,
tN-l
!v
*MI
*lkl
i
N
(P’
u/
$1
Fig. Al FEM discretization.
k:1,2,...,n p:l,2,...,m
(A.3
Short Notes where superscript (k) denotes the k th state variable (k: 1,2,..., n). Note that MT is the time derivative of NT(t). Then we define an augmented performance index Jo in terms of equations (A.1) and (A.2). After substitution of equations (A.3) and (A.4) into the new index Jo, we obtain the following discretized form of Jo:
81
equations of the following form: wu = WC&, w = lPT,s,p,+
(
w,= Jo = fx$§,x,
+ ; ;(x:D,x, r=,
+ n:@$x,
+ lu,),
(A.6)
where $
li E,=
R(t) dt, s &--I 4 L, = B(t) dt, I 11-I
D, =
N,(t)Q(t)Nf(t)
d&
(A.?
I 4-I K, =
’ [A(t)N:(t) - MT] dt (A.8) I h- I
and XN= [XC’,x!$, . . , , xl;l’]. TO establish the necessary conditions for a minimum, we require that the first variation in Jo vanishes for arbitrary variation in SU. This results in a system of linear algebraic
i P;D,P,+E c-1
-PT,§,P,-;
(
+ ufE,u,)
(A.9)
P,D,P, r=l
,
(A.10)
,
(A.ll)
> >
u = w-‘W,x,
(A.12)
X=(PW-‘w,+P,)x,
(A.13)
The matrices P,, Plk, PN and P, are formed from the row vectors of the P = K-IL and P, = K-‘06, matrices. Where K, ILand W, are the assembled matrices obtained by the summation of all element submatrices K, and 1: KX = LU + I&&). (A.14) For more details the reader should refer to the original work of Nakamichi and Washizu [l] or alternatively to Georgiou [5]. A general computer program has been developed capable of solving the linear optimal control problem with a quadratic criterion using the FEM. Computation of various element submatrices and overall matrices is straightforward. The dimensions bf the different matrices are given in the Nomenclature.