Identification of Markov Matrices of Milling Models - American

Jul 16, 2009 - Detailed modeling of a grinding mill can be achieved through Markov ... In this study we model a complex multiregion mill using Markov ...
0 downloads 0 Views 749KB Size
Download PDF
Ind. Eng. Chem. Res. 2009, 48, 9763–9771

9763

Identification of Markov Matrices of Milling Models† Jayanta Chakraborty and Doraiswami Ramkrishna* School of Chemical Engineering, Purdue UniVersity, West Lafayette, Indiana 47907

Detailed modeling of a grinding mill can be achieved through Markov chain models without involving lengthy computations. However, estimation of the key parameters of the model, elements of the Markov transition matrix, using observable quantities is not trivial. This powerful modeling tool can find wide applicability in operation and control of industrial mills if this set of parameters can be estimated using suitably observed quantities. In this study we model a complex multiregion mill using Markov chain and propose a general technique for estimation of the Markov transition matrix for breakage problems. This technique estimates the transition matrix from observed evolution of particle size distributions in various regions of the mill, based on extracting spectral information of the transition matrix from the data. It has been shown in this study that a specific grouping of states can lead to lower triangular block structure of the Markov transition matrix for a breakage problem. Detailed analysis of such a block matrix reveals that simple semilogarithmic plots of a set of observed quantities can be used in order to extract the spectral information of the transition matrix which in turn can be used to reconstruct the Markov matrix. A numerical example has been presented to illustrate various ideas which demonstrate that the proposed technique can be used successfully to estimate the Markov transition matrix very accurately from observations. Introduction Milling of solid particles to produce smaller particles is an important unit operation in many industries, including the pharmaceutical industry. As such operations are highly energy intensive, there has been considerable effort in understanding these operations quantitatively. The earliest work in this direction is due to Rittinger, Bond, and Kick, who related the degree of size reduction to the energy consumed.1 In many practical situations, prediction of the size distribution of the product particles is important, especially where the operation produces finished or close to finished product. Such a situation calls for the use of the population balance framework in modeling of solid grinding.2 Kapur3 was among the first to use population balance models (PBM) for describing comminution processes. The key issue in those early PBM was extraction of the breakage functions from experimental data, one that is referred to as the inverse problem of PBM.4 Similarity behavior of population balance equation (PBE) often facilitates such estimation, and Kapur used a certain form for the breakage function and obtained a similarity solution of the associated PBE. Ramkrishna5 demonstrated an elegant way of estimating the breakage function for power law breakage using the similarity behavior. Later, Sathyagal et al.4 suggested a way to find the similarity solution for continuous breakage PBE for a more general breakage kernel. In another effort, Herbst and Fuerstenau6 discussed scale-up of a combined mill classifier system using PBM where breakage functions have been extracted from experimental data. More recently, Berthiaux and Varinot7 simplified the analytical solution of the batch grinding equation given by Kapur8 under certain approximations and extracted the breakage parameters for a laboratory-scale stirred bead mill. Although compact and † We are pleased to contribute to a festschrift for Dr. B. D. Kulkarni who has made special contributions to stochastic modeling in the chemical engineering literature. In our concern with the application of Markov chains in this paper, it represents an appropriate topic with which to felicitate Dr. Kulkarni. * To whom correspondence should be addressed. Tel.: +1765-4944066. Fax: +1765-494-0805. E-mail: [email protected].

simple, such models do not contain enough detail about particle movement and region-specific grinding inside a mill and therefore have found limited application on an industrial scale. An emerging approach involving detailed modeling of a particle breakage process is the discrete element method (DEM), in which particles are located in a finite domain and various dynamic events among particles are conducted through dynamic equations for those “elements” and breakage occurs according to empirical breakage laws.9 While such models take a more physical view of the comminution process, their usefulness would depend upon the development of efficient computational strategies. Vogel and Peukert10 have discussed the possibility of predicting breakage parameters for a particular material by testing a single particle under impact. Such models are useful for a reasonably homogeneous material. But in general industrial materials pose considerable variation in terms of mechanical properties necessitating the development of phenomenological models such as the ones in this work. Markov chain models represent an attractive alternative for modeling particulate processes without extensive computations. Such a Markov chain model has been formulated by Duggirala and Fan11 for attrition processes, and Berthiaux12 formulated Markov chain for comminution. In a Markov chain model, particles inside a mill are represented by a set of probabilities, called probability-state Vector, which may also be viewed as mass fraction of particles in discrete bins. The comminution process leads to transitions at fixed discrete intervals due to which the probability-state vector of the system changes according to a set of “transition probabilities”. In other words, the comminution process changes the mass fractions of particles during discrete time intervals. As discussed by Berthiaux et al.,13 the Markov chain model can emulate a realistic mill with simplicity. For example, multiple regions inside a mill are described under the same framework as a single-region mill only by defining states by referring to the size as well as the location of a particle. In a Markov chain model, every dynamic phenomenon is described in terms of the transition matrix, and hence its determination is

