Hierarchical batch-to-batch optimization of cobalt oxalate synthesis process based on data-driven model

Chemical Engineering Research and Design 1 4 4 ( 2 0 1 9 ) 185–197

Contents lists available at ScienceDirect

Chemical Engineering Research and Design journal homepage: www.elsevier.com/locate/cherd

Hierarchical batch-to-batch optimization of cobalt oxalate synthesis process based on data-driven model Runda Jia a,b,∗ , Zhizhong Mao a,b , Dakuo He a,b , Fei Chu c a

School of Information Science & Engineering, Northeastern University, Shenyang 110004, China State Key Laboratory of Synthetical Automation for Process Industries, Northeastern University, Shenyang 110004, China c School of Information and Control Engineering, China University of Mining and Technology, Xuzhou 221116, China b

a r t i c l e

i n f o

a b s t r a c t

Article history:

The synthesis process has been widely used in cobalt hydrometallurgical industry. To better

Received 4 October 2018

operate the cobalt oxalate synthesis process, a data-driven model based hierarchical batch-

Received in revised form 17 January

to-batch optimization method is presented in this work. In the upper level of hierarchy, the

2019

proposed response surface model based modifier-adaptation (MA) strategy is used to calcu-

Accepted 31 January 2019

late the nominal control profile for the next level, and the design of dynamic experiment

Available online 8 February 2019

(DODE) method is also employed to symmetrically generate the dataset for response surface model building. In the lower level of hierarchy, the batch-wise unfolded PLS (BW-PLS) model

Keywords:

based self-tuning batch-to-batch optimization method is utilized to further refine the con-

Batch processes

trol profile on the basis of the result of the upper level. The main advantages of the proposed

Data-driven model

method are: (i) the size of the dataset for data-driven model building are rather modest, (ii)

Cobalt oxalate synthesis process

the control profile can be discretized into a large number of intervals to further improve the

Hierarchical batch-to-batch

optimization performances, and (iii) the unqualified batches is efficiently avoided during the

optimization

evolution of batch-to-batch optimization. The superior performances for the cobalt oxalate synthesis process are verified through simulation study. © 2019 Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

1.

Introduction

update, and (iii) optimization with model-based gradient. Since MFO (i) needs to perform the optimization on a trial and error basis, it typ-

Batch crystallization is increasingly important as they are widely used

ically converges slowly, and is more suitable for the batch processes

in many fields such as polymer, metal and pharmaceutical production (Corbett and Mhaskar, 2016). Compared to the continuous processes,

with a short cycle time and low operational cost (Kong et al., 2011; Zhu et al., 2014). When the model-based optimization methods (ii and iii) are used, the batch process model becomes critical.

batch processes exhibit a number of characteristics that lead to interesting control and optimization problems, some of which are: high nonlinearity, complex time dependent dynamics, and the end-product

During the last decade, several literatures have researched on first principal model based batch crystallization process control and optimization. For example, Shi et al. (2006) focused on predictive control

qualities are only available after the finish of the batches. In order to derive the maximum benefit from batch processes, the batch-to-batch

of particle size distribution (PSD) based on population balance models;

optimization method is based on iteratively updating the control profile using the information obtained from previous batch runs, and which can be grouped into the following different types (Srinivasan et al.,

multi-objective optimization was employed by Sarkar et al. (2006) to optimize the seeded batch crystallization processes for the best operation; Mesbah et al. (2010) proposed to apply dynamic optimization for

2001): (i) model-free optimization (MFO), (ii) optimization via model

maximizing the throughput of a semi-industrial batch crystallization

∗

Corresponding author at: School of Information Science & Engineering, Northeastern University, Shenyang, 110004, China. E-mail address: [email protected] (R. Jia). https://doi.org/10.1016/j.cherd.2019.01.032 0263-8762/© 2019 Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

186

Chemical Engineering Research and Design 1 4 4 ( 2 0 1 9 ) 185–197

process. However, the effort and time required to develop a comprehensive knowledge-drive model for batch crystallization generally carry a high cost, and almost all the first principle model based optimization methods are complex and computationally intensive (Fiordalis and Georgakis, 2013). For the end of minimizing the requirement of detailed first principle model, Hermanto et al. (2011) modeled the polymorphic crystallization process by combining a simplified first principle model with an updated PLS model. And the batch-to-batch and nonlinear model predictive control for polymorphic transformation was performed based on the proposed hybrid model. A similar batch-tobatch control of PSD in cobalt oxalate synthesis was also presented by Zhang et al. (2012). Although both of the above methods enable good performances, they are only suitable for systems with well built knowledge-driven model and it is still a challenge to implement these methods to the systems with little prior knowledge (Jia et al., 2015). Data-driven modeling methods address the issues of knowledgedriven or hybrid model based methods, because they can use the existing information routinely collected from batch crystallization processes and give an opportunity to easily calibrate and maintain these models (Flores-Cerrillo and MacGregor, 2003). To make the data-driven model more reliable, Zhang (2004) presented bootstrap aggregated neural networks based batch optimization method. However, since the control profile was optimized offline only once, the solution may result in suboptimal operation due to the plant-model mismatch. In order to cope with this issue, Zhang et al. (2009) further used an iterative learning control (ILC) method for maintaining a desired supersaturation profile based on batch-wise unfolded partial least squares (BW-PLS) models in a simulated crystallization process. A hierarchical ILC strategy for systematic design of the supersaturation controller for seeded batch cooling crystallizers was proposed by Sanzida and Nagy (2013). In their work, the upper level of hierarchy aimed to determine the constant supersaturation needed to produce the desired crystal, and the lower level of hierarchy was designed to maintain the supersaturation at a certain set-point. Although these results offer a substantial modification to the existing data-driven model based ILC methods, the desired control objectives are required to be determined in advance (Forgione et al., 2012), and the control profiles are usually discretized into a relative small number of intervals. Dong et al. (1996) first proposed a batch-to-batch optimization approach by using a neural net multi-way PLS (NNMPLS) model, and the process data was worked as feedback to correct the gradients computed from the NNMPLS model. Nevertheless, the NNMPLS model used in their approach was fixed, and the way to correct the gradients was not addressed in their work. Recently, batch-to-batch optimization of cobalt oxalate synthesis process is proposed by Jia et al. (2015). The latent variable model was used to describe the causal relationship between the control profile and final mean crystal size, and the modifier-adaptation (MA) strategy is employed to cope with the plantmodel mismatch. And they also proved that the BW-PLS model based ILC is equivalent to the proposed batch-to-batch optimization method by using only zero-order modifier (Jia et al., 2016). However, for some batch processes with high operational cost, it is impractical to generate a large number of initial samples for BW-PLS model building.

