PDF (2 MB) - American Chemical Society

Apr 17, 2019 - and end constraints can be reformulated as maxβG(x, u, ti) ≤ 0 ..... plus and minus standard deviations from their mean values based...
1 downloads 0 Views 1MB Size
Subscriber access provided by UNIV OF LOUISIANA

Process Systems Engineering

Robust Batch-to-Batch Optimization with Scenario Adaptation Boeun Kim, Jakob Kjøbsted Huusom, and Jay Lee Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b06233 • Publication Date (Web): 17 Apr 2019 Downloaded from http://pubs.acs.org on April 28, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Robust Batch-to-Batch Optimization with Scenario Adaptation Boeun Kim1, Jakob K. Huusom2, Jay H. Lee1* 1. Department of Chemical and Biomolecular Engineering, Korea Advanced Institute of Science and Technology, 291, Daehak-ro, Yuseong-gu, Daejeon, 34141, Republic of Korea 2. Process and Systems Engineering Centre (PROSYS), Department of Chemical and Biochemical Engineering, Technical University of Denmark, Building 229, DK-2800 Kgs. Lyngby, Denmark

ACS Paragon Plus Environment

1

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 45

ABSTRACT

A robust batch-to-batch optimization (RB2BO) is proposed for recipe optimization of batch processes in the presence of uncertainty. Robust optimization can result in performance loss due to its conservative nature, and measurement-based optimization approaches do not typically address the residual uncertainty after parameter identification. In the proposed approach, distribution and extreme-case scenarios of uncertain parameters are refined through Bayesian updating using the measurements, and then the robust optimization is performed with the updated scenarios. Due to the iterative scenario adaptation and robust optimization mechanism, the proposed RB2BO approach is intended to provide robust recipes less conservative compared to the conventional methods, resulting in high objective function values and minimal constraint violations. This is demonstrated through an example of a pectin extraction process by considering feedstock variabilities and terminal constraints for the desired product quality. In addition, the termination criteria of the batch-to-batch correction are suggested, and the strategies for the parameter estimation and measurement are discussed to improve the performance of the adaptation.

KEYWORDS

Optimization

Under

Uncertainty;

Robust

Optimization;

Batch-to-batch

Optimization; Scenario Adaptation; Parameter Estimation; Optimal Experimental Design

ACS Paragon Plus Environment

2

Page 3 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

1. INTRODUCTION A batch process is widely used in productions of specialty chemicals, pharmaceuticals, food products, and other commodities due to its flexibility in operation to meet diverse product specifications1, 2. Major operational challenges in batch processes come from their inherent nonstationary characteristics involving the nonlinear dynamics of chemical and bio-based processes. Moreover, several operational constraints should be satisfied for ensuring productivity (or yield), safety, and product quality. Given the complexity of a batch operation, a model-based optimization can be a systematic way to improve process performance while satisfying the required constraints of operation and product quality3, 4. In the presence of uncertainty, the recipe resulting from an open-loop optimization with given nominal parameters may be highly suboptimal and result in constraint violations. There are various types of uncertainties such as model errors (structural or parametric), variations in feedstock, and process disturbances, resulting in batch-to-batch variations in product yield and quality, and a loss of process reproducibility4, 5. In practice, the kinetic parameters in the model can be changed with a scale-up, degradation of catalysts, changes in the operating condition, etc. Especially, in biological cultivation processes such as fermentation, model parameters can change after a restart with new seeds of the same strain due to the nature of living organisms6-8. One of the main sources of batch-to-batch variations is the feedstock variability, e.g., composition and quality, as inconsistencies in the initial condition directly affect the consistencies of the final product quality. For example, in a pectin extraction process from peels of citrus fruits, the quality of raw materials can vary considerably with the types of fruits9, seasons and geographical locations10, 11. Various approaches, e.g., modifier adaption and direct input adaptation12, have been suggested to handle the plant-model mismatch during the operation

ACS Paragon Plus Environment

3

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 45

by updating the gradients of the objective function and the constraints but it is difficult to estimate the gradients accurately in practice13. In this research, it is assumed that there is no critical model structural error, and model adequacy is guaranteed. The approaches to take parametric uncertainty into account in dynamic optimization of batch processes can be classified by the type of available information (i.e., measurement availability)5. To ensure the constraint satisfaction in the absence of measurements, a robust optimization (RO) is typically employed as a conservative approach to cope with uncertainty based on prior knowledge about the uncertainty due to its relative simplicity and computational tractability14,

15.

Typical RO does not require accurate probability distribution of uncertain

parameters and finds a feasible solution over a specified uncertainty set in the worst-case sense5. However, since the occurrence of the worst-case scenario is likely to have a low probability, the resulting solution can be very conservative resulting in a high performance loss. This conservativeness problem becomes severer when the worst-case scenario used in the RO is exaggerated for safety. When measurements are available, a measurement-based optimization (MBO) approach is implemented as a two-step approach to provide a less conservative operation5. In the MBO, uncertain parameters are re-estimated with new measurements, and recipe optimization is repeated using the updated parameters in a deterministic manner (by so called the ‘certainly equivalence’ principle) as shown in Figure 1A. However, during the adaptation period, the updated parameters can still contain significant errors due to the ill-conditioning and noise16. Moreover, in the presence of feedstock variability (or other time-varying uncertainty), the initial condition can vary batch-to-batch so the deterministic optimization using the estimation result of a current batch may provide either infeasible or sub-optimal recipe for a subsequent batch17.

ACS Paragon Plus Environment

4

Page 5 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Thus, the MBO fails to account the residual uncertainties. Methods for MBO can also be classified depending on the type of measurement available (on-line/off-line) and correction strategy (batch-to-batch/mid-course)4, 5. Since on-line measurements and corrections are still rare in practice, only the off-line measurement-based optimization is considered in this manuscript. In this paper, a new robust batch-to-batch optimization (RB2BO) framework is proposed to address the above-mentioned limitations of the two conventional approaches: the conservativeness in the RO and the residual uncertainty in the MBO. In the proposed framework illustrated in Figure 1B, distribution of uncertain parameters is updated using the estimation result with new measurements, and then scenarios for the RO are adjusted based on the updated distributions. For the scenario adaptation, Bayesian estimation is employed to update the distribution of uncertain parameters. This update on information about uncertainty can significantly reduce the conservativeness of the solution from the RO. In the decision making step, the RO with the updated scenarios, instead of the nominal model-based deterministic optimization, is implemented to account for the residual uncertainty after the parameter identification. In addition, an optimal sampling point design based on the results of Bayesian updating and recipe optimization is proposed to increase the information content in the available experimental data at next batch run but not to affect the process performance18. To reduce time and cost for measurement and optimization, termination criteria of batch-to-batch adaptation are also suggested in this study. Performances of the proposed RB2BO with and without using the termination criteria are evaluated through an example of a pectin extraction process from peels of citrus fruit wherein two cases of raw materials are considered to compare the proposed RB2BO and the conventional RO and MBO (in terms of the averaged objective function values and

ACS Paragon Plus Environment

5

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 45

degrees of constraint violation). Strategies for improving the performance of the proposed framework will be suggested based on the simulation results.

2. PRELIMINARIES Consider the process of which dynamics is described by the following state-space system: x& t   f  x  t  ,u  t  ,β 

(2)

y t   h  x t 