10.1021/ie900456j CCC: $40.75  2009 American Chemical Society Published on Web 07/16/2009

9764

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

very critical in the success of the model. Therefore, an effective method is required for extracting the transition matrix of the Markov chain from observed data, such as the particle size distribution at different times. Berthiaux et al.13 have discussed a method of estimating the approximate transition matrix from the population balance model. This approach is of course natural, and one expects the derivation to be appropriate at suitably small times. Alternatively, one may postulate that the behavior of the process may be approximated by a Markov chain model without recourse to its evolution from the more rational approach of population balance. Whether or not such an approximation is allowed is to be established by experiment. Theoretical discussion on estimation of a transition matrix from observed data is very sparse. Most analytical studies on Markov chains are based on prior knowledge of the transition matrix. Anderson and Goodman14 discussed a method for extraction of transition matrix for a discrete Markov chain where every individual element can be observed distinctly. This method has been used for extraction of the transition matrix of disease propagation models.15 In this method the likelihood estimation of the transition matrix is given by the expression pˆij )

∑ ∑

T n (t) t)1 ij T-1 n (t) t)0 i

(1)

where nij is defined as the number of individuals in state i at t - 1 and in j at t. Clearly, this method can be used only when each individual element in the system is tracked over time, which is the case for clinical trials; it cannot be adopted for particle breakage processes unless experimental methods are devised to follow individual particles. In this paper we describe a systematic procedure for extraction of the transition matrix from the evolution of particle size distributions. This method is based on the extraction of the spectral information of the transition matrix. It also takes advantage of the structure of the transition matrix that results from the proposed grouping of states. In the next section we present a brief discussion about Markov chains and their essential properties and then discuss about the Markov chain model of the breakage problem. It is followed by a discussion of how the mathematical properties of the Markov transition matrix can be used in order to extract the transition probabilities. Then we demonstrate our method using a numerical example before drawing conclusions toward the applicability of the framework to real processes. Markov Chain Model There are certain types of processes where the future state of the system does not depend on its past history and depends only on the present state. Such processes can be described mathematically using the concept of a Markov chain. We shall provide only a very brief discussion about Markov chain here, and the reader is referred to the text by Gnedenko16 for additional details. A Markov chain model comes under the framework of stochastic models,17 where system and processes are described in terms of probabilities. A system with various possible states is therefore represented through a probability vector which changes with time as the system evolves. In a discrete Markov process, such changes in probability vector (transitions) occur only at discrete instants of time denoted by τi. Such transitions happen through a set of probabilities known as transition

probability. These are the probabilities of each possible transition at a given instant. There will be a probability associated with the transition between each pair of states. This set of transition probability thereby forms a matrix which is known as transition matrix or transition probability matrix or MarkoV matrix, and in this work we shall denote it by P. Hence, if the initial state of the system, denoted by the vector m(0), is known, the subsequent states are determined through the transition matrix in the following way: m(1) ) P·m(0) m(2) ) P·m(1) ) P·P·m(0) ) P2·m(0) l m(k) ) Pk·m(0)

(2)

where m(k) indicates the state vector of the system at the instant τk. In general, the transition matrix can be time-dependent. But many practical situations can be modeled sufficiently by assuming it to be independent of time. In the current problem we shall model the system through a time-independent transition matrix. Such a Markov chain is known as a homogeneous Markov chain. The elements of the matrix P are the probabilities, and therefore it is also termed as stochastic matrix. The stochastic matrix has an interesting property17 that, irrespective of the numerical values of various probabilities, it will have one eigenvalue as unity and all other eigenvalues will be less than one. Mathematically, λ1 ) 1.0 > |λi |

(3)

Another interesting property of many Markov chains is the existence of an ergodic state.17 After many transitions, the state vector of the system often approaches a unique stationary state independent of P, which is known as the ergodic state. Mathematically, lim m(k) ) m(∞)

kf∞

(4)

