Differential Dynamic Programming as a Tool for Dynamic Optimization in a Travelling-Grate Pelletizing Plant

Copyright (f) IFAC Auto m ation in Mining, Mineral and Metal Processing, Helsinki, Finland, 1983

DIFFERENTIAL DYNAMIC PROGRAMMING AS A TOOL FOR DYNAMIC OPTIMIZATION IN A TRAVELLING-GRATE PELLETIZING PLANT C. Hillberg Control Engineering Laboratory, Chalmers University of Technology, 5-41296 Gothenburg, Sweden

Abstract . Induration in a travelling- grate pelletizing plant is a nonlinear distributed parameter process . A high-order dynamic model of the induration zone has been simplified into a model of low order using the method of subdomains. The Differential Dynamic Programming (DDP) method is used for solving an optimal control problem containing the simplified model. The problem of minimizing the inlet gas energy necessary for the induration is considered. The purpose is to examine the feasibility of the DDP method as a tool for dynamic optimization in this kind of process. The approach of the method of subdomains used in the simplification is briefly that the distributed variables are approximated by truncated power series, in the space variable, with time-dependent coefficients. The approximations are inserted into the partial differential equations, forming residuals which are integrated over certain subdomains. The result of this procedure is a non-linear model of low order with time as the independent variable. The optimal control problem treated here has the inlet gas temperature as the control variable . The greater loss of energy at high temperatures is covered by using a quadratic func tion of the inlet gas temperature in the cost functional. Furthermore, the derivative of pellet temperature with respect to time is limited by a specified value in order to avoid heat shock. The influence of the energy obtained through the exothermal reaction of magnetite into hematite is considered. Variations in the model parameters bed height and pressure drop are performed to determine a defined efficiency of the energy of the inlet gas. Results of the optimization would be to use high beds and as high as possible a gas temperature. The results also show that the method of sub domains in combination with DDP forms a powerful tool for the discovery and verification of control strategies in processes of this type . Keywords. Distributed parameter system; non-linear equations; pelletizing plant; induration zone; optimal control . INTRODUCTION The induration of ~ron ore ~n a travellinggrate pelletizing plant is an example of a h i ghly non-linear process. The minimization of energy consumed per ton of product is important, as the process is highly energyconsuming . The mathematical models of the process available in the literature are often very complex . Hence, optimization from a theoretical point of view will be difficult and time consuming, if it is possible at all . The use of simplified models will be advantageous according to Voskamp and Vissers (1973) and Drugge (1975, 1981) among others. Optimization is often performed by using offline simulations and tests in pot furnaces as is done in Batterham et. a1. (1977). This paper presents the Differential Dynamic Programming (DDP) method, as a tool for dynamic optimization in non-linear processes

of this type. Some applications of DDP in a different field can be found in Hillberg and Jarmark (1983). The simplified model is obtained using the method of subdomains, which has already been successfully applied to the updraught drying zone by Breitholtz and Hillberg (1981). The main point of this paper is not the optimal result itself, but the fact that the DDP method in conjunction with the method of subdomains, provides a means of analyzing and optimizing this problem. tl0DELLING The notations used ~n this paragraph are explained i n Table 1. Assumptions: The stationary state of the process is represented in the model by a travelling

377

378

C. Hillberg segment as in Fig. 1, which moves at a constant speed through the machine. Consider a segment with the base area equal to one square meter. The bed dimensions are assumed to be constant during the process. Since we have assumed that the gas flow is downdraught, neither the hearth layer nor the grate bars are included in the calculations. It is quite possible to include the effect of internal temperature distributions of the pellets, as in the paper by Breitholtz and Hillberg (1980a). In this simple model, however, only the bulk temperature of the pellets at a specific level of the bed, z , is considered.

The parameters kO and kl can be found through measurements of the material used in the plant (Voskamp and Brasz, 1975). The heat released from the oxidation will then be q

ox