rate are generally preferably, which can significantly reduce the operational cost. Finally, to prevent the arising of the unqualified batches, the final cost of each batch should better be continuously reduced on the basis of a better start. To solve these issues, a hierarchical batchto-batch optimization method is proposed in this work. At the upper level of hierarchy, design of dynamic experiment (DODE) proposed by Georgakis (2013) and Klebanov and Georgakis (2016) is used to systematically generate the dataset for response surface model building. Then MA strategy is utilized to cope with the problem of plant-model mismatch, and the control profile is further improved at this level. At the lower level of hierarchy, the self-tuning batch-to-batch optimization method continues to reduce the final cost of batch process on the basis of a better nominal control profile. Not only can the convergence rate be greatly accelerated, but also the unqualified batches can be efficiently avoided. The layout of this work is organized as follows: In Section 2, the first principle model of cobalt oxalate synthesis process is introduced, and the DODE method and BW-PLS model are reviewed in brief. In Section 3, followed by the presentation of the hierarchical batch-to-batch optimization architecture, the MA strategy based on response surface model and self-tuning batch-to-batch optimization method based on BW-PLS model are elaborated in detail. Section 4 verifies the advantages and the effectiveness of the proposed method through comparisons among different batch-to-batch optimization methods. Finally, Section 5 draws some conclusions.

2.

Preliminary materials

First, a description of the cobalt oxalate synthesis process is provided below, followed by the description of the first principle model used in this work, and the issues of operating the synthesis process. Then the design of dynamic experiment (DODE) method and the batch-wised unfolded PLS (BW-PLS) model are reviewed briefly.

2.1.

Cobalt oxalate synthesis process

Metal cobalt powders are essential for the production of many types of products, such as abrasion strengthened composites and alkaline rechargeable batteries (Mao et al., 1999; Roure et al., 1999). Fed-batch cobalt oxalate synthesis process is a major operation unit commonly found in cobalt hydrometallurgy industrial, and the cobalt oxalate powders are produced in relative low-volume batch reactors. Then the following reaction takes place in the crystallizer. CoCl2 + (NH4 )2 C2 O4 → CoC2 O4 ↓ +2NH4 Cl

(1)

For batch processes with high operational cost, such as cobalt oxalate synthesis, the following issues should be seriously considered in designing the batch-to-batch optimization method. In order to min-

The cobalt chloride is reacted with the ammonium oxalate in liquid phase, and further leads to the desired cobalt oxalate crystals. Due to the complex characteristics of the synthesis process, a fed-batch mode is usually used to operate such unit. Firstly, the ammonium oxalate solution is stored in the dissolver and heated as shown in Fig. 1; then a fixed volume of the cobalt chloride solution is added into the crystallizer; finally, the ammonium oxalate solution is gradually fed and the reaction takes place until the batch termination. In order to maintain the constant temperature during the reaction, a heating jacket with PI controller is utilized. Since the reaction rate in the crystallizer is very quickly, the solution concentration C can be described from the view of mass conservation as follows:

imize the experiment cost, a systematic way should be firstly used to efficiently generate the dataset for data-driven model building. Secondly, a better initial (nominal) control profile and faster convergence

  FC dC c0 v0 F −  3G2 + Br03 − = 2 dt V (v0 + Ft)

In addition, moderate number of intervals should also be selected for discretizing the control profile to balance the final optimization result and the convergence rate. Recently, Camacho et al. (2007, 2015) proposed a self-tuning batch-to-batch optimization method for fed-batch processes by using unfolded PLS model. In their method, the feeding profile was discretized into much more intervals, while the number of trials did not require reaching the number of intervals in the feeding profile. However, compared to the MA strategy based batch-to-batch optimization method, the convergence rate was slowed down by gradually adapting the local data-driven model to the current operating point of the batch process.

(2)

Chemical Engineering Research and Design 1 4 4 ( 2 0 1 9 ) 185–197

187

Fig. 1 – The schematic diagram of the cobalt oxalate synthesis process. where c0 is the concentration of ammonium oxalate solution, and v0 is the initial volume of cobalt chloride solution, r0 is the initial crystal size at nucleation,  is the density of crystal, and  is the volumetric shape factor. G and B are the growth and nucleation rate respectively, which can be calculated using the following equations: G = G e−KG /T n (C)

G

solution unchanged during the batch evolution. Recently, we have discretized the feed rate profile into moderate number of intervals wherein it remains constant (Jia et al., 2015). Apparently, more intervals can further improve the optimization performances. However, because of the difficulty in estimating the experimental gradients, a large number of intervals are impractical in industrial synthesis process.

(3)

2.2. B = B e−KB /T (C)

B

where G and G are the growth rate and exponent coefficients, B and B are the nucleation rate and exponent coefficients. T is the temperature in the crystallizer, n is the agitation speed, and C is the supersaturation. KG , KB and  are three constant parameters. The suspension volume V in Eq. (2) can be simply modeled as follows: dV =F dt

(5)

where F is the feed rate of ammonium oxalate solution. In order to transform the population balance model into a set of ordinary differential equations (ODEs), the method of moments is used in this work. Based on the definitions for each moment of the PSD, the 2nd moment 2 in Eq. (2) can be calculated by solving the following set of ODEs: d0 F0 =B− dt V dj dt

(6)

j

= jGj−1 + Br0 −

Design of dynamic experiment (DODE)

(4)

Fj V

For batch processes, the control profile generally has a great effect on the end-product quality. Under this case, the classic design of experiment (DOE) method does not provide a systematic and efficient means to uncover the impact of a time-varying manipulated variable on the process performance. In order to solve this issue, Georgakis (2013) proposed a method to design of dynamic experiment (DODE). In addition, integrating DODE with response surface methodology (RSM) can describe the features of a batch process response through a small number of experiments. In order to generate the dataset for response surface modeling, the control profile is first defined by dimensionless variables in DODE as follows: N 

u() = u0 () + u()

(8)

where  = t/tf is the dimensionless time, and tf is the final time of a batch. u0 () and u() can be calculated using the following equations:

(7)

Since the characteristics of PSD can be described by using a function of the above four lowest order moments, j = 1, 2, 3 is considered in this work. Due to the absence of on-line instrument for principal variable, the plant operators usually set the feed rate of ammonium oxalate solution as invariant value. On the other hand, to ensure that all the cobalt chloride is exhausted, an excessive amount of ammonium oxalate is added into the crystallizer. However, the control profile has a significant impact on the end-product quality. Thus, to avoid solving a hybrid model based constrained dynamic optimization problem, Zhang et al. (2012) still kept the feed rate of ammonium oxalate

˛i ϕi ()

i=1

u0 () =

u() =

umax () + umin () 2 umax () − umin () 2

(9)

(10)

where umax () and umin () are the maximum and minimum bounds of the manipulated variable, respectively. ϕi () is an appropriate basis function defined between 0 and 1, where i = 1, 2, ..., N. The coefficient ˛i is the ith dynamic sub-factor (DSF) just as the coded variable in DOE, and can be used as input variables in the RSM. The control profile can be parameterized as an expansion of the shifted Legendre polynomials (Fiordalis and Georgakis,

188

Chemical Engineering Research and Design 1 4 4 ( 2 0 1 9 ) 185–197

2013). Since the values of the shifted Legendre polynomials are between −1 and 1 when 0 ≤  ≤ 1, the following inequality constraints should be guaranteed. −1 ≤ ˛1 ± ˛2 ± · · · ± ˛N ≤ +1

N 

ˇˆ i ˛i +

i=1

N N  

ˇˆ ij ˛i ˛j +

i=1 j=i+1

N 

ˇˆ ii ˛2i

(12)

i=1

Batch-wise unfolded PLS (BW-PLS) model