The Markov chain model appears conceptually appropriate for modeling milling problems. Under this framework, the mill is described through a probability vector. Insofar as we use the frequency argument of probability to associate the probability vector with the mass fraction vector, the resulting model may be viewed as a deterministic model for the comminution process. For a smaller mill with no spatial nonhomogeneity, this can be accomplished simply by taking the size (mass) fraction vector as the state probability vector. But industrial mills are not wellmixed, and various grinding regions are often present. Materials flow from one region to another in a complex fashion while grinding occurs by several mechanisms. For example, if a mill can be characterized by three regions, namely, entry zone, main grinding zone, and exit zone, the entry region of a mill may experience only the grinding of very big particles. The main grinding region may be characterized by much vigorous grinding, and the exit region may experience abrasion and thereby produce finer particles. The process of grinding of solid particles is also complicated by the phenomenon of multiparticle interaction. In this situation a group of particles are ground together18 inside a mill. Bilgili et al.18 have modeled the case by considering nonlinear PBM where every particle present inside the mill affects the breakage of a particle. While this may be true for a small well-mixed mill, it overemphasizes the multiparticle effect for a larger mill

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

9765

Figure 3. Various possible (shown by solid arrows) and forbidden (shown by dashed arrows) transitions.

Figure 1. Description of the particulate system inside a multiregion mill in terms of the location and size of the particles. Particles in various regions are designated by various shades.

Figure 2. Pictorial representation of probability-state vector for particles in a multiregion mill in terms of the location and size of the particles. The probability-state vector is represented through various blocks where individual elements in each block are particles of the same size but from different regions.

as such multiparticle interaction is unlikely to percolate more than a few particle diameters. Therefore, one must consider various regions inside a mill before considering such a nonlinear effect. In this study we shall exclude such multiparticle effect and discuss a model with region-specific grinding and transport of material among various regions. Development of the Model Figure 1 shows a mill which can be characterized by three different grinding zones and the particles inside the mill can be classified into four different bins according to their size. Particles in various regions have been designated by different shades as shown in Figure 1. A randomly picked particle can be found in any of the three regions and may belong to any of the four size classes. Therefore the probability-state vector for this case is composed of 12 states. In general, if a mill can be characterized by Nr regions and Nb bins, the complete-state vector will contain Nr × Nb elements. The probability vector of such a multiregion mill can be divided into smaller blocks. The probabilities of a particle being in a particular size (of ith size class), but in various regions can be taken as the ith block which can be denoted as mi. For a general case, each block will be of length Nr and Nb such vectors together will compose the state probability vector, denoted by m which can be written in its column form as T m ) [mT1 mT2 · · · mNb ]T

(5)

where we have used the symbol T as a superscript to represent row vectors. This situation is shown pictorially in Figure 2 for the mill shown in Figure 1. Now, the dynamics of the system is determined by the transition matrix, which is the probability of transition between any pair of states. As the system is described by N ) Nb × Nr states, the transition matrix will be an N × N matrix. We denoted this matrix by P, and the ijth element of this matrix (pij) is defined as the probability of transition from state j to state i.

Figure 4. Transition probability matrix for breakage of particles for the mill shown in Figure 1. The forbidden transitions are shaded in black. Various blocks are marked with thick lines

Because the process considered here is only breakage, smaller particles are produced from bigger particles and the reverse is not possible. Such impossible events are characterized by zero probability of transition, and we shall term such events as forbidden transitions. Figure 3 presents various possible and forbidden transitions for the mill shown in Figure 1. In Figure 3, change in the size of a circle represents a breakage event and change in shade indicates movement of particles from one region to another. In this figure, the transition p5,11 is a forbidden transition (represented by a dotted line in Figure 3) as state 5 contains bigger size particles than state 11 and aggregation is not permitted. However, p12,4 is a possible transition (represented by solid line in Figure 3) as state 4 contains much bigger particles than that of state 12. This particular transition signifies simultaneous breakage and movement of particles as evident from the change in both size and shading in a single transition. A particle is also allowed to move from one region to another without breaking. One such transition is shown in Figure 3 as p3,1. The transition matrix for the example considered above is shown in Figure 4 where the forbidden transitions have been shaded in black. It is clear from Figure 4 that the transition matrix is a block matrix where the blocks are arranged in a lower triangular fashion. The block structure of this matrix is originated from the description of the probability-state vector through various blocks. In a Markov chain model, states can be numbered in any suitable way. For example, it is possible to represent the state vector as shown in Figure 5, instead of that of Figure 2. In the description of state shown in Figure 5, the individual blocks in the state vector contain the particle size distribution in a single region. Such numbering of state has been used by Berthiaux et al.12 in formulating a Markov chain model for the same breakage problem. Figure 5 also shows various possible transitions within a block (represented by a thin arrow) and among blocks

9766

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

Figure 5. Various possible transitions for the same mill as shown in Figure 1 but the states are grouped according to Berthiaux et al.12 The thin unidirectional arrows signify the breakage process, while the thick arrow represented the movement of particles among various bins.