where x, y u, and β are the vectors of state, output, input and uncertain parameters, respectively. Functions f and h are assumed to be twice continuously differentiable. The state is initialized with x(t0)=x0. The initial condition is also considered as uncertain parameters, especially when the quality of raw material varies from batch to batch, in which case it could be included in β. A standard batch optimization problem can be stated as below:



u  arg max or min J   x  u, t f u



s.t. x& t   f  x  t  ,u  t   , x  0   x 0 y t   h  x t  G  x, u, t   0 t  0, t f 



T x t f

(3)

  0

ul  u  uu

where J is the scalar performance index to be maximized or minimized, and ϕ is a function of the terminal state defining the scalar performance index of the target batch operation. G(∙) and T(∙) represent the vectors of path and end constraints of the state, respectively, and decision variables u are bounded by lower bounds ul and upper bounds uu. tf is the final time of a batch operation,

ACS Paragon Plus Environment

6

Page 7 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

which can be either fixed or left as a free variable for optimization. Without loss of generality, we will assume a minimization of J from here on. Two conventional approaches to take the parametric uncertainty in β into account in the optimization are robust optimization (RO) and measurement-based optimization (MBO). In the absence of measurements, to ensure the constraint satisfaction, a scenario-based RO searches for an optimal and feasible recipe given ω scenarios of uncertain parameters Ω=[β1,…,βω] chosen based on a specified uncertainty set so that constraint satisfaction is guaranteed for all possible scenarios. Either the expectation or the worst-case of the objective function over the set is commonly used5. The expected value of the scalar performance function is given by



E  y  u, βl , t f 

  

(4)

and the mean value of productivity or cost can be computed as a weighted sum of the objective function values evaluated at given scenarios βl  Ω. Since the uncertainty set Ω is bounded, the worst-case objective function value and constraints can also be considered. Based on them, the optimization problem can be modified as a min-max problem:



min max  y  u, βl , t f u

βl Ω



(5)

where the resulting solution minimizes the maximal effect of the uncertainty on the performance index. In addition, the path and end constraints can be reformulated as max β G  x, u, ti   0 and





max β T x  t f   0 , respectively. For certain well-behaved system types, the worst-case scenario lies on the extreme points (e.g., those on the boundary or a vertex) of a given uncertainty set5. For general problems, a rigorous determination of the worst case scenario for a given recipe is very difficult, let alone optimizing it. For the ease of discussion, we adopt a discrete scenariobased formulation in this paper.

ACS Paragon Plus Environment

7

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 45

In the MBO, uncertain parameters are refined using measurements to reduce the effect of uncertainty and then the decision variables u are re-optimized with the updated parameters5. This two-step procedure is repeatedly executed with the advent of new measurements by the following two optimization problems, i.e., parameter estimation and recipe optimization. The MBO in the presence of uncertainty and measurements can be expressed as below:



max or min J   x k  u k , β k , t kf  k k k



u  tc , t f   





s.t. x&k  t   f  x k  t  ,u k  t  ,β k  , x k  0   x 0k y k  t   h  xk  t   G  x , u , β , t   0 t  t , t  k

k

k



k c

(6)

k f



T x k  β k , t kf   0 ul  u k  t   uu

where superscript k means a variable for the kth batch run and βk is the vector of uncertain parameters adjusted prior to the kth adaptation of the recipe. When the MBO is performed in a batch-to-batch manner, tck becomes zero and the optimal recipe for the entire kth batch uk 0, t kf  is obtained by solving eq (6). If the recipe optimization with mid-course correction is conducted, βk is identified using the data obtained up to the correction time tck and the recipe for the remaining batch time tck , t kf  of the kth batch is optimized. In the MBO, using the collected data, a typical way to update the uncertain parameters is by solving the following optimization problem, which is to minimize the weighted squared sum

ˆ and Ym, of the errors between the predicted and measured outputs corresponding to Y respectively:

ACS Paragon Plus Environment

8

Page 9 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

βˆ  arg min   u, β    u, β  



1 ˆ Y  u, β   Ym 2



T



ˆ  u, β   Y Cy1 Y m



(7)

where Ym ∈ ℝNm×Ny is the vector of all measured variables with Nm sampling instances. The measured outputs with additive noise ( y mj ∈ℝNm and j=1,…,Ny) are all collected in the vector of

 

m Ym: y mj  y j  ε j and Ym=  y1 

T

T

T ,K ,  y mNy   and the measurement noise εj is assumed to be 

independent and normally distributed with zero mean and constant variance σj2. The weighting matrix Cy1 is the inverse of the output experimental error covariance matrix Cy, which is a diagonal matrix containing the respective noise variances in such a case. If the estimation problem is well-conditioned, the parameter estimate βˆ will be uniquely determined with reasonably small variances. However, estimating all uncertain parameters often leads to an ill-conditioning problem resulting in large variances of the estimates19. For nonlinear models, according to the Cramér-Rao bound and the Gaussian approximation, the parameter variance-covariance matrix C can be approximated as below20:

C  FIM 1   ST Cy1S   VΖ 2 V T 1

(8)

where C is proportional to the inverse of the Fisher information matrix (FIM, or Hessian matrix). S is the sensitivity matrix obtained from partial derivatives of the outputs with respect to the parameters, and V is the matrix containing the orthogonal right singular vectors of the scaled sensitivity matrix S% Cy1/2S  UΖV T corresponding to the eigenvectors of C. Ζ is the diagonal matrix containing the singular values of S%, and the eigenvalues of C are the inverse of the squares of the singular values of S%. Therefore, ill-conditioning of an estimation problem can be

ACS Paragon Plus Environment

9

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 45

assessed by the singular value decomposition (SVD) of the sensitivity matrix16, 19. When S is a nearly singular (or singular) matrix with one or more nearly zero singular values, the estimation problem is ill-conditioned giving large variances (i.e., large eigenvalues of C) of the parameter estimates. Also, a high level of measurement noise leads to a large variance. To regularize the ill-conditioned estimation problems in practice, parameter subset selection (PSS) has been widely used due to its simplicity and interpretability21. In this approach, more influential and less correlated parameters are selected for estimation, while the rest of the parameters are fixed at their nominal values. For example, the orthogonalization method ranks the parameters based on the norm and linear dependency of each sensitivity vector22. In each iteration, it selects the parameter giving the largest norm of sensitivity vector and then subtracts the remaining parameters’ effects on the outputs covered by the previously selected parameter using Gram-Schmidt orthogonalization. Various methods for deciding a well-conditioned and meaningful parameter subset have been suggested16, 21, and a method for estimating a subset of transformed parameters has recently been proposed to increase the accuracy of parameter estimates in such problems23, 24. In this method, the original parameters are transformed to the directions of the principal components of the parameter covariance matrix (i.e., α  V T β ), and a subset of the uncorrelated transformed parameters are estimated. It doesn’t require the other ranking procedures since the transformed parameters are ranked in descending order of the parameters’ variances. In addition, the unselected transformed parameters are fixed at their nominal values to retain the prior physical meanings of the parameters, so the constraints on the original parameters can still be applied to the transformed parameters as follow:

ACS Paragon Plus Environment

10

Page 11 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

βl  V2α 2   V1α1  βu  V2α 2 

(9)

where βl and βu are lower and upper bounds on the original parameters, respectively.