Generally, a completely knowledge-driven model for batch industrial crystallization is not available for process control or optimization. In a set of pioneering researches, Nomikos and MacGregor (1995) proposed the batch-wise unfolding approach, which was the most suitable means for describing the differences among batches. In their approach, a 2-dimensional matrix X ∈ RI×JK was constructed out of the three-way batch process data X ∈ RI×J×K . The matrix X is formed by I batches, J manipulated variables, and K sample points. And the corresponding output data y ∈ RI is composed of I end-product qualities. After batch-wise unfolding, all the measurements obtained from a single batch are treated as discretized variables. Due to the correlationship in the input data matrix X and noisy, PLS algorithm is usually employed to build the data-driven model instead of OLS algorithm. In its basic form, nonlinear iterative partial least squares (NIPALS) method sequentially extracts the latent variables and weight vectors (Wold et al., 2001). Assuming the observations for each set of variables in X and y are scaled to zero mean and unit variance, the PLS ∼

algorithm decomposes the scaled X ∈ RI×JK and y˜ ∈ RI into the following equations: ˜ = TPT + E X

(13)

y˜ = Tq + f

(14)

where the superscript “∼” means scaled, T ∈ RI×A is the matrix composed of A latent variables, P ∈ RJK×A and q ∈ RA are the loadings, and E ∈ RI×JK and f ∈ RI are the residuals. The number of latent variable A is usually determined by cross-validation (CV), which select the one with the least square prediction error PRESSA associated (Xu and Liang, 2001). Based on the underlying model in Eqs. (13) and (14), we can rewrite the PLS model into a regression equation as follows: T

˜ yˆ PLS = y + y · (Bˆ x)

(16)

where W ∈ RJK×A is the matrix of the weight vectors.

where yˆ RSM is the predicted end-product quality of the response surface model, and ˇˆ i , ˇˆ ij and ˇˆ ii are the regression coefficients that can be obtained by ordinary least squares (OLS) algorithm.

2.3.

−1 Bˆ = W(PT W) q

(11)

Based on the results of DODE, a second-order response surface model can be built as follows:

yˆ RSM = ˇˆ 0 +

input data vector. The regression coefficient vector Bˆ ∈ RJK can be calculated using the following equation:

(15)

where yˆ PLS the predicted end-product quality of the PLS model, y and y are the mean and the standard deviation of the output data y respectively, and x˜ ∈ RJK is the new available scaled

3. Hierarchical batch-to-batch optimization method The hierarchical batch-to-batch optimization proposed in this work is a two-layer optimization architecture integrating the MA strategy based on response surface model and self-tuning batch-to-batch optimization method based on BW-PLS model, which enhances the optimization performances via combining their advantages. In this section, the architecture of the proposed hierarchical batch-to-batch optimization is firstly described, and then the two batch-to-batch optimization methods in the hierarchy are elaborated.

3.1. Architecture of hierarchical batch-to-batch optimization method The proposed hierarchical batch-to-batch optimization is a systematic method to determine the optimal control profile by two consecutive batch-to-batch optimization methods (MA strategy based on response surface model and self-tuning batch-to-batch optimization method based on BW-PLS model). MA strategy based on response surface model is applied at the upper level of hierarchy, aimed to determine the nominal control profile for the lower level of hierarchy. To this end, DODE method is firstly used to generate the dataset for building the response surface model. However, due to the plant-model mismatch, the optimal control profile obtained from the response surface model based optimization problem can be further improved by MA strategy. Taking the advantages of DODE method, the experiments required to generate the modeling dataset is quite moderate. At the same time, since a small number of DSFs is utilized and MA strategy directly employs the modifiers to compensate for the plant-model mismatch, only less iteration number is required to finish the procedure of batch-to-batch optimization at this level. Thus, the costs for the whole upper level of hierarchy are generally affordable, even for the batch process with high operational cost. The self-tuning batch-to-batch optimization method based on BW-PLS model is applied at the lower level of hierarchy. Since the self-tuning batch-to-batch optimization method at the lower level directly discretized the control profile into a large number of intervals, another dataset should be generated to calibrate the BW-PLS model. To avoid interrupting the regular production, the magnitude of the dither signal cannot be too large. Fortunately, only local BW-PLS model is required in self-tuning batch-to-batch optimization method, thus the costs of such dataset will not significantly deviate from that of the nominal control profile. Since the final control profile obtained from the upper level of hierarchy is selected as the nominal one at this level, all the costs of the modeling dataset should also be acceptable. On the other hand, although self-tuning batch-to-batch optimization method needs a large number of iterations to reach the optimum, the cost of batch process will be continually reduced. In this course, the regular production could not be interrupted. The flow chart of the proposed hierarchical batch-to-batch optimization method is

Chemical Engineering Research and Design 1 4 4 ( 2 0 1 9 ) 185–197

189

Fig. 2 – The flowchart of the data-driven model based hierarchical batch-to-batch optimization method. illustrated in Fig. 2, and the details will be elaborated in the following sections.

Hence, the response surface model based optimization problem can be formulated as follows:

Remark 1. The proposed hierarchical method makes the batch-to-batch optimization more practicable, especially for some batch processes with high operational cost. Since the MA strategy can efficiently resolve the plant-model mismatch, the DODE can also be performed in laboratory. Then, the optimal control profile of the upper level of hierarchy could be acquired after several iterations. In the light of the characteristics of the self-tuning batch-to-batch optimization method, some disturbances of the batch process can also be handled at the lower level. In most the time, the lower level of hierarchy operates the batch process. When some remarkable variability happens, the upper level of hierarchy can be redeployed to speed up the procedure of batch-to-batch optimization.

min ˚(˛) : = (˛, yˆ RSM )

3.2. Modifier-adaptation method based on response surface model In practical application, the mapping relating the control profile u(t) and the measured end-product quality yP is not known precisely. Thus the second-order response surface model illustrated in last section can be used to approximate the functional relationship, and the dynamic optimization problem can be transformed into a static one. In this work, the control profile is also parameterized as an expansion of the shifted Legendre polynomials. And assume that the following first three polynomials are used. ϕ0 () = 1,

ϕ1 () = −1 + 2,

ϕ2 () = 1 − 6 + 6 2

(17)

˛

s.t. G(˛) : = g(˛, yˆ RSM ) ≤ 0

(18)

− 1 ≤ ˛1 ± ˛2 ± ˛3 ≤ +1 where ˛ = [˛1 , ˛1 , ˛3 ]T is the vector of decision variables, and the predicted end-product quality yˆ RSM can be estimated using the response surface model in Eq. (12). Generally, the objective of OLS algorithm is only to minimize the mean square errors (MSE). Therefore even if the response surface model exactly predicts the end-product qualities, the gradients obtained from the model may not equal to the experimental ones. Then the operating point obtained from optimization problem (18) may not necessarily result in optimal plant operation (the control profile can be formulated by reusing Eq. (8)). In order to compensate for the plant-model mismatch, the forgoing modification procedure can be applied to determine the new operating point (Marchetti et al., 2009, 2010). Letting ˛(l) denote the lth operating point, the response surface model output can be modified as follows: yM (˛) = yˆ RSM (˛) + ε(l) + T (˛ − ˛(l) ) (l)