where S is the modal matrix of P comprising the latter’s eigenvectors and Ξ is the diagonal matrix composed of the eigenvalues of P. Generally a mill can be observed at discrete time intervals, and hence the probability-state vector m(k) at any time τk can be observed. This can be done in an invasive way by drawing particles from each region through an inspection hole or in a noninvasive manner using advanced imaging technique such as X-ray tomography.19 Now, substituting eq 6 in eq 2 m(k) ) S·Ξk·S-1·m(0)

(7)

which can be written as N

m(k) i )

∑ s ·λ ·(S ij

k j

-1

·m(0))j

(8)

j)1

If the eigenvalues of the transition matrix are well-separated, i.e., |λ1 | . |λ2 | . |λ3 | . · · · . |λN |

(9)

the evolution at long time can be written in terms of the dominant eigenvalue as Figure 6. Location of forbidden transitions (shown in black) in the transition probability matrix when the states are numbered according to Berthiaux et al.12 as shown in Figure 5.

(represented by thick arrows). It can be noted that within each block, particles can transit only to their right due to breakage, but a particle can traverse between any two blocks in both direction as to-and-fro motion is possible within various regions in a mill. The corresponding transition matrix has been shown in Figure 6, marking the forbidden transitions in black. It is very interesting to note that in this case the lower triangular arrangement of blocks has disappeared. Instead, individual blocks have assumed a lower triangular shape. The reason for such difference is clear from Figures 3 and 5. The lower triangular feature of the matrix shown in Figure 4 resulted from the one-way transition between various blocks as shown in Figure 3, and as two-way transition within a block was permitted, each of the blocks were filled. On the other hand, the arrangement of blocks shown in Figure 5 is such that in each block there is only one-way transition possible, but bothways transition is possible among elements of two different blocks. The one-way transition in each block is reflected in the lower triangular nature of each of the blocks. But as transition among various blocks is possible in both ways, such lower triangular blocks appear everywhere in the matrix. In the next section we shall show that the transition matrix with lower triangular blocks is particularly helpful in extracting the transition matrix from observed data. Identification Strategy As discussed earlier, we shall show that the spectral information of the breakage transition matrix can be obtained from the observed quantity rather easily. Such information can be used to regenerate the transition matrix. In this section we shall describe how the spectral information of the transition matrix of the breakage problem can be obtained from the experimental data. The transition matrix can be written in terms of its eigenvalues and eigenvectors as P ) S·Ξ·S-1

(6)

k k -1 (0) m(k) i ) s1i·λ1·(S ·m )1 ) R1iλ1

(10)

Hence, a semilogarithmic plot of m(k) vs k will yield the i dominant eigenvalue as slope and R1i as intercept. Such plots for various components of m(k) will yield the same eigenvalue as the long-term limit, but different values of R1i. All such values of R1i can be represented through a vector which is proportional to the eigenvector corresponding to λ1. We shall denote this vector as r1. At this point, we notice that as the left-hand side of eq 10 is a probability; it is always positive for any value of k. Therefore, the first dominant eigenvalue λ1 and the corresponding r1 both must be positive. After we extract the first dominant eigenvalue, the next prominent eigenvalue can be obtained by peeling off the dominant eigenvalue using k k -1 (0) -1 (0) m(k) i - s1i·λ1·(S ·m )1 ) s2i·λ2·(S ·m )2

(11)

We shall call the quantity on the left-hand side of eq 11 the first residue and denote it by r(k) 1i . Hence, in terms of this notation, eq 11 can be written as k r(k) 1i ) R2iλ2

(12)

Since the above residue can be negative, its absolute value may be used to plot against k on a semilogarithmic scale. k |r(k) 1i | ) |R2i |·|λ2 |

(13)

Hence, using eq 13, the absolute value of the second dominant eigenvalue and the corresponding eigenvector can be extracted. It may appear from eq 13 that the sign of the eigenvalue and the corresponding R cannot be extracted through this method. But the signs of the corresponding quantities can be recovered by simple observation of the data. We notice that if the residue (k) has the same sign irrespective of k, λ2 must be positive and r1i R2i has the sign of the residue. On the other hand, if this residue is showing flipping sign with each increment in k, λ2 must be negative and sign of R2i is the same as the sign of the residue for even values of k.

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

Such a peeling-off procedure can be continued until all the eigenvalues and the corresponding ri are extracted. All such ri will produce a matrix, and we shall denote this matrix by Α. It can be shown that due to biorthogonality of eigenvectors, scalar multiples of the eigenvectors can be used instead of eigenvectors in eq 6, and hence P can be reconstructed as P ) A·Ξ·A-1