3. ROBUST BATCH-TO-BATCH OPTIMIZATION (RB2BO) The proposed robust batch-to-batch optimization (RB2BO) comprises an update on the information about uncertain parameters and scenarios using new information that have available prior to a batch run (e.g., the result of a previous batch run), and an update of the recipe based on the updated information. An adaptation of the uncertain parameters and the scenario set used for the robust optimization is suggested as a way to update the recipe using new measurements in Section 3.1. Based on the updated information, the model-based optimizations of the recipe and sampling strategy are explained in Section 3.2. The overall algorithm of the proposed RB2BO is described in Section 3.3.

3.1 Update of distribution and scenarios Since uncertain parameters can seldom be directly measured, they should be estimated using measured outputs. For this, the parameter estimation problem can be formulated as in eq (7). In the RB2BO, based on the estimation results D=( βˆ 1 ,…, βˆ k ), the probability distribution of the uncertain parameters is recomputed to adjust the scenarios for the RO. The information about uncertainty can be conveniently and recursively updated in the Bayesian framework where the posterior distribution p(β|D) is proportional to the multiplication of the prior distribution p(β) and likelihood p(D|β). The prior distribution represents the existing knowledge before the measurements and D is considered as the new observed data that can be used to update the

ACS Paragon Plus Environment

11

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 45

distribution. To express the posterior distribution in a closed form, a conjugate prior is generally used if the posterior distribution has the same shape as the prior probability distribution. The corresponding conjugate priors for various distribution types have been investigated25. For example, Bayesian inference for a mixture of Gaussians can be found in 26, 27. In this paper, we provide the specifics for the Bayesian updating for uncertain parameters which are represented by a multivariate normal distribution with constant mean μ ∈ ℝn and covariance matrix Σ∈ℝn×n. Σ is the symmetric positive-definite matrix and n is the number of the uncertain parameters. To define proper scenarios for the RO, the covariance matrix as well as the mean values of the uncertain parameters need to be updated. This is especially true in the presence of batch-to-batch variations due to feedstock variability. The normal-inverse-Wishart (NIW) distribution p(μ,Σ)∝NIW(μ,Σ|μ0,κ0,Σ0,ν0) is used as the conjugate prior of a multivariate normal distribution with unknown mean and covariance matrix25, 26. Using the NIW prior, the closed form of the posterior p(μ,Σ|D,μ0,κ0,Σ0,ν0) after the kth batch operation is given by μk 

κ0 k μ0  β κ0  k κ0  k

Λk  Λ0  P 

T κ0 k β  μ 0  β  μ 0   κ0  k



P   j 1 β  βˆ j k

(10)

 β  βˆ 

(11)

T

j

where κ0 and ν0 (>n-1) are the degrees of freedom indicating the strengths of the prior mean and covariance matrix, respectively. k is the batch index, κk=κ0+k, and νk=ν0+k. If κ0 and ν0 are large values, the prior mean and covariance matrix strongly affect the posterior distribution. P is the scatter matrix, i.e., a sample sum of squares matrix, and β is the mean of the estimates

1 k ˆ  βj k j 1

or observed data. In eq (10), the posterior mean μk is calculated from the convex combination of

ACS Paragon Plus Environment

12

Page 13 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

the prior mean μ0 and maximum likelihood estimator (MLE), and it is adjusted towards the MLE as k increases. Hyperparameter Λ0 is defined as the ν0Σ0, and the posterior hyperparameter Λk is obtained by the contributions from the prior Λ0, the observations P, and the uncertainty in β as presented in eq (11). The posterior mode of the covariance matrix Σk can be derived from Λk divided by (ν0+k+n+2)26. In the case study in Section 4, κ0 and ν0 are chosen as 1 and 4, respectively, to accelerate the convergences of μ and Σ to their MLE. In the presence of uncertainty in model parameters, the same closed form as in eqs (10) and (11) can also be employed for updating the information about the uncertainty. When the uncertain model parameters are time-invariant, say, multivariate normal with known covariance matrix Σ corresponding to the parameter covariance matrix, the information is used as the conjugate prior and thus the closed form of the posterior after the kth batch operation is given by

Σ k   Σ 01  kΣ 1  and 1

μ k  Σ k  Σ 01μ 0  kΣ 1β  24,25. To address model structural error,

augmenting state and/or adding state noise can be considered in the proposed RB2BO context similar to the state augmentation used in the state estimators (e.g., moving horizon estimator28, extended Kalman filter29, etc.) to cover the plant-model mismatch and/or unmodeled disturbances. However, to be rigorous, significant a priori knowledge about structural mismatch is required, which is most likely to be not available, and the estimates can become biased if unanticipated mismatch occurs 30. After the Bayesian updating, a scenario set Ωk is adapted to determine a less conservative recipe in the RO. The scenario selection for RO highly depends on the aim of RO and the target system’s behavior with respect to a given uncertainty set, e.g., the manner in which the set affects

ACS Paragon Plus Environment

13

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 45

the feasibility or objective function value of the optimization31. Since the choice of the uncertainty set directly affects the computational tractability of the resulting optimization as well as the performance of the process, it must be carefully designed for each case. As mentioned before, for certain well-behaved problem types (e.g., convex ones), the worst-case scenario lies at the extreme points of an uncertainty set, and therefore a scenario set for RO can be chosen to include such points5, e.g., plus and minus standard deviations from their mean values based on prior knowledge or historical data about uncertainty (i.e., μk±σk, μk±2σk or μk±3σk depending on the level of conservativeness desired). For the classical norm-based bounded uncertainty sets (i.e., box, ellipsoidal and polyhedral uncertainty sets), those points on the boundary or vertex can be searched using given adjustable parameters and nominal values of uncertain parameters. The size of the uncertainty set depends on the selected norm and the corresponding adjustable parameter, affecting the conservativeness and robustness of a solution from the RO32. On the other hand, it is very difficult to rigorously guarantee an inclusion of the worst-case scenarios in the scenario set in the general nonlinear case5. It should also be combined with a method for uncertainty propagation, e.g., Monte Carlo simulations, polynomial chaos expansions33, unscented transformation34, etc. Despite the difficulty, in many practical problems, certain physical considerations may lead to a small finite set of parameter values that can be checked for a worst-case performance or constraint violation and such a set may be adopted in the RO.

3.2 Decision making: recipe optimization and optimal sampling design There is a conflict between objectives of exploration and exploitation in experimental design. Maximizing a performance index (exploitation) may not be equivalent to or even conflict

ACS Paragon Plus Environment

14

Page 15 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

with maximizing the information content in collected data for the parameter identification (exploration)18. Operating variables (e.g., temperature and feeding rate) may affect the process performance and the information content of data gathered, whereas decisions related to measurement (e.g., measured variables, type, and location of sensors and sampling schedules) may influence only the information content. Therefore, in the proposed RB2BO, sequential optimization problems are formulated to determine the major operating variables without regard to the information content first and then to determine the sampling points to maximize the information content, i.e., the exploration is subordinated to the exploitation, since the optimization of process performance is the most important goal in model-based operation optimization and considering both simultaneously in an optimally balanced way is not trivial 18. Since the reduction in process uncertainty affects the control performance35, approaches to solve the dual control problem can also be applied to RB2BO36. However, there are still significant challenges in optimally designing the probing signal to balance the objectives of learning and control. Firstly, with the updated scenario set Ωk, the following robust optimization is solved to decide the recipe for the next batch operation. The optimization problem in eq (6) can be rewritten as below:









u k 1 0, t kf 1   arg max or min J   wll xlk  u, βlk , t f   



u 0,t f 



l 1



s.t. x&lk  t   f  xlk  t  ,u  t  ,βlk  , x k  0   x 0k y lk  t   h  xlk  t   G  xlk , u, βlk , t   0 t  0, t f 



T xlk  βlk , t f

   0 β

k l

(12)

 Ωk

ul  u  t   uu

ACS Paragon Plus Environment

15

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 45

where wl is the weight coefficient for the lth scenario. These weight coefficients are normalized to sum up to 1. J is expressed as a weighted sum of the objective function values calculated for the given scenarios βlk  Ωk. Secondly, sampling points during the k+1th batch operation are optimized with uk+1 to maximize the information content in the available experimental data. A model-based optimal design of the sampling points can be formulated as the following optimization problem: st k 1  arg min ΨOED  Ck 1 

(13)

where stk+1 is the decision variables giving the minimum objective function value of ΨOED. Since the confidence region of the parameter estimates can be evaluated using the parameter variancecovariance matrix C, ΨOED is a scalar function of Ck+1 calculated by using eq (8) with μk and uk+1. Among various optimality criteria, the D-criterion ΨD=det(C)1/p is the most widely used owing to its ease of geometrical interpretation and scaling invariance37. To alleviate the combinatorial problem of the optimal experimental design in eq (13), the sampling time interval can be fixed at a constant value18.

3.2 Overall algorithm of the proposed RB2BO The overall algorithm of the proposed RB2BO is illustrated in Figure 2, and proceeds as the following steps: Step 0 – Initialize: Start with μi,0, Σi,0, and Ωi,0 and k = 1. Step 1 – Robust optimization: Perform the robust optimization in eq (12) to find the optimal and feasible recipe for given Ωi,k-1. Step 2 – Optimal sampling point design: Search the optimal sampling points by solving eq (13).

ACS Paragon Plus Environment

16

Page 17 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(the kth batch operation and measurement collection are executed using ui,k and sti,k) Step 3 – Parameter estimation: Identify the uncertain parameters using Ymi ,k . Step 4 – Bayesian updating: Update the distribution of the uncertain parameters through the Bayesian update as represented in eqs (10) and (11). Step 5 – Scenario adaptation: Update the scenario set based on μi,k and Σi,k. Step 6 – Termination of the inner loop: Return to Step 1 with k=k+1 until the following criteria are satisfied (if k>5):

 ij,q 1   ij,q  ò, q  k , , k  5and j  1,..., n  ij,q 1 Ωij,q 1  Ωij,q Ωij,q 1

 ò, q  k , , k  5and j  1,..., 

(14)

(15)

Step 7: Return to Step 0 and reset the prior information about uncertainty with i=i+1. where i is the event index for the uncertainty propagation. Equations (14) and (15) represent the relative errors of two successive mean values of the uncertain parameter and extreme-case scenarios, respectively. The inner loop is continued until the relative errors for five batches in a row become smaller than a tolerance ( ò). Also, once a major change occurs in the uncertain parameters (e.g., the type of raw material is changed), reset the prior information about the uncertainty (i.e., the mean, covariance matrix, and scenarios) and then restart the robust batch-tobatch adaptation (i=i+1).

4. CASE STUDY 4.1 Process description – pectin extraction process

ACS Paragon Plus Environment

17

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 45

A batch pectin extraction process from citrus fruits such as lemon, orange, and lime is adopted as a case study for comparing performances of the conventional methods and the proposed RB2BO method. The extraction process is the first processing step in commercial pectin manufacturing and crucial to the yield and product quality attributes. To extract pectin (water-soluble) and protopectin (water-insoluble), which is a combination of pectin and cellulose, the protopectin is hydrolyzed first using acid and then the pectin is released from the peels by water38. The quality of the extracted pectin highly depends on its physical and chemical characteristics such as intrinsic viscosity (IV) and degree of esterification (%DE)11. A high temperature and a low pH yield a high pectin concentration but IV and %DE significantly decrease at these conditions, which promote the hydrolysis, de-polymerization, and deesterification. Therefore, these operating variables should be carefully decided to achieve a high pectin yield with desired quality of the extracted pectin. Due to the natural variability of raw material quality in this as in many bio-based processes, there are uncertainties in the initial conditions for even a same type of raw material leading to the variations in the quality of end product. In addition, the direct measure of the initial conditions is impractical before each batch run, so they should be estimated based on measured data. To consider the inherent uncertainty in the feedstock and satisfy the required quality of the extracted pectin, Caroço et al. (2018) carried out an analysis of raw material variability in a pectin extraction process and performed a model-based robust optimization using prior knowledge about the feedstock uncertainty obtained from historical data9. In practice, the distribution of the initial condition varies by season, location of the supplied dried citrus peels, and pretreatment10, 11. Therefore, the uncertainty scenarios constructed from the historical data

ACS Paragon Plus Environment

18

Page 19 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

may be exaggerated or not include some possible scenarios for new peels, resulting in a performance loss or a constraint violation.

4.1.1 Dynamic model A dynamic model of the pectin extraction process developed by Andersen et al. (2017) is employed in the case study since this model can predict IV, %DE and pectin concentration with respect to varying initial conditions of raw material and major operating variables, i.e., temperature (Temp.) and pH. For the modeling, they assumed that (i) the peels have a same shape of flakes, characterized by a one-dimensional slab geometry, (ii) protopectin and pectin are uniformly distributed in the peel initially, and (iii) the effects of viscosity on diffusion and mass transfer can be neglected. A detailed description of the assumptions and model parameters can be found in Andersen et al. (2017)39 and Caroço et al. (2018)9.

In the pectin extraction process, the pectin concentrations in the peel and bulk solution, %DE, and IV are described by the following equations involving the acidic hydrolysis of the protopectin, diffusion of the pectin in the peel, mass transfer between the peel surface and the bulk solution, and degradation and de-esterification of the pectin in the bulk solution: peel C protopectin t, x 

t peel C pectin t, x 

t

C bulk pectin  t , L  t



 D pectin

Atotal kmasstransfer Vtotal

peel  khydrolysis C protopectin t, x 

peel  2C pectin t, x 

C

peel  khydrolysis C protopectin t, x 

(17)

bulk  t , L   C bulk pectin  t , L    k degradation C pectin  t , L 

(18)

x peel pectin

(16)

2

ACS Paragon Plus Environment

19

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

bulk Cdegraded pectin  t 

t bulk Cester t, L 

t

 DE 0

Atotal kmasstransfer Vtotal %DE  t  

C

peel pectin



Page 20 of 45

 kdegradationC bulk pectin  t 

(19)

bulk  t , L   C bulk pectin  t , L    k de  esterification Cester  t , L 

(20)

bulk Cester  t  100

bulk fGA C bulk pectin  t   Cdegraded pectin  t 

IV  t   IV0 exp  k IV t 



(21)

(22)