(19)

where ε(l) and (l) are the lth model modifiers, which can be calculated by using the following equations: ε(l) = yP (˛(l) ) − yˆ RSM (˛(l) )

(20)

190

(l) =

Chemical Engineering Research and Design 1 4 4 ( 2 0 1 9 ) 185–197

∂yP ∂yˆ RSM (˛ ) − (˛(l) ) ∂˛ (l) ∂˛

(21)

According to Eq. (12), the elements in ∂yˆ /∂˛ can be computed as follows: ∂yˆ RSM = ˇˆ i + ∂˛i

N 

ˇˆ ij ˛j + 2ˇˆ ii ˛i

(22)

j=1,j = / i T

and ∂yˆ RSM /∂˛ = [∂yˆ /∂˛1 , ∂yˆ /∂˛2 , ∂yˆ /∂˛3 ] . Hence, the next optimal operating point ˛∗(l+1) is calculated by solving the modified optimization problem in Eq. (23). min ˚M (˛) : = (˛, yˆ M ) ˛

s.t. GM (˛) : = g(˛, yˆ M ) ≤ 0

(23)

− 1 ≤ ˛1 ± ˛2 ± ˛3 ≤ +1 To avoid excessive correction, the next operating point ˛(l+1) can be filtered with a first-order exponent filter, i.e. ˛(l+1) = (I − K)˛(l) + ˛∗(l+1)

(24)

where K is a diagonal gain matrix. Example 1. A simulated batch reactor is first considered in this work. The reactions are A → B → C, which can be described by the following equations: 3 /348u

x˙ 1 = −4 × 103 e−2.5×10

3 /348u

x˙ 2 = 4 × 103 e−2.5×10

x12

(25) 3 /348u

x12 − 6.2 × 105 e−5×10

x2

(26)

where x1 and x2 are the dimensionless concentrations of the species A and B, and u is the dimensionless control profile. Suppose the initial conditions are x1 (0) = 1 and x2 (0) = 0, and the batch terminal time tf = 1. Then the objective of such batch process is to maximum the final concentration of the species B, i.e. x2 (tf ). And the hard constraint on control profile is 0.85 ≤ u ≤ 1.15. First, DODE is used to generate the dataset for building the response surface model. u0 () and u() are selected as 1 and 0.15, respectively. Since three DSFs are used in this example, a D-optimal with the point exchange algorithm defines the 13 experiments, which includes the 10 samples for response surface model building and 3 additional experiments to estimate the lack-of-fit (LOF) statistic. After the identification of the model parameters, the optimization problem (18) is directly solved. And the final concentration of species B has reached to 0.6082. Based on this result, the MA strategy is used to further improve the end-product quality, and the Broyden’s formula (Gao and Engell, 2005) is used to estimate the experimental gradients. After 6 iterations, x2 (tf ) ends at 0.6089 as shown in Fig. 3 (the batch number 0 means the best quality in the training dataset). Since only three DSFs are used in optimization problem (18), the end-product quality has slightly improved from 0.6082 to 0.6089. However, benefit from the fast convergence rate of the MA strategy, such growth takes merely 6 iterations. Remark 2. The quality of the solution is strongly dependent on the parameterization of the control profile. Thus, to further improve the batch-to-batch optimization performances,

Fig. 3 – Evolution of the cost for the proposed response surface model based MA strategy at the upper level of hierarchy in toy example. more DSFs (˛i ) can be used to define the control profile. However, as the number of DSFs grows, much more experiments should be performed to generate the dataset for response surface model training. For some industrial batch processes with high operational cost, such a large number of experiments are generally forbidden, even if these experiments are carried out in the laboratory (e.g. the cobalt oxalate synthesis process in this work). Therefore, although a suboptimal control profile is firstly implemented to operate the batch process, the optimization performances can gradually be improved by the following batch-to-batch optimization method in the lower level of hierarchy. Remark 3. Recently, we have directly discretized the control profile into moderate number of intervals, and performed MA strategy based batch-to-batch optimization (Jia et al., 2015). However, as the increasement of the intervals, such strategy may fail to find the optimal control profile due to the incapability of estimating the experimental gradients. Different from the batch-to-batch optimization method mentioned above, DODE method parameterizes the control profile only into several DSFs. Thus, the issue of experimental gradients estimation is more easily in this work. In addition, the response surface model used here is a global data-driven model. Therefore, it is unnecessary to update the response surface model after each batch run.

3.3. Self-tuning batch-to-batch optimization method based on BW-PLS model Since relative small number of DSFs is used to parameterize the control profile, the quality of the optimal solution can be further enhanced. Recently, a batch-to-batch optimization method is proposed by Camacho et al. (2007, 2015), and the BW-PLS model is used to capture the relationship between the discretized control profile u and the end-product quality yP . Since the gradients obtained from the BW-PLS model are only valid for the average control profile, a dither signal is appended to the batch process to guarantee the persistent excitation in their work. At the kth iteration, the control profile can be formulated as follows:

 ¯(k) + c(k) · u(k+1) = u

dˆyPLS du

 + g(k) · (k+1) (k)

(27)

Chemical Engineering Research and Design 1 4 4 ( 2 0 1 9 ) 185–197

191

¯ is the scaled average control profile that is calculated where u by the dataset for building the local BW-PLS model. dˆyPLS /du is the gradients obtained from the BW-PLS model, and the parameter c is the step length that can be tuned to further improve the performances of the optimizer, which can be tuned by using the following adaptive rule (Camacho et al., 2007). c(k) = 1 −

 PRESS  A

PRESS0

(28) (k)

where PRESSA is the square prediction error of the CV with A latent variables, and PRESS0 is the square prediction error computed only for the average output.  is a vector of excitation values corresponding to the control profile, and g is the excitation gain, which can also be tuned as follows: g(k) =

 PRESS  A

PRESS0

(k)

·g

(29)

and g is the constant parameter. Since the gradients are calculated by the BW-PLS model in their work, much more intervals can be uses to discretize the control profile (e.g. more than 100 intervals). However, the local BW-PLS model is only assumed to be valid in a relative small region. To ensure the model only used in this region, the gradients are constrained to move three times the standard deviation with the following definition (Camacho et al., 2015). dˆyPLS = 3 ·  u(k)  Bˆ (k) / Bˆ (k) 2 du(k)

(30)

where  u(k) is the standard deviation of the discretized control profile. Although a tight control was used in hydrometallurgical plant to guarantee the final product quality, the initial condition variability was difficult to control. The source of this variability can usually be measured before each batch run. Under this situation, the BW-PLS model can be enhanced by incorporating this information along with the manipulated variables. With this enhancement, the variability in the performance due to the initial conditions can be distinguished, and the feeding policy can be better corrected. The augmented input data matrix is defined as follows Xa = [Z X]

(31)

where Z ∈ RI×M is the matrix of M initial conditions over I batches. And PLS algorithm is performed on Xa to obtain the regression coefficient vector. As in Camacho et al. (2007), the contribution of the initial conditions for batch (k + 1) can be calculated by using the following equation. T