9767

This set of equations can be written as

(14)

This procedure depends on the separation of the eigenvalues. Unfortunately, this possibility is a very fortunate circumstance and it would seem that this method would fail if the eigenvalues are not suitably separated. In the following analysis we shall show that the above procedure can be utilized even if the eigenvalues of the transition matrix are not separated well. Analysis. As shown in Figure 4, the transition matrix for pure breakage process can be written as

Now, if λ is an eigenvalue of LNbNb, it will result in a trivial solution for x1. As x1 is a null vector, and λ is not an eigenvalue of L22 either, x2 will also be a null vector. Similarly, all other xis except xNb will be null. Clearly, xNb will be the eigenvector of LNbNb corresponding to eigenvalue λ. Hence, the eigenvector of P corresponding to λ can be written as [0T 0T · · ·

T · · · 0T xNb ]T

(19)

Here 0 denotes the null vector of length Nr. As LNbNb will have Nr eigenvalues, there will be Nr eigenvectors with such structure. In a similar way, if λ is any of the eigenvalues of LNb-1,Nb-1, x1,x2,..., xNb-2 will be null and xNb-1 will be the eigenvector of LNb-1,Nb-1 and xNb is given by LNb,Nb-1xNb-1 + (LNb,Nb - λI)xNb ) O

(20)

Hence, the matrix of eigenvectors can be written as where Lij(Nr × Nr) is a matrix representing the ijth block of the transition matrix. The characteristic polynomial for the above matrix can be written through the determinant equation:

Here Xii are the matrices of the eigenvectors of Lii. It can be easily shown that the inverse of S assumes the same shape as S, and we shall denote this matrix as Clearly, this equation is equivalent to |L11 - λI|·|L22 - λI|· · · · ·|LNbNb - λI| ) 0 And, hence, the eigenvalues of the matrix P forms smaller groups where each group corresponds to each diagonal block. There will be Nr eigenvalues in each group, and Nb such groups will exist. Now, the equation for eigenvectors can be written as The matrix of eigenvalues is written as

where the eigenvector corresponding to the eigenvalue λ has been represented through Nb smaller blocks of length Nr.

9768

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

where Λi is a diagonal matrix composed of the eigenvalues of Lii. Now, the general dynamic equation m(n) ) SΞnS-1m(0)

(24)

can be written as

Figure 7. Various possible transitions responsible for the elements in the below diagonal blocks of the transition matrix shown in Figure 4.

composed of the scalar multiples of the eigenvector. Hence, Ljj can be reconstructed using Ljj ) Ajj·Λj·Ajj-1 where the last three terms of eq 24 have been combined into the vector γ given by i

γi ) Λni

∑Y m ij

(0) j

(26)

j)1

It can be seen from eq 25 and eq 26 that, due to the lower triangular block feature of the matrix S, only Λ1 is involved in the evolution of m1. Hence, as long as the eigenvalues of the associated matrix (L11) are well-separated, they can be extracted with accuracy along with the corresponding r using the evolution of m1 in the way discussed earlier. In this case we need not to refer to any other eigenvalues of P. We notice that in general as the number of regions in a mill is much smaller than the number of bins, the individual blocks are much smaller in dimension than the overall matrix. It is much more probable that the eigenvalues of each block are separated enough. Now, because both Λ1 and Λ2 are involved in the evolution of m2, they cannot be extracted in that way unless they are separated enough. At this point, we notice that if γ1 can be set to zero, we can eliminate the part that contains Λ1 from the dynamics of m2 and only Λ2 remains. As evident from eq 26, γi depends on the initial feed and can be set to zero by manipulating the feed. For example, γ1 can be set to zero simply by setting m1(0) to zero and then Λ2 can be extracted accurately by plotting the dynamics of m2(k). In general, to set γi to zero, we need to set m1(0) through mi(0) to zero. Mathematically, if the highest size present in a feed is mj(0), γk ) 0 k < i

A series of such feed then can be used to extract all the diagonal blocks of the transition matrix. After we evaluate all the diagonal blocks, we need to evaluate the matrices situated below the diagonal. We first examine the physical processes responsible for various transition probabilities in the below diagonal matrices. Such transitions are shown in Figure 7. The transition shown by the solid arrow represents processes in which the daughter particles formed through breakage are found in the same region as the parent particle. Such transitions are located on the diagonals of the below diagonal blocks. The off-diagonal elements of those blocks represent simultaneous breakage and movement of particles in a single transition, as shown by the dashed arrow in Figure 7. Such transitions are possible only at the boundary of the regions and for most practical purposes can be neglected. Hence, the block matrices below the diagonal blocks can be considered as a diagonal matrix. Now, we examine the first transition using a feed that contains particles of the jth size class only. This transition can be written as