where C with superscript peel represents the concentration of protopectin or pectin in the peel, and C with superscript bulk indicates the concentration of pectin, degraded pectin or ester groups on the pectin in the bulk solution. Inside the peel, the insoluble protopectin is converted to pectin through the acid hydrolysis as described in eqs (16) and (17), and the internal diffusion of the pectin is depicted by Fick’s first law in eq (17). At the boundary of the peel (x=L), the mass transfer of the pectin occurs due to the difference between the concentrations of pectin at the peel surface and in the bulk solution, as shown in eqs (18) and (20). The reactions of degradation and de-esterification of the pectin, and the change in IV are described as first-order reactions as presented in eqs (18) – (20) and (22). All reaction rate constants for the hydrolysis (khydorlysis), degradation (kdegradation), de-esterification (kde-esterification) and change in the IV (kIV), and the diffusion coefficient (Dpectin) are expressed using modified Arrhenius type functions since they depend on the process temperature and/or pH9, 39. Model parameter values are summarized in Table 1.

4.1.2 Optimization problem

ACS Paragon Plus Environment

20

Page 21 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

In a model-based optimization of the pectin extraction process, temperature and pH are optimized to maximize the final pectin concentration in the bulk solution, and the desired quality specifications of the extracted pectin are enforced by inequality constraints on the IV and %DE at the final time (tf). Thus, the conventional robust optimization problem can be stated as below: max J  u

1





C l 1

bulk pectin

 u, β , t  l

f

s.t. x&l  t   f  xl  t  ,u, q  , x  0   βl



y l  t   h  xl  t  

T xl  βl , t f

  0

(23)

βl  Ω

ul  u  uu peel peel where the uncertain parameters β are the initial values of IV, %DE, C pectin and C protopectin . The

vector of final states xl(tf) is derived from the simulation of the nonlinear process model described in Section 4.1.1 initialized with given βl, and the mean value of the final pectin concentration for all βl (l=1,…,ω) is maximized as represented in eq (23). There are no path constraints on the state variables. Other design variables q are fixed a priori based on the standard conditions used in commercial-scale pectin production: the amounts of peel and water are set as 900 and 150,000 kg with the extraction time tf of 9 hrs. In addition, the constraints on the final IV and %DE are given according to product applications. In this case study, the required characteristics of the pectin for the jelly production (termed “jelly case” hereafter) are incorporated into the optimization: IV(tf) ≥ 6 dl/g and 58 ≤ %DE(tf) ≤ 659.

4.2 Comparison of performances of the conventional methods and the proposed RB2BO method In this section, performances of the conventional RO and MBO methods and the proposed RB2BO method are compared in terms of the averaged objective function value and

ACS Paragon Plus Environment

21

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 45

degree of constraint violation during 100 batches using two cases of peels. The number of batches (k) is set as 100 since about 100 batch extractions are executed for a single type of peel in actual industrial practice. Compared to the historical data of the lime peels, Case 1 has a higher mean value of IV0 and a lower mean value of %DE0, whereas Case 2 has higher mean values of both IV0 and %DE0, and both cases show smaller variances of IV0 and %DE0 (see Table 2). In the simulations for each case, actual initial conditions are sampled from multivariate normal distributions with their mean values reported in Table 2. Since the final values of IV and %DE for a given operating condition are monotonically affected by their initial values, the worst-case scenarios were shown to lie on the boundary of the uncertainty set in this example9. Therefore, based on this prior knowledge about the process, the extreme points of the uncertainty set are included as scenarios for the conventional RO and the proposed RB2BO methods in this case study. For the conventional RO, the scenario set Ω is chosen to be the same (i.e., with μ±σ of the initial conditions from the historical data in Table 2) as in the previous study9. In step 5 of the RB2BO, μi,k(IV0)-2×σi,k(IV0) and μi,k(%DE0)±2×σi,k(%DE0) are used as the extreme-case scenarios for the k+1th robust optimization in consideration of the limits of the final IV and %DE in the jelly case. In the next subsections, the adaptation results of the two scenarios, i.e., μi,k(IV0)2×σi,k(IV0) and μi,k(%DE0)+2×σi,k(%DE0), will be represented and analyzed in the RB2BO as they represent the worst-case scenarios (within the corresponding level set) in terms of the end quality constraint satisfaction. For the optimal sampling point design, the D-criterion is used to determine sti,k in step 2 of the RB2BO, and the number of sampling points and sampling time interval are set as 5 and 1 hr, respectively. During the kth batch operation, IV, %DE, and pectin concentration are assumed

ACS Paragon Plus Environment

22

Page 23 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

to be sampled and analyzed at sti,k, and the error variances for them are fixed at 0.052, 0.52 and 0.12, respectively. However, according to the sensitivity matrix of the dynamic model, the final peel peel peel value of %DE is affected by the changes in C pectin ,0 as well as %DE0. C pectin ,0 and C protopectin ,0

inevitably influence the final pectin concentration. Therefore, the estimation of the initial conditions of the pectin extraction process is an ill-conditioned estimation problem. To increase the accuracy of the estimates (especially those of %DE0), a method estimating a subset of transformed parameters23 described in Section 2 is employed. Thus, three of four transformed parameters are estimated with Ymi ,k in step 3 of the RB2BO, and then the distribution of the initial conditions are updated using eqs (10) and (11). Note that, to improve the estimates of peel peel C pectin ,0 and C protopectin ,0 , water extraction experiment needs to be executed before the acid

extraction but this is thought to be impractical. Fortunately, in this case study, precise estimates of IV0 and %DE0 are required given the constraints on their final values, but the variations in peel peel C pectin ,0 and C protopectin ,0 are not critical so just their mean values are used in the recipe

optimization. The nonlinear constrained optimizations for the parameter estimation and recipe optimization are solved using fmincon function in MATLAB R2017a Optimization Toolbox. The peel peel bounds on IV0, %DE0, C pectin ,0 and C protopectin ,0 for the estimation are given by 7 – 9 dL/g, 70 – 80

%, 50 – 150 kg/m3 and 20 – 230 kg/m3, respectively. In the simulation of the dynamic model given in Section 4.1.1, eq (17) is reformulated as ordinary differential equations (ODEs) through spatial discretization (into 10 points) using the central difference method and thus 15 ODEs are integrated by using ode15s in MATLAB R2017a as in the previously published modeling paper 39.

For the recipe optimization, all three methods determine optimal temperature and pH while

ACS Paragon Plus Environment

23

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 45

satisfying the two inequality constraints for the jelly case. The conventional RO and proposed RB2BO methods optimize them for the two given initial condition scenarios (ω=2). Average computational time for the conventional RO and MBO methods and the proposed RB2BO method are 113.1, 41.1 and 45.5 secs for this case.

4.2.1 Case 1 In Case 1, the uncertainty set based on prior knowledge is exaggerated; the true worstcase scenario of IV0 is higher than that from the historical data, though the true worst-case scenario of %DE0 is similar to that from the historical data. Consequently, the conventional RO may give a very conservative recipe leading to a high performance loss in this case. To demonstrate the performance of the proposed RB2BO, 100 batches of extractions are simulated using the recipes resulting from the MBO, RO, and the proposed RB2BO. In addition, the RB2BO is tested with and without using the termination criteria in eqs (14) and (15). Table 3 represents the simulation results of Case 1. As a result, 47 % of the 100 batches violate the constraint for the jelly case when the conventional MBO is employed. Using the recipe from the conventional RO, the constraints are violated only 1 % of the time; however, the averaged final pectin concentration decreases by about 10.6 %. On the other hand, the proposed RB2BO without using the termination criteria (i.e., the adaptation is performed after every batch) provides more robust handling of the constraints than the conventional MBO and also gives about 7.3 % higher averaged final pectin concentration than the conventional RO while showing a similar % of constraint violation as the conventional RO. In the proposed RB2BO, since the mean, variances and scenarios of the uncertain parameters are updated based on the measurements resulting from batch runs as shown in Figure 3, the conservativeness of the recipe