z(k) )  z(k) ) · Bˆ Z(k) yˆ Z(k+1) = ((z(k+1) − ¯

(32)

where z(k+1) is the input data vector corresponding to initial z(k) and  z(k) are the average and the standard deviconditions, ¯ ation of initial conditions, and Bˆ Z(k) is the regression coefficient vector corresponding to initial condition matrix Z. Then Eq. (30) can be modified as follows

dˆyPLS = du(k)

(3 − yˆ Z(k+1) ) ·  u(k)  Bˆ u(k) / Bˆ u(k) 2 0

yˆ Z(k+1) ≤ 3

otherwise

(33)

Fig. 4 – Evolution of the cost for the BW-PLS model based self-tuning batch-to-batch optimization method at the lower level of hierarchy in toy example. where Bˆ u(k) is the regression coefficient vector corresponding to the manipulated variables. In this work, the final control profile obtained from the last section is selected as the nominal one, and a pseudo random binary signal (PRBS) is also added to generate the initial dataset for BW-PLS model building. Then the BW-PLS model based self-tuning batch-to-batch optimization algorithm can be summarized as follows (Camacho et al., 2007): Step 1. Operate the new batch k with the control profile u(k) . Step 2. Update the BW-PLS model by incorporating the control profile and the end-product quality obtained from the last batch. Step 3. Estimate the gradients by using the new BW-PLS model, and calculate the next control profile u(k+1) in Eq. (27). Step 4. Check the termination conditions. If it is not met, increase the index k → k + 1, and go back to Step 1. Example 2. Also consider the process model in Example 1. We employ the self-tuning batch-to-batch optimization method to further improve the optimization performances. The optimal control profile obtained from the last example is used as the nominal one in this example. The control profile is discretized into 100 intervals, and the number of intervals is much larger than other batch-to-batch optimization methods. In order to generate the dataset for building the initial BW-PLS model, a PRBS with a magnitude of ±0.01 is added to excite the process dynamic. The average quality of the 10 samples for model building is 0.6086, which is not much worse than the nominal end-product quality (0.6089). After 100 iterations, the end-product quality is ended at 0.6105 as shown in Fig. 4, and the batch number 0 is the nominal end-product quality obtained from the upper level of hierarchy. In Fig. 5, the optimal control profile obtained from the proposed hierarchical batch-to-batch optimization method is compared with the one acquired by using simultaneous approach (Srinivasan et al., 2003). The initial control profile obtained by solving optimization problem (18) is also drawn in Fig. 5. It can be seen that the contours of the control profiles for these three methods are quite similar, especially for the first two profiles. And the optimal end-product quality is only slightly increased to 0.6107 for simultaneous approach. Remark 4. In self-tuning batch-to-batch optimization method, the BW-PLS model is rebuilt after each batch run

192

Chemical Engineering Research and Design 1 4 4 ( 2 0 1 9 ) 185–197

4. Application to cobalt oxalate synthesis process In this section, a simulation study for batch-to-batch optimization of cobalt oxalate synthesis process is presented. Firstly, the optimization problem of synthesis process is formulated, and the optimization results of the upper level of hierarchy are illustrated. Then the lower level of hierarchy continues to improve the optimization performances, which is further compared with other batch-to-batch optimization methods. To make the comparison fairer, all the data are scaled to zero mean and unit variance before data-driven model building. All the algorithms are implemented in Matlab 2016b, and the PLS toolbox is used to scale the dataset and train the BW-PLS model. Fig. 5 – Comparison of the control profiles for three different methods.

by using the moving window strategy. This is because the gradients are only locally valid for the average control profile of the dataset of the batches employed to build the BW-PLS model. In order to describe the local region precisely, the magnitude of the dither signal should not be too large. Under this circumstance, the final costs of the dataset for BW-PLS model building will not significantly deviate from the nominal one, and such loss is generally affordable even for batch process with high operational cost. Moreover, the new control profile is constrained by Eq. (30) during the evolution of batch-to-batch optimization, thus the risk of the unqualified batches is greatly reduced. Remark 5. Although the final result is slightly improved by simultaneous approach, an accurate first principle model should be calibrated in advance. However, the proposed hierarchical batch-to-batch optimization method needs little or no prior knowledge about the batch process. On the other hand, different from MFO method, carefully designed experiments are performed to generate the dataset for data-driven model building instead of the trial and error process. Thus, the proposed method is more suitable for the batch process with high operational cost. Since the input and output data reflect the current state of the batch process, the practical optimal control profile can finally be achieved with the evolution of batch-to-batch optimization. Remark 6. Due to the problem of variable duration batches, the proposed optimization method cannot be directly used. Under this circumstance, the measurements of manipulated variables for BW-PLS model building should be equalized. Work has been carried out into various methods for equalizing batch lengths (Kourti, 2003). Recently, Corbett and Mhaskar (2016, 2017) employed subspace identification methods and data-driven model to cope with the batch data with varying durations. However, due to the limitations of the BW-PLS model used in this work, the desired batch duration should be predetermined for batch-to-batch optimization. Then dynamic time warping (DTW) can be used to equalize the measurements of manipulated variables (Lu et al., 2004). Since the prediction accuracy of BW-PLS model is affected by the equalization, the proposed method would take more iteration to reach the optimal solution.

4.1. Analysis of the optimization results of the upper level of hierarchy The first principle model of the cobalt oxalate synthesis process described in Section 2.1 is used to test the proposed hierarchical batch-to-batch optimization method. The mean crystal size for each batch is selected as the end-product quality, which can be calculated as follows:

Ln =

1 (tf ) 0 (tf )

(34)

In this work, the objective function is to maximize the mean crystal size defined in Eq. (34), and leading to the following dynamic optimization problem (35).

min Ln F(t)



tf

c0 F(t) ≥ qmin

s.t.

(35)

0

0 ≤ F(t) ≤ Fmax

where Fmax is the upper bound on the feed rate of ammonium oxalate, and qmin is the minimum quantity of ammonium oxalate that is required to be added into the crystallizer. The first constraint is to ensure that almost all the cobalt chloride should be reacted with ammonium oxalate, and the second one is the hard constraint on the manipulate variable. At the upper level of hierarchy, the initial feeding profile is selected as F(t) ≡ 0.002m3 /s to imitate manual operating, and 0.5% white noise is appended to the end-product qualities to simulate the measurement noise. In this study, the feed rate of ammonium oxalate is also parameterized as an expansion of the first three shifted Legendre polynomials. And the feed rate is selected to satisfy the first constraint in optimization problem (35). The shifted Legendre polynomials are chosen as they represent a complete set of orthogonal polynomial functions, and this allows for generation of the time-variant feeding profiles over the dimensionless time between 0 and 1. The expansion of the feeding profile function is:

F() = F0 + F

3  i=1

˛i ϕi ()

(36)

Chemical Engineering Research and Design 1 4 4 ( 2 0 1 9 ) 185–197

Fig. 6 – Design space defined by Eqs. (11) and (35) for the three DSFs case in cobalt oxalate synthesis process.

Fig. 7 – Feeding profiles of ammonium oxalate generated by using the ten design points from within the design space shown in Fig. 6. where F0 and F are all selected as f0 = 0.002m3 /s according to the practical situation. For the case of three DSFs, the dimensionless feeding policy is parameterized as follows:





F() = (1 + ˛1 − ˛2 + ˛3 ) + (2˛2 − 6˛3 ) + 6˛3  2 · f0

(37)

In order to satisfy the first constraint in optimization problem (35), the following inequality can be obtained. q ˛1 ≥ min − 1 co f0

(38)

Then the design space defined by Eqs. (11) and (38) is shown in Fig. 6, and we also utilize D-optimal design with the point exchange algorithm to generate the 13 experiments. The first 10 samples are used for building the response surface model, and the other ones are employed to estimate the LOF statistic. Fig. 7 shows the dimensionless feeding profiles for the experimental design points of model building, and all the profiles satisfy the constraints in optimization problem (35). Second-order response surface model is fitted to the final mean crystal size using the OLS algorithm. In order to find the optimal values of the DSFs, the optimization problem (18) is solved firstly. Optimal feeding profile is generated by reusing

193

Fig. 8 – Optimization performance of the proposed response surface model based MA strategy at the upper level of hierarchy in cobalt oxalate synthesis process.

Eq. (36), and the end-product quality is 2.29 ␮m, which is much better than that of the manual operating (1.71 ␮m). Due to the plant-model mismatch, the response surface model predicted the final mean crystal size 2.39 ␮m. Hence, the optimal feeding profile obtained by offline optimization may result in suboptimal operation. Based on the optimal control profile of DODE, the MA strategy proposed in this work is also used to further increase the final mean crystal size. Since only three DSFs are selected as decision variables, the convergence rate is very fast. Although the final mean crystal size slightly increases from 2.29 ␮m to 2.36 ␮m (see Fig. 8), the proposed response surface model based MA strategy takes only 9 iterations. Moreover, during the progress of batch-to-batch optimization, the end-product qualities do not deteriorate. It should be noted that, Broyden’s formula (Gao and Engell, 2005) is also used to estimate the experimental gradients. Although other gradient estimation methods can be used to exactly identify the experimental gradients, it lies beyond the scope of this work.

4.2. Analysis of the optimization performances of the lower level of hierarchy The aim of the proposed hierarchical batch-to-batch optimization at the lower level is to further maximize the final mean crystal size. To this end, the feeding profile is directly discretized into 110 intervals. After the discretization, the dynamic optimization problem is transformed into a seeking problem in a 110-dimensional space, and the nominal feeding profile is the final optimal one obtained from the upper level of hierarchy. In this work, a much more ambitious optimization problem for synthesis process is pursued than many other works (Zhang et al., 2012; Jia et al., 2015). The feeding profile is also restricted by the two constraints in optimization problem (35), which directly limits some particular profiles before they are implemented to the synthesis process. If the feed rate of a certain interval is negative or above the upper bound Fmax , then it is set to 0 or Fmax , respectively. If the quantity of ammonium oxalate added into the crystal-

194

Chemical Engineering Research and Design 1 4 4 ( 2 0 1 9 ) 185–197

Fig. 9 – Comparison of the evolutions of the final mean crystal size for the four different optimization methods in cobalt oxalate synthesis process. lizer is less than the minimum quantity qmin , the following optimization problem is firstly solved.

ı∗ = argminJ(ı) := ı

tf  c0 ıui(k+1) − qmin N N

2 (39)

i=1

where ui(k+1) is the ith element of u(k+1) , and ı > 1. After obtaining ı∗ , the corrected control profile is: u(k+1) = ı∗ u(k+1) 

(40)

The goal of application to synthesis process is to evaluate the optimization performances of the proposed hierarchical method. To this end, the following four batch-to-batch optimization methods with different combination are compared in this work. • MA: MA batch-to-batch optimization method based on BW-PLS model. The feeding profile is discretized into 11 intervals, and F(t) ≡ 0.002m3 /s is selected as the nominal feeding profile. • ST(A): Self-tuning batch-to-batch optimization method based on BW-PLS model. The feeding profile is discretized into 110 intervals, and F(t) ≡ 0.002m3 /s is also selected as the nominal feeding profile. • ST(B): Self-tuning batch-to-batch optimization method based on BW-PLS model. The feeding profile is also discretized into 110 intervals, and DODE method is used to generate the nominal feeding profile. • PH: The proposed hierarchical batch-to-batch optimization method. At the upper level of hierarchy, MA batch-to-batch optimization method based on response surface model is used to generate the nominal feeding profile for the next level. At the lower level of hierarchy, self-tuning batch-tobatch optimization method based on BW-PLS model is used to further refine the feeding profile, which is also discretized into 110 intervals. The simulation study is performed to compare the three different batch-to-batch optimization methods and the one proposed in this work. For this comparison, the evolution of final mean crystal size for each method is drawn in Fig. 9. Although the convergence rate of MA method is very fast,

Fig. 10 – The HSD interval of the analysis of variance of ACC (a) and FNL (b) for the four different optimization methods in cobalt oxalate synthesis process.

the final mean crystal size merely achieves at 2.27 ␮m. This is because the feeding profile is only discretized into 11 intervals, which seriously impacts the batch-to-batch optimization performances. Compared with MA method, the final mean crystal size of ST(A) method is improved (2.39 ␮m), while the convergence rate is not better than that of the MA method as shown in Fig. 9. It should be noted that both of the methods select the F(t) ≡ 0.002m3 /s as the nominal feeding profile. For ST(B) method, the DODE method is used to generate the initial profile, and the optimization performance is greatly improved (see Fig. 9). Finally, the proposed hierarchical method is applied, and both the convergence rate and final mean crystal size (2.87 ␮m) are further enhanced. This is because a better nominal control profile can greatly improve the performances of self-tuning batch-to-batch optimization method. To make the comparison more clearly, the optimization results are also evaluated by two performance indices (Camacho et al., 2015). The first one calculates the accumulated final mean crystal size during the 100 batches, namely ACC. And the second one measures the average of the final mean crystal sizes of the last 5 batches, namely FNL. ACC distinguishes the faster optimization method, while FNL evaluates the effectiveness of the optimal solution when each method runs to convergence. The honestly significant difference (HSD) tests are illustrated in Fig. 10 with the multi-compare command in MATLAB, and the HSD plots are an option for multi-comparisons of a set of tests. 95% confidence level is used in such plots, which indicates that the compared methods are significantly different if their intervals are disjoint. Both the results of FNL and ACC show significant different statistical differences in the mean when different batch-to-batch optimization methods are applied. The MA method and ST(A) method, in cobalt oxalate synthesis process, only yield a mean FNL equal to 2.26 ␮m and 2.38 ␮m respectively, and a mean ACC equal to 109.92 ␮m and 104.41 ␮m respectively. Although the ST(B) method obtains a mean FNL equal to 2.78 ␮m and a mean ACC equal to 126.47 ␮m, the PH method attains a mean FNL equal to 2.87 ␮m and a mean ACC equal to 129.90 ␮m. The conclusions from Fig. 10 are that the proposed hierarchical batch-to-batch optimization method has higher performances.

195

Chemical Engineering Research and Design 1 4 4 ( 2 0 1 9 ) 185–197

Table 1 – The comparison of optimization performances with different batch-to-batch optimization methods.

Initial solution (␮m) Final solution (␮m) Number of iterations Number of intervals

FPM

Hybrid model

JIT-BW-PLS model

PH (1)

PH (2)

2.22 2.36 7 22

2.23 2.31 16 22

2.15 2.34 133 22

2.27 2.35 64 22

2.29 2.87 91 110

Fig. 11 – Optimization performance of the proposed hierarchical batch-to-batch optimization method under variable duration and initial conditions in cobalt oxalate synthesis process. To further evaluate the ability of the proposed method to handle variable duration and initial conditions, a total of 30 batches, whose terminal times range from 600 s to 720 s, are simulated to build the BW-PLS model. And the initial conditions are simulated by adding 2% of variability to the concentration of cobalt chloride solution and ammonium oxalate solution. In the lower level of hierarchy, DTW is first used to equalize these uneven-length training data set. In Fig. 11, the optimization performances are draw for this scene, and the average performance increases from 2.29 ␮m to 2.44 ␮m. Although the prediction accuracy is affected by variable duration batches, with the updating of BW-PLS model, the plant-model mismatch can be gradually diminished by the self-tuning batch-to-batch optimization method in the lower level of hierarchy. It should be noted that, due to the limitation of BW-PLS model used in this work, the batch duration should be fixed during the evolution of batch-to-batch optimization, and the desired batch duration is also selected as 660 s in this example. Finally, the proposed hierarchical method is compared with several other batch-to-batch optimization methods used in batch crystallization processes to further evaluate the optimization performances. The first one is the first principle model (FPM) based batch-to-batch optimization method. Paengjuntuek et al. (2008) proposed the optimization method for the estimation of uncertain kinetic parameters, and the updated kinetic parameters are applied for determining an optimal operating policy for the next batch run. The second method is based on hybrid model, which is composed of a nominal FPM and an updated PLS model (Hermanto et al., 2011; Zhang et al., 2012). And the updated PLS model is used to compensate for the mismatch between the nominal FPM and practical batch crystallization process. The third one is a data-driven model

based derivative-free batch-to-batch optimization method (Jia et al., 2018). The just-in-time (JIT) modeling method and BW-PLS (JIT-BW-PLS) model are integrated into the trust-region framework to perform batch-to-batch optimization. The batch process model described in Section 2 is also used to simulate the practical cobalt oxalate synthesis process. In case of a mismatch in kinetic parameters, the 5% increase of KG and 5% decrease of KB are introduced in the batch process model. The feeding profile is divided into 22 intervals to further enhance the optimization performances. For hybrid model based method, the moving window strategy is adopted to update the PLS model, while KG and KB are unchanged during the evolution of batch-to-batch optimization. The number of latent variables retained in the PLS model is determined by CV. In JIT-BW-PLS model based method, the size of relevant dataset is equal to 10, and CV is also employed to identify the number of latent variables. The optimization performances of these different batchto-batch optimization methods are compared in Table 1. After calibrating the initial (nominal) models, the initial solutions of these methods are reordered. The worse method is JIT-BWPLS model based method, since the JIT modeling method only selects the samples around the nominal feeding policy to build the initial local PLS model. Due to the existence of compensated model in hybrid model based method, the initial solution is slightly better than that of FPM based method. Since DODE is used in the proposed method to generate the dataset for RSM, a much better initial solution can be obtained. Each batch-to-batch optimization method is stopped when five consecutive batches do not outperform the average performance plus a threshold, and the average mean crystal size of the five batches is utilized as the final solution. At the same time, the numbers of iterations for different methods are also kept in Table 1. The convergence rates of FPM based method and hybrid model based methods are very fast. This is because a nominal FPM is used in both of the optimization methods. However, only updating the PLS model cannot completely compensate for the plant-model mismatch. For JIT-BW-PLS model based method and the proposed method, much more iterations should be required to reach the convergence, since there is no prior knowledge in both of the methods. Compared with the JIT-BW-PLS model based method, the proposed PH (1) method takes less iteration and obtains similar final mean crystal size due to the hierarchy structure. Moreover, since the proposed method PH (2) can also discretize the control profile into a large number of intervals (110), the final mean crystal size can be further improved, which is much better than that of other methods. It should be noted that, only two parameters are modified to introduce the plantmodel mismatch in this instance. In practical cobalt oxalate synthesis processes, much more parameters cannot be determined preciously, under this situation, the FPM based method may fail to improve the mean crystal size from batch to batch.

196

Chemical Engineering Research and Design 1 4 4 ( 2 0 1 9 ) 185–197

generate the dataset for response surface model building. Since the feeding profile was firstly parameterized into a small number of DSFs at the upper level of hierarchy, the MA strategy was proposed to directly compensate for the plant-model mismatch in this work. Then the self-tuning batch-to-batch optimization method was used to discretize the control profile into a large number of intervals, and the final cost of each batch can be further improved during the evolution of batch-to-batch optimization. The efficiency of the proposed hierarchical batch-to-batch optimization method was illustrated by using a toy example and a simulated cobalt oxalate synthesis process. Considering the advantages of the proposed hierarchical method, it is also attractive to extend the optimization method to other industrial batch processes with high operational costs. Fig. 12 – The results for optimization of shape of the particle size distribution by used the proposed hierarchical batch-to-batch optimization method.

4.3. Batch-to-batch optimization of shape of the particle size distribution Since the properties of cobalt oxalate are determined by the shape of PSD, a desired shape of the PSD may also be achieved by the proposed hierarchical batch-to-batch optimization method. Instead of using lumped values of the distribution such as moment, the population balance equations (PBE) together with the mass balances are employed in this work (Zhang et al., 2012). Then the objective of optimization is to obtain a desired distribution by the end of each batch while respecting the constraints in the original optimization problem. Following the idea in Garg and Mhaskar (2018), the optimization problem (35) can be reformulated. min J = (r − rdes )Q(r − rdes ) F(t)



tf

c0 F(t) ≥ qmin

s.t.

T

(41)

0

0 ≤ F(t) ≤ Fmax where r is the terminal distribution profile, rdes is the desired distribution profile, and Q is the weighting matrix. As shown in last example, the response surface model is firstly built, and the proposed MA strategy in the upper level of hierarchy is used to perform batch-to-batch optimization. And the matrix Q is selected as identity matrix. Then the lower level of hierarchy is employed to further refine the object of optimization. In Fig. 12, the optimization results are illustrated, where the initial shape is greatly different from the desired one. After 7 iterations, the shape is gradually improved by the upper level of hierarchy. After performing the self-tuning batch-to-batch optimization strategy in the lower level of hierarchy, the final shape of distribution is nearly the same as the desired one.

5.

Conclusions

This work presented operating the cobalt oxalate synthesis process in the architecture of hierarchical batch-to-batch optimization. The proposed optimization method was composed of MA strategy based on response surface model and selftuning batch-to-batch optimization method based on BW-PLS model. The DODE method was also used to systematically

Acknowledgements The authors would like to acknowledge the anonymous reviewers for their helpful comments. This work was supported by the National Natural Science Foundation of China (Nos. 61873049, 61733003, 61773105 and 61503384), the Fundamental Research Funds for the Central Universities (No. N160404004).

References Camacho, J., Picó, J., Ferrer, A., 2007. Self-tuning run to run optimization of fed-batch processes using unfold-PLS. AICHE J. 53, 1789–1804. Camacho, J., Lauri, D., Lennox, B., Escabias, M., Valderrama, M., 2015. Evaluation of smoothing techniques in the run to run optimization of fed-batch processes with u-PLS. J. Chemom. 29, 338–348. Corbett, B., Mhaskar, P., 2016. Subspace identification for data-driven modeling and quality control of batch processes. AIChE J. 62, 1581–1601. Corbett, B., Mhaskar, P., 2017. Data-driven modeling and quality control of variable duration batch process with discrete inputs. Ind. Eng. Chem. Res. 56, 6962–6980. Dong, D., McAvoy, T.J., Zafiriou, E., 1996. Batch-to-batch optimization using neural network models. Ind. Eng. Chem. Res. 35, 2269–2276. Fiordalis, A., Georgakis, C., 2013. Data-driven, using design of dynamic experiments, versus model-driven optimization of batch crystallization process. J. Process Control 23, 179–188. Flores-Cerrillo, J., MacGregor, J.F., 2003. Within-batch and batch-to-batch inferential adaptive control of semibatch reactors: a partial least squares approach. Ind. Eng. Chem. Res. 42, 3334–3344. Forgione, M., Mesbah, A., Bombois, X., Van den Hof, P.M.J., 2012. Iterative learning control of supersaturation in batch cooling crystallization. Proceedings of the American Control Conference (ACC), 6455–6460. Gao, W., Engell, S., 2005. Iterative set-point optimization of batch chromatography. Comput. Chem. Eng. 29, 1401–1409. Garg, A., Mhaskar, P., 2018. Utilizing big data for batch process modeling and control. Comput. Chem. Eng. 119, 228–236. Georgakis, C., 2013. Design of dynamic experiment: a data-driven methodology for the optimization of time-varying processes. Ind. Eng. Chem. Res. 52, 12369–12382. Hermanto, M.W., Braatz, R.D., Chiu, M., 2011. Integrated batch-to-batch and nonlinear model predictive control for polymorphic transformation in pharmaceutical crystallization. AICHE J. 54, 1008–1019. Jia, R., Mao, Z., Wang, F., He, D., 2015. Batch-to-batch optimization of cobalt oxalate synthesis process using modifier-adaptation strategy with latent variable model. Chemom. Intell. Lab. Syst. 140, 73–85.

Chemical Engineering Research and Design 1 4 4 ( 2 0 1 9 ) 185–197

Jia, R., Mao, Z., Wang, F., 2016. Self-correcting modifier-adaptation strategy for batch-to-batch optimization based on batch-wise unfolded PLS model. Can. J. Chem. Eng. 94, 1770–1782. Jia, R., Mao, Z., Wang, F., 2018. Combining just-in-time modelling and batch-wise unfolded PLS model for the derivative-free batch-to-batch optimization. Can. J. Chem. Eng. 96, 1156–1167. Klebanov, N., Georgakis, C., 2016. Dynamic response surface models: a data-driven approach for the analysis of time-varying process outputs. Ind. Eng. Chem. Res. 55, 4022–4034. Kong, X., Yang, Y., Chen, X., Shao, Z., Gao, F., 2011. Quality control via model-free optimization for a type of batch process with a short cycle time and low operational cost. Ind. Eng. Chem. Res. 50, 2994–3003. Kourti, T., 2003. Multivariate dynamic data modeling for analysis and statistical process control of batch processes, start-ups and grade transitions. J. Chemom. 17, 93–109. Lu, N., Gao, F., Yang, Y., Wang, F., 2004. PCA-based modeling and on-line monitoring strategy for uneven-length batch processes. Ind. Eng. Chem. Res. 43, 3343–3352. Mao, L., Shan, S., Yin, S., Liu, B., Wu, F., 1999. Effect of cobalt powder on the inner pressure of Ni/MH batteries. J. Alloys Compd. 293, 825–828. Marchetti, A., Chachuat, B., Bonvin, D., 2009. Modifier-adaption methodology for real-time optimization. Ind. Eng. Chem. Res. 48, 6022–6033. Marchetti, A., Chachuat, B., Bonvin, D., 2010. A dual modifier-adaptation approach for real-time optimization. J. Process Control 20, 1027–1037. Mesbah, A., Landlust, J., Huesman, A.E.M., Kramer, H.J.M., Jansens, P.J., Van den Hof, P.M.J., 2010. A model-based control framework for industrial batch crystallization processes. Chem. Eng. Res. Des. 88, 1223–1233. Nomikos, P., MacGregor, J.F., 1995. Multiway partial least squares in monitoring batch process. Chemom. Intell. Lab. Syst. 30, 97–108. Paengjuntuek, W., Kittisupakorn, P., Arpornwichanop, A., 2008. Batch-to-batch optimization of batch crystallization processes. Chin. J. Chem. Eng. 16, 26–29.

197

Roure, S., Bouvard, D., Dorémus, P., Pavier, E., 1999. Analysis of die compaction of tungsten carbide and cobalt power mixture. Powder Metall. 42, 164–170. Sanzida, N., Nagy, Z.K., 2013. Iterative learning control for the systematic design of supersaturation controlled batch cooling crystallization processes. Comput. Chem. Eng. 59, 111–121. Sarkar, D., Rohani, S., Jutan, A., 2006. Multi-objective optimization of seed batch crystallization processes. Chem. Eng. Sci. 61, 5282–5295. Shi, D., El-Farra, N.H., Li, M., Mhaskar, P., Christofides, P.D., 2006. Predictive control of particle size distribution in particulate processes. Chem. Eng. Sci. 61, 268–281. Srinivasan, B., Primus, C.J., Bonvin, D., Ricker, N.L., 2001. Run-to-run optimization via control of generalized constraints. Control Eng. Pract. 9, 911–919. Srinivasan, B., Palanki, S., Bonvin, D., 2003. Dynamic optimization of batch processes I. Characterization of the nominal solution. Comput. Chem. Eng. 27, 1–26. Wold, S., Sjötröm, M., Eriksson, L., 2001. PLS-regression: a basic tool of chemometrics. Chemom. Intell. Lab. Syst. 58, 109–130. Xu, Q.S., Liang, Y.Z., 2001. Monte Carlo cross validation. Chemom. Intell. Lab. Syst. 56, 1–11. Zhang, J., 2004. A reliable neural network model based optimal control strategy for a batch polymerization reactor. Ind. Eng. Chem. Res. 43, 1030–1038. Zhang, J., Nguyan, J., Xiong, J., Morris, J., 2009. Iterative learning control of crystallization process using batch wise updated linearised models identified using PLS. 19th European Symposium on Computer Aided Process Engineering (ESCAPE), 14–17. Zhang, S., Wang, F., He, D., Jia, R., 2012. Batch-to-batch control of particle size distribution in cobalt oxalate synthesis process based on hybrid model. Power Technol. 224, 253–259. Zhu, S., Yang, Y., Yang, B., Shao, Z., Chen, X., 2014. Model-free quality optimization strategy for a batch process with short cycle time and low operational cost. Ind. Eng. Chem. Res. 53, 16384–16396.