(27)

k

γk ) Λnk

∑Y m kj

(30)

(0) j

Simplifying,

kgi

j)i

This situation is further simplified if a feed of only one size is used. For a feed of size mj(0) γk ) 0 k < j γk ) Λnk Ykjm(0) kgj j

(28)

Hence, with this feed, the dynamic equation reduces to

[][ ]

(32)

mj+i ) Lj+i·m(0) j

(33)

mj mj+1 l mNb

Ljj·m(0) j

(1)

)

Lj+1j·m(0) j l LNbj·m(0) j

Writing for the ith block,

Because Lj+i is a diagonal matrix, this can be written as (0) Lj+i, j(k, k) ) m(1) j (k)/mj (k)

and therefore, the evolution of mj can be followed in the same way as discussed in the previous section in order to extract the eigenvalues Λj and the corresponding matrix Ajj which is

(34)

Therefore, all of the off-diagonal blocks below the diagonal block Ljj can be estimated using the same feed that was used for estimation of Ljj. Hence, the recommended strategy for

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

9769

Figure 8. Example of a breakage transition matrix where region-specific breakage has been considered. Particles have been assumed to move only to the neighboring region in a single transition.

extraction of the transition matrix is the following: (1) Feed the mill with particles of any particular bin size (mj(0)); (2) follow the evolution of mj and estimate Ljj; (3) use the first transition of the previous observation to estimate Lj+1,j through Ln,j; (4) use a series of such feed to estimate the entire matrix. Numerical Example. In this section we shall demonstrate the strategy discussed in the previous section with numerical example. We consider the transition matrix as shown in Figure 8 written for a mill with three regions and particles classified into three bins. In this case, the transfer of particle only to the neighboring bins has been considered which give rise to the tridiagonal matrices on the diagonal blocks. We also considered region-specific grinding in this example, and for this reason all the below diagonal blocks are different. This particular matrix has all positive, negative, and zero eigenvalues and also it has few closely spaced eigenvalues. The eigenvalues of this matrix are [0.6500; 0.2266;

Figure 9. Evolution of the particle-state vector toward the ergodic state according to the transition matrix shown in Figure 8 and the initial-state vector given by eq 35.

- 0.1766; 0.7000; 0.4000; 0.0000; 1.0000; 0.6822; 0.1246]

For this being the underlying transition matrix, the observed evolution can be constructed by following eq 2 for any given initial probability vector. For example, if we consider feed of the largest size particle distributed among the three regions with equal probability, as given by eq 35, m(0) ) [0.3333; 0.3333; 0.3333; 0.00; 0.00; 0.00; 0.00; 0.00; 0.00]

(35)

the evolution leads to production of the finest particle distributed among various regions as shown in Figure 9. We notice that the system reaches a stationary state at the later stage of evolution which is known as the ergodic state for a Markov chain. The state vector for the ergodic state is given by m(∞) ) [0.00; 0.00; 0.00; 0.00; 0.00; 0.00; 0.143; 0.428; 0.428]

(36)

Now, according to the proposed strategy, we need to feed the mill with particles of various sizes in order to extract the transition matrix. In this demonstration, we multiply the appropriate initial-state vector with the transition matrix using eq 2 to generate the observed evolution. This evolution can be used to regenerate the transition matrix. We begin with particles of smallest size distributed among various regions. The initial-state probability vector for this case is given by m(0) ) [0.00; 0.00; 0.00; 0.00; 0.00; 0.00; 0.33; 0.33; 0.34]

(37)

Figure 10. Evolution of the particle-state vector toward the ergodic state according to the transition matrix shown in Figure 8 and the initial-state vector given by eq 37.

Since this is the feed in the smallest possible size, the corresponding evolution (shown in Figure 10) will involve only the redistribution of the particles among various regions. Such redistribution reaches an ergodic state as shown in Figure 10 and is given by m(∞) ) [0.00; 0.00; 0.00; 0.00; 0.00; 0.00; 0.143; 0.428; 0.428]

(38)

We notice that the ergodic state given by eq 38 is exactly the same as that obtained with a different feed (eq 36) previously. This is an inherent property of Markov chain that the ergodic state does not depend on the initial feed and is a characteristic of the transition matrix. As discussed previously, the eigenvalue corresponding to ergodic states is 1.0 and the corresponding R values are given by the third block of the ergodic state vector. Hence, Λ3(1, 1) ) 1.0 r3(1) ) [0.143; 0.428; 0.428 ] Now we need to peel off this part from the overall dynamics using eq 11 and plot the residue in logarithmic scale in order to extract the next prominent eigenvalue as in eq 13. A typical