ACS Paragon Plus Environment

24

Page 25 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

decreases as the batch runs are repeated. In Case 1, when employing the proposed RB2BO, the scenario of IV0 gets closer to the actual scenario as k increases (Figure 3A) and this adaptation of the scenarios considered in the robust optimization contributes to the increase in the averaged final pectin concentration. The batch-to-batch correction can be stopped when the termination criteria in eqs (14) and (15) are satisfied. In Case 1, the adaptation is performed until the 39th batch run, and then temperature and pH become fixed for the rest of batches. As reported in Table 3, the proposed RB2BO with using the termination criteria provides almost identical performance to that without using the termination criteria. Thus, the time and cost for the sampling and analysis and solving the optimization problems can be saved by using the termination criteria.

4.2.2 Case 2 In Case 2, actual values of IV0 and %DE0 are higher than those seen from the historical data so the upper limit of the final %DE is violated when using the recipe obtained from the conventional RO; the prior worst-case scenario of IV0 is exaggerated and that of %DE0 is poor. The 100 batch runs of the pectin extraction process are simulated as in Case 1, and the performances of the four optimization strategies are summarized in Table 4. With the use of the termination criteria, the resulting recipe after the 35th adaptation is used without further change for the remaining batches in the proposed RB2BO. Using the recipe from the MBO, the required constraint for the jelly case is violated 89 % of the time. The conventional RO shows the violations of the constraint on final %DE 86 % of the time with 12.7 % reduction in the averaged final pectin concentration compared to the MBO.

ACS Paragon Plus Environment

25

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 45

As shown in Table 4 and Figure 4B, when employing the proposed RB2BO, the possibility of the undesired constraint violation significantly decreases to 4 % since the updated scenario of %DE0 approaches the true one (=76.13), which is higher than the prior assumed (=74.93) as k increases. Moreover, the extreme-case scenario of IV0 also increases from the exaggerated prior to its true value after a few batches (Figure 4A). Thus, with or without the termination criteria, the averaged final pectin concentration from the RB2BO increases by about 8.8 % compared to that from the conventional RO (Table 4). These results are consistent with those from Case 1.

4.3 Discussion on strategies for parameter estimation and measurement In both Case 1 and Case 2, the adaptation results for %DE0 show a little bias between the true mean value and updated mean value or the true worst-case scenario and updated extremecase scenario. This is because of the ill-conditioning with respect to the estimation of %DE0 and the relatively high noise level for measuring %DE. Accuracy of the parameter estimates is important to improve the fidelity of the extreme-case scenarios and ultimately the quality of the resulting recipe in the proposed RB2BO. Therefore, strategies for the parameter estimation and measurement (e.g., estimation method, decision and measured variables, sampling schedules, etc.) should be carefully designed, especially for ill-conditioned estimation problems. For example, to increase the accuracy of the parameter estimates in the case study, the estimation technique and the number of sampling points during one batch are decided through Monte Carlo simulations. Figure 5 shows the mean squared errors in IV0 and %DE0 resulting from 5,000 runs of Monte Carlo simulations using two different parameter subset selection (PSS) methods with nine different numbers of sampling points; three (or transformed) parameters are

ACS Paragon Plus Environment

26

Page 27 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

selected for estimation using the PSS by orthogonalization (PSS-O)22 and by transformation peel (PSS-T)23. Using the PSS-O, %DE0 is the top-ranked parameter followed by IV0, C pectin ,0 , and

peel C protopectin ,0 in this case study. For the simulations, true initial values are sampled assuming

multivariate normal distributions with their mean values obtained from the historical data reported in Table 2. The measurement noise variances, sampling time interval and batch time are the same as in the case study (Section 4.2.), and the additional sampling time is determined by eq (13) with the D-criterion. For all cases, the two methods have identical mean squared error in IV0, which is not correlated with other uncertain parameters. On the other hand, the PSS-T method gives about 51 % smaller mean squared error in %DE0 than the PSS-O since the PSS-T can provide smaller variances of the estimates by estimating uncorrelated transformed parameters23. Thus, the PSS-T is employed to estimate the uncertain parameters in the case study. As shown in Figure 5, mean squared errors in IV0 and %DE0 decrease with the increasing number of sampling points, but this is at the cost of increased experimental burden. To reduce the cost and effort for the sampling in the case study, the number of sampling points is chosen as 5, at which 90 % of the total improvement in the averaged errors of both IV0 and %DE0 is achieved. In the proposed RB2BO, the sampling times are also optimized to improve the information content in the available data, and the accuracy of the parameter estimates and extreme-case scenarios. For example, in Case 2, the sum of the squared errors between the true

ACS Paragon Plus Environment

27

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 45

and updated extreme-case scenarios of %DE0 is reduced by 30 % with the optimal sampling point design as shown in Figure 6. This also contributes to the reduction in the conservativeness of the resulting recipe: RB2BO with optimal sampling point design shows about 0.2 % of improvement in the average objective function value with the same level of constraint violations compared to RB2BO without it. The optimal sampling point design is essential to enhance the quality of the estimates and the recipe in practice when measured variables and sampling points are limited or when measurement noise levels are high.

5. CONCLUSIONS A framework of robust batch-to-batch optimization (RB2BO) has been proposed to overcome the limitations of the conventional optimization approaches under uncertainty, i.e., robust optimization and measurement-based optimization. The novelty of this framework is that scenario adaption and robust optimization are iteratively performed to alleviate the conservativeness of the recipe resulting from the conventional robust optimization while still considering the residual uncertainty after the estimation of the uncertain parameters. In the case study shown, the proposed framework provided less conservative and more robust recipes than the conventional RO and MBO, respectively. Furthermore, degree of constraint satisfaction as well as process performance was able to be improved even despite poor priors. The proposed RB2BO using the suggested termination criteria showed almost identical performance to that without the termination criteria, so the cost for measurements and optimizations could be reduced by adopting the suggested termination criteria. In this study, the proposed robust batch-to-batch optimization was applied to the case in the food industry where batch-to-batch uncertainty, e.g., feedstock variability, is considered

ACS Paragon Plus Environment

28

Page 29 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

along with off-line product quality measurements. This framework can be potentially used to handle

natural

variabilities

of

renewable

feedstocks,

e.g.,

in

biorefinery

and

cultivation/fermentation17. The suggested iterative scenario adaptation and robust optimization approach can be applied to other processes, of which performance and product quality are susceptible to parametric uncertainty, such as polymerization and pharmaceutical production processes, and can help to increase the process reproducibility after a scale-up or changes in the operating condition. Furthermore, it can be extended to more general cases such as those with on-line measurements for mid-course correction(s) of the recipe parameters. In addition, other techniques for estimating uncertain parameters or updating uncertainty distribution as well as for updating the scenario set can be employed to improve the feasibility and optimality of a solution in the proposed framework. When available data are limited and the estimation problem is illconditioned, careful scheduling of sampling is required to enhance the quality of the adaptation.

6. AUTHOR INFORMATION * Corresponding author Prof. Jay H. Lee: E-mail address: [email protected]; Tel.: 82-42-350-3926; Fax: 82-42-350-3910