The effect of the nonisothermal gas flow through the bed is accounted for by using the extended Ergun relation (Szekely and Carr, 1968; Breitholtz and Hillberg, 1980b). 2 2 G / (2E:) log( p (z h

! (AG~ (T (t,z»

+

o

Let Ts(t,z) and Tg(t,z) be the temperatures of the segment and gas, and let qox be the heat released from the oxidation of magnetite. The heat balance for the segment modelled, as in Breitholtz and Hillberg (1980a), becomes

!

+ (1)

and the heat balance for the gas aT (t, z) g

g

az

= - a A(T (t,z) - T (t,z» g

s

p

p(O)

2

+ BG ) dz

(t, z»

(T

dp

(7)

g ' g

(2)

1. 75(1 - E)

B

q, d

E3

p

p (T ) = 348./(T g

g

3

g

+ 273.) kg/m , T g (8)

~

Gc

g

+

A

The speed of the travelling segment is constant. Hence, the x axis corresponds to a time axis.

s

g

h»

p(h)

The water content of the segment on entering the zone is assumed to be negligible.

g

0) / p (z g

g

Only convective heat transfer between gas and pellets is considered.

aA(T (t,z) - T (t,z»

(6)

(t, z)

g

(T ) g

(20.4 + 0.029 T ) g

. 10- 6 Pa s, T

~n °c

g

(9)

If we use the above approximations for the gas density p and the gas viscosity ~ g' and assume tha~ we can express the gas temperature as N

g

The heat transfer coefficient, a , can be changed into an overall heat transfer coefficient, since it is often used for adjustment in the model to fit the practical conditions of a plant.

T (t, z) g

The exothermal reaction of magnetite into hematite

The above model will now be simplified by the method of subdomains, (Breitholtz and Hillberg, 19n1; Finlayson, 1972). We thus assume solutions in the following forms:

is modelled ac cording to Voskamp and Brasz (1975). The following equation is assumed to hold.

~

i=O

(10)

then we can calculate a flow G, at each instant t , for a given pressure drop ~p

N

s

a. (z/h)

T (t, z) g

g b. (z/h)

l:

c

l:

i=O

i

~

i=O N

i

(11)

~

i=O N

C (t ,z)

s

l:

T (t, z)

(4)

where C(t,z) is the proportion of magnetite by weight in a pellet. The reaction rate KB is approximated by

b.(t)(z/h)i

l:

c. (z/h)

i

~

and expand the reaction rate power series

(12)

(13) in the

Differential Dynamic Programming NB f. Ti KB(T) = E 1 i=O

(14)

We form the residuals Vp aT s s Vp c - aA(T -T ) - (- t> H) s s at g s M KBC

R s

(lS)

379

more non-linear functional can be used as well. Furthermore, we require that the pellet temperature reach a specified temperature at the end of the zone. That is, we want the bulktemperature Tsb(t) to remain at the level Tf during the interval tl < t < t f Tsb is defined as

aT

~- aA(T

R

- Qc

R c

dC a t + KBC

g

g az

g

- T ) s

N

(16)

s a./(i + 1) T (t,z)dz=E s i=O 1

( 24)

(17)

Hence, we will add the cost

and then integrate the residuals Rs' Rg and Rc over the subdomains. The result is two n on-linear differential equations with time as the independent variable, and one algebraic equation.

(2S)

to the cost functional.

i

0,1, .. ,N

b.

i

1

(18)

s

1,2, .. , N g

(19)

To avoid heat shock the heating rate of the pellets must not exceed a specified value Yp °C/s . This value is, of course, dependent on the material used. For example the value 200 °C/min applies to LKAB magnetite concentrate (Drugge, 1975). Hence, the constraint g ( u, T s; ) t =

m~x

{ aTsat (t, z) } _ Y p _ < 0 , 0
i

=

(20)

Q,l, .. , N c

Furthermore, let the gas temperature at z = 0 be the control variable u(t)

Eq. (19) can be used to eliminate bi in Eq. (18). Hence, the system is modelled by Ns + Nc + 2 ordinary differential equations. OPTI1~

For the sake of simplicity let the vector x be concatenated by the vectors a and c and let xo be the value of x at t = 0 , (see Eqs. 11,12,13). The optimal control problem is then formulated by the following. !linimize the cost functional t

CONTROL PROBLEM

g g

(22)

2

o

Since the pelletizing process consumes much energy, we have tried to minimize the inlet gas energy to the zone. The inlet gas energy c an be expressed as

Gc T (t,O)dt

f

ql f O.S(Gu) dt +

V(u)

E. 1n

(27)

g(u,Ts;t) = 0 (21)

THE

must be satisfied. The limit value u for u , can be calculated at each instant from the constraint condition

( 28)

subject to x. = K.(x.,u) 1

In order to compensate to some extent for the greater loss of energy at high temperatures, we have used the cost functional t

1

i,j = 1,2, ... ,N g (x.,u;t)

f

O.S f

J

2

(Gu) dt

(23)

o Another alternative could be that Cl is the cost of the additional energy used. There is no requirement by the optimizing program that a quadratic cost functional must be used. A

x. (0)

J

~

1

s

0 , j

+ Nc + 2

1,2, ... ,N

(29)

S

+ 1

(30)

where K. stands for the non-linear functions that resGlt from the integration of the residuals (Eqs. 18 and 20).

380

C. Hillberg !ffiTHOD OF SOLUTION

The problem stated above was now solved by the DDP method, complemented by a convergence control technique, CCP, which is described in Jacobson and Hayne (1970) and Hillberg and Jarmark (1983). The DDP method, which will give us an open loop optimal control, works iteratively, starting from a guessed control history li(t) which is successively improved. Using Ns ~ Nc ~ N ~ 3 , the process considered here is m03elled by eight differential equations. Furthermore, the DDP method requires the adjoint equations

vx

(31)

where H is the Hamiltonian and is a multiplier defined as

~(x,u;t)

(32)

if g(x,u;t) is active, else zero. Hence, 2 . 8 ~ 16 differential equations have to be formulated and integrated. Using a secondorder DDP method, which usually gives a better convergence when approaching the optimum, the number of differential equations to be implemented would be 52. A first-order method will in most applications give a satisfactory optimal solution. In this study a first-order DDP method was used. The computational procedure followed for DDP is as follows: First we chose a reasonable control history li(t), which applied to Eq. (29) and the cost functional will give us a ~ominal trajectory x(t) and a nominal cost V . The boundary conditions of Eq. (31) is now known and we integrate Eq. (31) backwards in time, ,,,hile we determine a new control u*(t) by mini~ization of the Ha~iltonian

u*~argminH(~,u,V ;t), g(~,u;t) u

x

< 0

(33)

Also we integrate aCt)

~

- * ,Vx;t) + H(x,u,Vx;t),a(t ) -H(x,u f

~

0, (34)

where aCt) is the predicted cost change obtained when u*(t) is applied instead of ti(t) . Eq. (29) is now integrated forward in time using u*(t) and a new cost value v*is obtained. A comparison of the real cost change V* - V and the predicted one, a(O) , will give us information about the convergence status. If the status is satisfactory let ti ~ u* and x ~ x* and repeat the iteration procedure until a(O) is less than a specified value. The DDP ~ethod is more thoroughly explained in the references cited above. The fact that the flow rate G is dependent on T (t,z) was solved in the following way. When ~he algorithm was started a flow distribution suiting the nominal control ti(t) was calculated by Eq. (7). After the optimum was reached, we calculated a new distribution and started the optimization procedure again. This method worked well and the number of flow rate

calculations needed here is two or three if the li(t) is far away from the optimal control. The model and the optimization program could be extended to cover a bigger part of the machine, but in the present study the final time tf was set at 7.5 min. The pressure drop over the bed is a parameter to the program and several values were tried, for instance: a pressure drop of 2500 Pa for the first 4.5 minutes and 4200 Pa for the remaining time to 7.5 minutes, which applies to the Halmberget plant, Sweden (see Drugge, 1975). RESULTS AND DISCUSSION To determine the degree of the polynomials to be used, we can start with a low degree and then raise it until the model gives a satisfactory description of the system. Further information on error indication can be found in Breitholtz (1982). Test simulations were done and the result using cubic profiles is depicted in Figs. 2, 3 and 4. In these first figures an overly high reaction rate for the oxidation and a reasonable control were used to show the performance of the model in an extreme situation. The top temperature in Fig. 2 shows an overshoot when t is about 6 min. The reason for this is found in Fig. 4, where we can see that the concentration of magnetite is almost zero at z ~ 0 after six minutes, and therefore no heat contribution is obtained from the oxidation process. Hence, the top temperature at this point will descend. Fig. 4 shows us that a polynomial of the third degree, modelling the concentration of magnetite, can not model C(t,z) ~ 0 for O~z,;:zc and C(t,z)iO for zc~z.:::h properly. The result when a polynomial of the fifth degree is used can be found in Fig. 5. Effects on the pellet temperature from this change in degree were small. The increase in computing time when the order of the concentration polynomial was raised from three to five was about 20 %. This indicates that a simulation using all polynomials of the fifth order is not out of the question with respect to computing time. The optimization program, in which the simulation part can be extracted, was written in FORTRAN and run on a PDP 11-23 computer. The computing time for the simulation in Fig. 2, using a step length of five seconds, was 1.5 min. including output printing. This figure can be lowered if a program designed only for simulation is used. It is the opinion of the author that the approximative model reproduces the process state in a satisfactory way and that this kind of models are well suited for implementation in a simulator used as an operator~s guide. The input temperature of the gas was now optimized with the first order DDP method for

Differential Dynamic Programming several pressure drops. The results are found in Figs. 6 - 9 and in Table 2. The following properties have to be defined. Let the outlet energy, Eout ' be defined in a way similar to that of the inlet energy t

E

out

f

f c GT (t,h)dt

o

g

(35)

g

and ~et r be a measure of the efficiency of the ~nlet energy used for the oxidation of magnetite. The efficiency r is formulated as (-6H)Vp r

=

s

h

0.35

(36)

where CTot denotes the ratio of magnetite left in the bed at t = t . f 1 h

h

f C(tf,z)dz

(37)

o

The optimal control policy is to let the inlet gas temperature be as high as possible. In this example the gas temperature was decreased when the pellet bulk temperature of 1000 °c was reached, to a magnitude that was enough to keep the bulk temperature at constant level, as was formulated in the cost functional. The limit value u for u, was as we have said, determined by the maximal heating rate allowed. The maximal temperature derivative of the pellets with respect to time occurs at z = 0 for t < 5 , and at o < z < h when the process proceeds; this is dependent of given process parameters (see Fig. 8). T~e n~mber

of iterations needed for the optiprogram was about five, see Fig. 9, and each iteration reauired about three minutes' computing tim~. The total number of iterations was dependent on how far away from optimum the initial quessed control history was, and on the accuracy demands made on the soluti on. m~ z at~on

In Table 1 we can see that the efficiency r increases for falling pressure drop when the be d height is constant. It also increases with increasing bed height, at a c onstant pressure drop. Hence, the control policy would be.to use high beds and a low pressure drop, wh~ch would mean a small gas flow. This also leads t o that the loss of energy in the outlet gas is decreased. This optimization made here is perhaps not suitable for on-line implementation or as an on-line operator~s guide if further simplifications in model and program can not be made to decrease the computing time. It could, however, still be a powerful tool for developing and verifying control strategies or for finding the optimal way of operating this kind of process.

381

Some of the model equations used, Eqs. (1- ~ may have some minor defiencies, but the main aim is to present and demonstrate the DDP and the subdomain methods as tools for the dynamic optimization of a process like this. Also by performing this optimization we can learn much about the process. A suggestion for further work would be to improve the accuracy of the model, by for example, having a more accurate model for the oxidation process, taking into consideration a non-zero water content for the segment when it enters the zone, and including heat transport by radiation. CONCLUSIONS The first-order DDP method used has solved the optimal control problem formulated in a satisfactory way and can be rec ommended as a tool for developing and verif y ing control strategies in this kind of pro c esses. The original model is simplified by the method of subdomains. The simplified non-linear mo del reproduces the states quite satisfact orily and can be used in a simulator as an operator~s guide thanks to its short computing time. In order to increase efficienc y and de c re a se heat losses, a hi gh bed in conjunction with a gas flow of small magnitude and a s high a gas temperature as possible must be used. The gas temperature is limited by the ma ximal heating rate for avoiding heat sh ock. REFE"ZENCES Batterham, R.J., J.A. Thurlby and G.J. Thornton (1977). Optimization of an iron-ore indurator. Chem. Eng. (GB) Sept. Breitholtz, C. (1982). Control models f o r distributed parameter proc esses by the method of subdomains. Te chni c al Re po rt No. 127 Chalmers University of Technology, Gothenbur g , Sweden. Breitholtz, C. and C. Hillberg (1980a). Development of a r e f e r e nc e mo del for the drying zone of a travelling grate pelletizing plant. Proc eed ings of 3rd IFAC Symposium on Automati on in Mining Mineral and Hetal Processing, Montrea l , Au g . 1980. Breitholt z , C. and C. Hillberg (1980b). Nonisothermal gas flow through the updraught drying zone of a travelling grate pelleti z ing plant. Rep ort R 80-01, Chalmers University of Technology, Gothenburg, Sweden, Feb. 1980 . Breitholtz, C. and C. Hillberg (1981). A simplified model of a heat- and massexchange process. Proceedings of BIAS 1981, International Conference on Control of Industrial Processes, Hilano , Oct. 1981 . Drugge, R. (1975). Modelling of Heat Exchange Phenomena During the Induration of Magnetite Pellet on a Travelling Grate. Preprint 75-6-51, Annual Meeting AIllli SME, New York, 1975. Drugge, R. (1981). Optimization o f pelletizing processes in sintering of iron ore concentrates at LKAB, Sweden. Preprint 3rd

382

C. Hillberg

International Symposium on Agglomerations, Nlirnberg, May 1981. Finlayson, B.A. (1972). The Hethod of Weighted Residuals and Variational Principles. Academic Press, New York. Hillberg, C. and B. Jarmark (1983). On pursuit-evasion between two realistic aircraft. SAAB TN Ae 75, Saab-Scania, Aerospace division, library, S-58l 88 LinkCiping, Sweden. Jacobson, D.H. and D.Q. Hayne (1970). Differential Dynamic Programming. Elsevier, New York. Szekely, J. and R.G. Carr (1968). On nonisothermal flow of gases through packed beds. Trans. of AIHE, Vo1. 242, pp. 918 - 921. Voskamp, J.H. and J. Brasz (1975). Digital simulation of the steady state behaviour of moving bed processes. Heasurement and Control, Vol. 8, Jan. 1975. Voskamp, J.H. and H. Vissers (1973). Optimization of a pellet indurating machine by means of simplified mathematical models. Proceedings of the Conference on 11athematical Process 110dels in Iron and Steelmaking, Amsterdam 19 - 21 Feb. 1973. TABLE 1

Notations Used Definition

I

Value

Uni t

Variables and functions Time-dependent parameters in series expansions

ai,bi,c i C(t,z)

Proportion of magnetite by weight in a pellet

C Tot Ein

C(t

t ) f Inlet energy

J/m

E

Outlet energy

J/m

ut Fi,Gi,Hi,K

i

=

2

2

Non-linear functions -1

KB(T ) s qox(t,z)

Reaction rate for oxidation

s

Heat released from oxidation

J m

R ,R ,R

Residuals

s

g

c

-3 -1 s

t

Time

r

Eff iciency factor

T (t,z) g

Temperature of gas

K

Ts(t,z)

Temperature of segment

K

T sb (t)

Bulk temperature of segment

K

u( t)

Temperature of gas at

u

Limi t value for

z

Height coordinate of the bed

m

Density of gas

-3 kg m

Viscosity of gas

Pa s

z

=

0 , control variable

u

K K

383

Differential Dyanmic Programming Definition

Unit

Heat transfer area between gas and pellet per unit volume of the bed

270.

-1 m

Heat capacitivity of pellets

880.

J kg

1193.

J kg

Process parameters

11

A c

Value

s

Heat capacitivity of gas Diameter of pellet

0.012

Gas flow through the zone

G

Height of bed

h =

o o

Constant in the reaction rate for oxidation Constant in the reaction rate for oxidation Molar weight for

=

1.8 0.35 0.5 6000. 0.231

Fe 0

3 4

450.

s

0.55

Heat of oxidation per

92.

Yp

Maximum heating rate

Fe 0

3 4

3.33

Ps

Density of pellet

- " PI

Pressure drop over the bed

0 < t < 315

-"P2

Pressure drop over the bed

315

E:

Void fraction of the bed

0.457

Shape factor of pellet

1.

4

Pressure drops (Pa)

0.35

" PI

2500

" P2

4200

" PI

2500

" P2

2500

" PI

3300

" P2

3300

" PI

4200

" P2

4200

" PI .6P 2

3300

" PI

3300

"P 2

3300

.6P l

2500

"P 2

4200

0.35

0.25

5

0.45

6

7 x)

x)

Bed height (m)

0.35

3

3350.

<

t

~

t

-2 -1 -1 J m K s -1 kJ mole -1 °c s -3 kg m Pa Pa

f

The Efficiencl r , the quotient of outlet-inlet energl, and the remaining proportion of magnetite for several bed heights and pressure drops

0.35

2

kg

Final time

-"H

1

K s

40.

Case No.

m -1 s

270.

Pellet volume per unit volume of the bed

TABLE 2

m -2 -1 kg m s

Time in the cost functional

Heat transfer coefficient between gas and pellets

Cl

-1 K- l -1 K- l

0.35

Eou/Ein

C Tot

0.54

0.65

0.55

0.59

0.60

0.55

0.55

0.66

0.54

0.51

0.70

0.52

0.38

0.80

0.47

0.64

0.49

0.60

0.44

0 . 60

0.62

r

3300

Same as Case 1, but the result applies to nominal control.

u(t)

min(u(t), 1050)

which is used as

c.

384

Hillb e r g

la

Fig. 1.

The segment .

.0.5

10

Magne tite cone. Fi g. 4.

234567 Ti me (m in.) Fig . 2 .

Pellet- temperature response to reasonable gas temperature vs . time .

o Fig . 5 .

..c

---..c" N

0.

<» "0

"0

.3 05 <»

>

ro

~

o o

~~

Fig . 3 .

__

- L_ _ _ _~_ _~_ _ _ _~~_ _~~~

500 1000 Pellet temperaturz (OC ) Relative bed depth vs. pellet temp temperature at various time lapses .

Re l ative bed depth vs . pellet temperature at various time lapses .

Q5

~ .O

Ma~net i t e con c.

Relative bed depth vs. magnetite concentration at various time lapses . A polynomial of the fifth degree is used .

385

Differential Dynamic Programming

10

--.c

.r. N

Case 4

0. <:>I

"0

Case 2

"0 <:>I

oD

0.5

<:>I

>

~

& 00

2 3 4 5 6 7 Time (min)

Posit i on in bed where maximal rate of Ts(t,z) occurs, vs. time. See also Table 2 .

Fig. 8.

01234567 Time (min) Fi g . 6 .

Solved contro l and pellet t emperatu r es vs . time. Case 3 .

10

10, ;'r

I / /

I '

/

/

~, --·-~W_"_III

900

1/

I; / / / /

j 11/ /J>~ 05

J t ! ! / I. I"~ 0234567 Time (min.)

Fig . 7.

\

3)=0 sca 7,90

Depths having the same pel let temperature . Case 3 .

__._

0.5 .Q

...~ III

o

U

o ~........-r-r""T"""T--r-r-r-~ o 1 2 3 4 5 6 7 8 9 10 Iteration number

Fig. 9.

Convergence pi c ture. Cos t /i niti al cost vs. itera ti on no .