9770

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009 Table 2. Comparison of the Actual and Estimated Eigenvalues for L33 actual EV

estimated EV

1.0 0.6828 0.1172

1.0 0.6824 0.1332

Table 3. Comparison of the estimated L33with the actual matrix. The actual values are given in bold 0.42 (0.4) 0.60 (0.6) 0.01 (0)

0.20 (0.2) 0.60 (0.6) 0.20 (0.2)

0.00 (0) 0.20 (0.2) 0.80 (0.8)

Table 4. Comparison of the Actual and Estimated Eigenvalues for L11

Figure 11. Plot of Ln(m(k)(7) - m(∞)(7)) vs k for extraction of the second dominant eigenvalue of L33. Table 1. Eigenvalues and the Corresponding Intercepts (r) Obtained Using the Evolution Shown in Figure 10 intercepts (R) size used class used F F F

7 8 9

EV obtained

1

2

3

1.0/0.6822/0.1246 0.1428 0.0595 0.1217 1.0/0.6828/0.1140 0.4286 0.0831 -0.1811 1.0/0.6825/0.1353 0.4286 -0.1425 0.0522

actual EV

estimated EV

0.65 0.2266 -0.1766

0.65 0.2257 -0.1749

The comparison of the estimated and exact eigenvalues are shown in Table 2, and we note an excellent agreement between them. Now the eigenvalues (Λ3) and the corresponding matrix A33 together can be used to reconstruct L33 using eq 30, which has been shown in Table 3. The exact values of the elements are given in the bracket for easy reference. We notice that the matrix could be reproduced with remarkable accuracy. Other diagonal blocks can also be extracted using the same procedure. The estimated eigenvalues for the diagonal block L11 has been shown with the exact values in Table 4. As expected, the estimated EVs show excellent agreement with the exact values. We notice that the third eigenvalue in this case is negative and the corresponding residual should show a sign flip with

Figure 12. Plot of Ln(m(k)(7) - m(∞)(7) - R2λk2) vs k for extraction of the third eigenvalue of L33.

plot of Ln(|m(k)(i) - m(∞)(i)|) vs k is shown in Figure 11. We note that the plot shows a remarkable linear feature almost over the entire range of k demonstrating dominance of an eigenvalue in that region. After we extract the eigenvalue as the slope and corresponding R as intercept from Figure 11, the values of which are shown in Table 1, we need to peel off this part again in order to extract the last eigenvalue. The corresponding plot is shown in Figure 12. Again we notice an excellent linear feature here demonstrating the success of the peeling-off process. We also note that the third and hence the smallest eigenvalue die very quickly, as also evident from Figure 11. The extracted eigenvalues and the corresponding R for various size fractions are summarized in Table 1.

Figure 13. Characteristic oscillation in the residual observed for the negative eigenvalue of L11. Table 5. Comparison of the Estimated L11 with the Actual Matrixa 0.12 (0.10) 0.52 (0.55) 0.01 (0) a

0.10 (0.10) 0.14 (0.15) 0.41 (0.40)

0.00 (0) 0.21 (0.20) 0.44 (0.45)

The actual values are given in bold.

Table 6. Comparison of the Estimated L22 with the Actual Matrixa 0.20 (0.2) 0.50 (0.5) 00 (0.0) a

0.10 (0.1) 0.30 (0.3) 0.30 (0.3)

The actual values are given in bold.

0.00 (0.0) 0.10 (0.1) 0.60 (0.6)

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

9771

homogeneous Markov chain model can describe grinding in the mill. In this regard, it would be of interest to see whether a specific transition matrix can be applied to a uniform spacing of discrete times of particle size measurements. It seems reasonable to conjecture that particle size distributions from a mill could resist a homogeneous Markov chain model because of the varying environment in which comminution occurs as the process progresses in time. For example, as Bilgili18 suggests, break-up of a particle could result from interaction with its neighbors. It is possible that such a process could be intertwined within a Markov chain model formulation by allowing the transition matrix to depend on the state vector. Addressing this and other issues of formulation and model identification represent the broad direction of our future effort. Figure 14. Estimated transition matrix.

changing values of k. Such behavior could be observed distinctively, as shown in Figure 13, and hence negative eigenvalues can be identified easily with the current strategy. The reconstructed L11 has been shown in Table 5, and we notice an excellent agreement in this case also. The diagonal block L22 has a zero eigenvalue which can also be estimated using the current strategy. Such an eigenvalue can be easily identified as the residual becomes zero before all the eigenvalues are extracted. But extraction of the corresponding R needed some attention. Because the first two eigenvalues are nonzero, they can be extracted in a regular way and the residual becomes (k)