7. ACKNOWLEDGEMENTS This work was supported by the Innovation Fund Denmark through the BIOPRO2 strategic research centre [grant number 4105-00020B] and by the Advanced Biomass

ACS Paragon Plus Environment

29

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 45

R&D Center (ABC) of Global Frontier Project funded by the Ministry of Science and ICT (ABC-2011-0031354).

References 1.

Bonvin, D.; Srinivasan, B.; Hunkeler, D., Control and optimization of batch

processes. Ieee Control Systems Magazine 2006, 26, (6), 34-45. 2.

Benavides, P. T.; Salazar, J.; Diwekar, U., Economic Comparison of Continuous

and Batch Production of Biodiesel Using Soybean Oil. Environmental Progress &

Sustainable Energy 2013, 32, (1), 11-24. 3.

Srinivasan, B.; Palanki, S.; Bonvin, D., Dynamic optimization of batch processes:

I. Characterization of the nominal solution. Computers & Chemical Engineering 2003, 27, (1), 1-26. 4.

Bonvin, D.; Srinivasan, B.; Ruppen, D. Dynamic optimization in the batch

chemical industry; 2001. 5.

Srinivasan, B.; Bonvin, D.; Visser, E.; Palanki, S., Dynamic optimization of batch

processes: II. Role of measurements in handling uncertainty. Computers & chemical

engineering 2003, 27, (1), 27-44. 6.

Niu, D. P.; Zhang, L.; Wang, F. L., Modeling and Parameter Updating for

Nosiheptide Fed-Batch Fermentation Process. Industrial & Engineering Chemistry

Research 2016, 55, (30), 8395-8402.

ACS Paragon Plus Environment

30

Page 31 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

7.

Kim, B.; Jang, H.; Eom, M. H.; Lee, J. H., Model-Based Optimization of Cyclic

Operation of Acetone-Butanol-Ethanol (ABE) Fermentation Process with ex Situ Butanol Recovery (ESBR) for Continuous Biobutanol Production. Industrial & Engineering

Chemistry Research 2017, 56, (8), 2071-2082. 8.

Jenzsch, M.; Gnoth, S.; Kleinschmidt, M.; Simutis, R.; Lubbert, A., Improving the

batch-to-batch reproducibility in microbial cultures during recombinant protein production by guiding the process along a predefined total biomass profile. Bioprocess

Biosyst Eng 2006, 29, (5-6), 315-21. 9.

Caroço, R. F.; Kim, B.; Santacoloma, P. A.; Abildskov, J.; Lee, J. H.; Huusom, J.

K., Analysis and model-based optimization of a pectin extraction process. Journal of

Food Engineering 2018, 192, 61-71. 10.

Caroço, R. F.; Bevilacqua, M.; Armagan, I.; Santacoloma, P. A.; Abildskov, J.;

Skov, T.; Huusom, J. K., Raw material quality assessment approaches comparison in pectin production. Biotechnol Prog 2019, 35, (2), e2762. 11.

May, C. D., Pectins. In Thickening and gelling agents for food, Imeson, A. P., Ed.

Springer Science & Business Media: 2012. 12.

Chachuat, B.; Srinivasan, B.; Bonvin, D., Adaptation strategies for real-time

optimization. Computers & Chemical Engineering 2009, 33, (10), 1557-1567. 13.

Mendoza, D. F.; Graciano, J. E. A.; dos Santos Liporace, F.; Le Roux, G. A. C.,

Assessing the reliability of different real‐time optimization methodologies. The Canadian

Journal of Chemical Engineering 2016, 94, (3), 485-497. 14.

Geletu, A.; Li, P., Recent developments in computational approaches to

optimization under uncertainty and application in process systems engineering.

ChemBioEng Reviews 2014, 1, (4), 170-190. 15.

Gorissen, B. L.; Yanikoglu, I.; den Hertog, D., A practical guide to robust

optimization. Omega-International Journal of Management Science 2015, 53, 124-137. 16.

McLean, K. A. P.; McAuley, K. B., Mathematical modelling of chemical processes

 obtaining the best model predictions and parameter estimates using identifiability and estimability procedures. Canadian Journal of Chemical Engineering 2012, 90, (2), 351366.

ACS Paragon Plus Environment

31

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

17.

Page 32 of 45

Mitsos, A.; Asprion, N.; Floudas, C. A.; Bortz, M.; Baldea, M.; Bonvin, D.;

Caspari, A.; Schafer, P., Challenges in process optimization for new feedstocks and energy sources. Computers & Chemical Engineering 2018, 113, 209-221. 18.

Luna, M. F.; Martinez, E. C., Iterative modeling and optimization of biomass

production using experimental feedback. Computers & Chemical Engineering 2017, 104, 151-163. 19.

Lopez, D. C.; Barz, T.; Korkel, S.; Wozny, G., Nonlinear ill-posed problem

analysis in model-based parameter estimation and experimental design. Computers &

Chemical Engineering 2015, 77, 24-42. 20.

Bard, Y., Nonlinear parameter estimation. New York, NY : Academic Press,

1974. 21.

Kravaris, C.; Hahn, J.; Chu, Y. F., Advances and selected recent developments

in state and parameter estimation. Computers & Chemical Engineering 2013, 51, 111123. 22.

Yao, K. Z.; Shaw, B. M.; Kou, B.; McAuley, K. B.; Bacon, D., Modeling

ethylene/butene copolymerization with multi‐site catalysts: parameter estimability and experimental design. Polymer Reaction Engineering 2003, 11, (3), 563-588. 23.

Kim, B.; Lee, J. H., Improved Parameter Estimation of Ill-Conditioned Problems.

In Proceedings of the 57th IEEE Conference on Decision and Control, Florida, USA, 2018. 24.

Kim, B.; Lee, J. H., Parameter Subset Selection in a Class of Ill-Conditioned

Estimation Problems. Journal of Process Control Submitted. 25.

Murphy, K. P. Conjugate Bayesian analysis of the Gaussian distribution; 2007.

26.

Murphy, K. P., Machine Learning: A Probabilistic Perspective (Adaptive

Computation and Machine Learning series). In The MIT Press: London, UK: 2018. 27.

Lee, K.; Marin, J.-M.; Mengersen, K.; Robert, C., Bayesian inference on finite

mixtures of distributions. In Perspectives in mathematical sciences I: Probability and

statistics, World Scientific: 2009; pp 165-202. 28.

Robertson, D. G.; Lee, J. H.; Rawlings, J. B., A moving horizon‐based approach

for least‐squares estimation. AIChE Journal 1996, 42, (8), 2209-2224.

ACS Paragon Plus Environment

32

Page 33 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

29.

Pannocchia, G.; Rawlings, J. B., Disturbance models for offset‐free

model‐predictive control. AIChE journal 2003, 49, (2), 426-437. 30.

Patwardhan, S. C.; Narasimhan, S.; Jagadeesan, P.; Gopaluni, B.; Shah, S. L.,

Nonlinear Bayesian state estimation: A review of recent developments. Control

Engineering Practice 2012, 20, (10), 933-953. 31.

Gabrel, V.; Murat, C.; Thiele, A., Recent advances in robust optimization: An

overview. European journal of operational research 2014, 235, (3), 471-483. 32.