(m (i) - R1λk1 - R2λk2) ) R3λk3

(39)

Taking limits, for k ) 0, (0)

lim (m (i) - R1λ01 - R2λ02) ) lim R3λ03

λ3f0

(0)

λ3f0

(40)

m (i) - R1 - R2 ) R3

Hence, the value of R corresponding to the zero eigenvalue is simply given by the value of the residual for k ) 0 and L22 is obtained as given in Table 6. The below diagonal blocks can also be extracted with accuracy using eq 34 and the complete matrix has been produced in Figure 14, which can be compared with Figure 8. Conclusions and Outlook for Future We have demonstrated a strategy for determining the transition matrix of a Markov chain model toward its application to complex mills assuming that its environment can be represented by multiple grinding zones. By grouping the mass distribution vector in terms of its distribution in different regions for each bin size, the Markov transition matrix assumes a triangular block structure. Thorough analysis of such a transition matrix reveals that such lower triangular shape can be exploited in order to estimate the transition matrix from observed data. Simple logarithmic plots of a set of transient data have been shown to yield the spectral information of the transition matrix which in turn has been used to reconstruct the matrix. A numerical example has been considered in order to illustrate the methodology and its capacity for accurate identification of the transition matrix with accuracy. Whether or not the methodology can deal with real experimental data from a milling process depends upon how well a

Acknowledgment We gratefully acknowledge funding for this research from the NSF Engineering Research Center for Structured Organic Particulate Systems. Literature Cited (1) Austin, L. G. A commentary on the Kick, Bond and Rittinger laws of grinding. Powder Technol. 1973, 7, 315. (2) Ramkrishna, D. Population Balances: Theory and Applications to Particulate Systems in Engineering; Academic Press: San Diego, 2000. (3) Kapur, P. C. The energy-size reduction relationships in comminution of solids. Chemical Eng. Sci. 1971, 26, 11. (4) Sathyagal, A. N.; Ramkrishna, D.; Narsimhan, G. Solution of inverse problems in population balances-II. Particle break-up. Comput. Chem. Eng. 1995, 19, 437. (5) Ramkrishna, D. Drop breakage in agitated liquid liquid dispersions. Chem. Eng. Sci. 1974, 29, 987. (6) Herbst, J. A.; Fuerstenau, D. W. Scale-up procedure for continuous grinding mill design using population balance models. Int. J. Min. Proc. 1980, 7, 1. (7) Berthiaux, H.; Varinot, C. Approximate calculation of breakage parameters from batch grinding tests. Chemical Eng. Sci. 1996, 51, 4509. (8) Kapur, P. C. Kinetics of batch grinding: part-B: An approximate solution to the grinding equation. Trans. AIME 1970, 247, 309. (9) Morrison, R. D.; Cleary, P. W. Using DEM to model ore breakage within a pilot scale SAG mill. Miner. Eng. 2004, 17, 1117. (10) Vogel, L.; Peukert, W. From single particle impact behavior to modeling of impact mills. Chem. Eng. Sci. 2005, 60, 5164. (11) Duggirala, S. K.; Fan, L. T. Stochastic analysis of attritionsA general cell model. Powder Technol. 1989, 57, 1. (12) Berthiaux, H. Analysis of grinding process by Markov chains. Chemical Eng. Sci. 2000, 55, 4117. (13) Berthiaux, H.; Mizonov, V.; Zhukov, V. Application of the theory of Markov chains to model different processes in particle technology. Powder Technol. 2005, 157, 128. (14) Anderson, T. W.; Goodman, L. E. Statistical inference about Markov chains. Ann. Math. Stat. 1957, 28, 89. (15) Craig, B. A.; Sendi, P. P. Estimation of the transition matrix of a discrete-time Markov chain. Health Econom. 2002, 11, 33. (16) Gnedenko, B. V. The Theory of Probability; Mir: Moscow, 1969. (17) Cox D. R.; Miller H. D. The Theory of Stochastic Processes; John Wiley & Sons: New York, 1965. (18) Bilgili, E.; Scarlett, B. Population balance modeling of non-linear effects in milling process. Powder Technol. 2005, 153, 59. (19) Wang, L. B.; Frost, J. D.; Lai, J. S. Three-dimensional digital representation of granular material microstructure from X-ray tomography imaging. J. Comput. CiVil Eng. 2004, 18, 28.

ReceiVed for reView March 20, 2009 ReVised manuscript receiVed June 2, 2009 Accepted June 11, 2009 IE900456J