Zhang, Y.; Feng, Y.; Rong, G., New Robust Optimization Approach Induced by

Flexible Uncertainty Set: Optimization under Continuous Uncertainty. Industrial &

Engineering Chemistry Research 2016, 56, (1), 270-287. 33.

Nagy, Z. K.; Braatz, R. D., Distributional uncertainty analysis using power series

and polynomial chaos expansions. Journal of Process Control 2007, 17, (3), 229-240. 34.

Telen, D.; Vallerio, M.; Cabianca, L.; Houska, B.; Van Impe, J.; Logist, F.,

Approximate robust optimization of nonlinear systems under parametric uncertainty and process noise. Journal of Process Control 2015, 33, 140-154. 35.

Bonvin, D.; Georgakis, C.; Pantelides, C.; Barolo, M.; Grover, M.; Rodrigues, D.;

Schneider, R.; Dochain, D., Linking models and experiments. Industrial & Engineering

Chemistry Research 2016, 55, (25), 6891-6903. 36.

Mesbah, A., Stochastic model predictive control with active uncertainty learning:

A survey on dual control. Annual Reviews in Control 2018, 45, 107-117. 37.

Franceschini, G.; Macchietto, S., Model-based design of experiments for

parameter precision: State of the art. Chemical Engineering Science 2008, 63, (19), 4846-4872. 38.

Kelco, C., GENU Pectin Book. CP Kelco U.S.

39.

Andersen, N. M.; Cognet, T.; Santacoloma, P. A.; Larsen, J.; Armagan, I.;

Larsen, F. H.; Gernaey, K. V.; Abildskov, J.; Huusom, J. K., Dynamic modelling of pectin extraction describing yield and functional characteristics. Journal of Food Engineering 2017, 192, 61-71.

ACS Paragon Plus Environment

33

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 45

List of Table Table 1. Model parameters of the dynamic model for pectin extraction via acid hydrolysis Table 2. Mean values and variances of the initial conditions obtained from the historical data of the lime 9, 10, and for Case 1 and Case 2 Table 3. Operation results of 100 batches using the recipes obtained from the MBO, RO and the proposed RB2BO for Case 1 Table 4. Operation results of 100 batches using the recipes obtained from the MBO, RO and the proposed RB2BO for Case 2

List of Figures Figure 1. Comparison between (A) conventional measurement-based optimization (two-step approach) and (B) the proposed robust batch-to-batch optimization (RB2BO) Figure 2. Overall algorithm of the proposed robust batch-to-batch optimization (RB2BO) Figure 3. Adaptation results of the mean values and extreme-case scenarios of IV0 and %DE0 in Case 1 Figure 4. Adaptation results of the mean values and extreme-case scenarios of IV0 and %DE0 in Case 2 Figure 5. The mean squared errors in IV0 and %DE0 resulting from 5,000 runs of Monte Carlo simulations using the parameter subset selection by orthogonalization (PSS-O) and by transformation (PSS-T) with nine different choices of number of sampling points during one batch; temperature=59.3 ℃ and pH=1.64. Figure 6. Adaptation results of the extreme-case scenario of %DE0 with and without optimal sampling point design in Case 2

ACS Paragon Plus Environment

34

Page 35 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Table 1. Model parameters of the dynamic model for pectin extraction via acid hydrolysis 9, 39 Parameters

Unit

Value

kmasstransfer

m/s

1.10×10-3

Da

J/mol

4.81×104

D0

m2/s

9.17×10-4

Ea,hydorlysis

J/mol

8.11×103

αhydrolysis

L/mol·s

3.25×1011

βhydrolysis

s-1

-2.6×108

Ea,degradation

J/mol

7.48×103

αdegradation

L/mol·s

2.30×10-3

βdegradation

s-1

1.80×10-4

Ea,de-esterification

J/mol

2.91×104

αde-esterification

L/mol·s

8.46

βde-esterification

s-1

4.7×10-2

Ea,IV

J/mol

4.63×104

αIV

L/mol·s

2.1×103

βIV

s-1

8.13×10

ACS Paragon Plus Environment

35

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 45

Table 2. Mean values and variances of the initial conditions obtained from the historical data of the lime 9, 10, and for Case 1 and Case 2

Mean

Historical data

peel C pectin ,0

peel 3 C protopectin ,0 (kg/m )

IV0 (g/L)

%DE0

7.88

73.93

96.04

177.23

0.1387

-0.1776

-1.440

0.7313

0.9947

-0.8976

6.713

1107

-662.0

Covariance matrix

(kg/m3)

585.5 Mean

8.17

73.04

96.89

172.09

Variance

0.0556

0.4925

591.0

325.1

Mean

8.15

75.38

94.73

192.07

Variance

0.0210

0.0937

579.0

421.3

Case 1

Case 2

ACS Paragon Plus Environment

36

Page 37 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Table 3. Operation results of 100 batches using the recipes obtained from the MBO, RO and the proposed RB2BO for Case 1 % of the constraint Expected objective function violation value excluding the (final IV/%DE) constraint violation cases

Optimization strategy

Temp. (℃)

pH

MBO

67.7 (avg.)

1.95 (avg.)

47 (18/35)

10.02

RO (using historical data)

55

1.55

1 (0/1)

8.96

RB2BO (adaptation for every batch)

61 (avg.)

1.68 (avg.)

2 (2/1)

9.61

RB2BO (adaptation until the 39th batch)

61.9

1.70

3 (3/1)

9.66

ACS Paragon Plus Environment

37

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 45

Table 4. Operation results of 100 batches using the recipes obtained from the MBO, RO and the proposed RB2BO for Case 2 % of the constraint Expected objective function violation value excluding the (final IV/%DE) constraint violation cases

Optimization

Temp. (℃)

pH

MBO

67.3 (avg.)

1.85 (avg.)

89 (33/84)

10.73

RO (using historical data)

55

1.55

86 (0/86)

9.37

RB2BO (adaptation for every batch)

61 (avg.)

1.62 (avg.)

4 (1/3)

10.15

RB2BO (adaptation until the 35th batch)

62.9

1.66

4 (1/3)

10.24

ACS Paragon Plus Environment

38

Page 39 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 1. Comparison between (A) conventional measurement-based optimization (two-step approach) and (B) the proposed robust batch-to-batch optimization (RB2BO)

ACS Paragon Plus Environment

39

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 45

Figure 2. Overall algorithm of the proposed robust batch-to-batch optimization (RB2BO)

ACS Paragon Plus Environment

40

Page 41 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 3. Adaptation results of the mean values and extreme-case scenarios of IV0 and %DE0 in Case 1

ACS Paragon Plus Environment

41

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 45

Figure 4. Adaptation results of the mean values and extreme-case scenarios of IV0 and %DE0 in Case 2

ACS Paragon Plus Environment

42

Page 43 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 5. The mean squared errors in IV0 and %DE0 resulting from 5,000 runs of Monte Carlo simulations using the parameter subset selection by orthogonalization (PSS-O) and by transformation (PSS-T) with nine different choices of number of sampling points during one batch; temperature=59.3 ℃ and pH=1.64.

ACS Paragon Plus Environment

43

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 45

Figure 6. Adaptation results of the extreme-case scenario of %DE0 with and without optimal sampling point design in Case 2

ACS Paragon Plus Environment

44

Page 45 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Table of Contents graphic

ACS Paragon Plus Environment

45