Model-based approach to the design and scale-up of wheat milling operations — Proof of concept

Accepted Manuscript Title: Model-based approach to the design and scale-up of wheat milling operations — Proof of concept Authors: Filippo Dal-Pastro, Pierantonio Facco, Eliana Zamprogna, Fabrizio Bezzo, Massimiliano Barolo PII: DOI: Reference:

S0960-3085(17)30120-7 https://doi.org/10.1016/j.fbp.2017.09.005 FBP 905

To appear in:

Food and Bioproducts Processing

Received date: Revised date: Accepted date:

16-5-2017 11-9-2017 18-9-2017

Please cite this article as: Dal-Pastro, Filippo, Facco, Pierantonio, Zamprogna, Eliana, Bezzo, Fabrizio, Barolo, Massimiliano, Model-based approach to the design and scaleup of wheat milling operations — Proof of concept.Food and Bioproducts Processing https://doi.org/10.1016/j.fbp.2017.09.005 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Model-based approach to the design and scale-up of wheat milling operations – Proof of concept Filippo Dal-Pastroa, Pierantonio Faccoa, Eliana Zamprognab, Fabrizio Bezzoa, Massimiliano Baroloa,* CAPE-Lab – Computer-Aided Process Engineering Laboratory Department of Industrial Engineering University of Padova via Marzolo, 9 - 35131 Padova PD (Italy) a

b

Bühler AG Corporate Technology Department Gupfenstrasse, 5 - 9240 Uzwil (Switerland)

_____________________ * Author to whom all correspondence should be addressed ([email protected])

Submitted to Food and Bioproducts Processing

Copyright © 2017 F. Dal-Pastro, P. Facco, F. Bezzo, E. Zamprogna, M. Barolo UNPUBLISHED No parts of this paper may be reproduced or elsewhere used without the prior written permission of the authors

Revised manuscript submitted on: September 10, 2017 Original manuscript submitted on: May 16, 2017

2

Highlights 

Mathematical model approach proposed to design/scale-up wheat milling operations



Near infrared spectra used to characterized the wheat feed material



Required experimentation moved from the industrial scale to a smaller scale



Proof of concept on industrial roller mill confirmed the potential of the approach

Abstract In this proof-of-concept study, we consider the problem of finding the conditions at which an industrial wheat milling process should be operated to provide a product with assigned quality from a given wheat variety. We propose an approach where mathematical models are used to assist the miller in determining the required operating conditions. The basic idea we explore is to move most of the experimentation from the industrial-scale equipment to a small-scale one where data can be obtained more quickly and at an inferior cost. The operation of an industrial roller mill, whose milling gap is the operating variable to be determined, is used as a test bed to assess the feasibility of the proposed methodology. An extended set of data deriving from experimentation in a small-scale mill are used jointly with a dataset from the industrial mill to design a multivariate statistical model based on latent variables, which is then inverted to find the desired operating conditions. The results show that the proposed approach is viable and is particularly effective when the industrial-scale dataset is limited. It is also shown that near infrared spectra can be effectively used to characterize the wheat feed.

_________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

3

Keywords: wheat milling; roller mill; scale-up; PLS; joint-Y PLS.

1. Introduction A wheat kernel is composed of three parts: endosperm (the starchy part from which flour is produced), germ (the embryonic plant) and bran (protective layers that cover the endosperm and the germ). The main objective of a wheat milling process is to recover, at minimum cost, as much endosperm as possible with an appropriate particle size distribution and with minimum amount of bran and germ. To achieve this goal, industrial wheat milling processes employ a gradual size reduction approach, through repeated milling operations (by roller mills) followed by sifting operations (by plansifters). As an example, a milling process may exploit up to seventeen different steps, arranged in a complex configuration (Campbell, 2007). In the operation of industrial wheat milling plants, two problems often need to be faced. When a plant is to be operated for the first time, very limited knowledge about the process is available, often deriving from preliminary experimental campaigns carried out in smaller-scale (e.g., laboratory or pilot) equipment. The issue therefore arises about whether the small-scale plant data can be used to find the conditions (“new operating conditions”) at which the industrialscale plant should be operated to obtain a product of assigned quality from a given wheat variety. In this study, we will refer to this problem as to a process scale-up problem. On the other side, if the industrial-scale plant has long been in operation, it may be required to determine the (new) operating conditions that are needed to manufacture a product with a quality profile that is different from that of products already obtained, or to manufacture a given product from a wheat variety that is different from those used in previous production campaigns. This second problem, which is in essence similar to the scale-up one, may be denoted as a process design problem. _________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

4

In this study we assess the potential of alternative methodologies to address the process scaleup or process design problems in industrial wheat milling processes. Given the proof-of-concept of the study, the analysis is not extended to an entire milling plant, but is limited to a roller mill. Namely, the objective is estimating the milling gap to be used in an industrial-scale roller mill to obtain a flour product with a desired mass distribution from an assigned wheat variety. Additionally, we assess the possibility of using near infra-red spectra to characterize a given wheat feed. From a practical standpoint, the process scale-up problem, as well as the process design one, can be tackled by carrying out an extensive experimental campaign on the industrial-scale plant, in such a way that the required new operating conditions can be estimated based on the experiments and on the miller’s experience. However, besides being time consuming, this approach may be very expensive due to the potential off-specification of the resulting flour product if the new operating conditions are specified incorrectly, as well as to the significant cost of the industrial experiments. As a matter of fact, the cost of conducting a single experimental run in an industrial roller mill is about four times greater than conducting an experimental run in a small scale mill (Bühler AG, 2016; Molino Quaglia, 2016). Therefore, developing a process scale-up/design methodology that allows one to reduce the number of experiments that need to be carried out at the industrial level would be greatly beneficial. The basic idea we explore in this study is moving most of the experimentation from the industrial-scale equipment to the small-scale one. Particularly, we assess the possibility to reduce the industrial-scale experimental effort by creating an extended dataset at the small-scale level that can be used to complement the industrial-scale one. Under this perspective, the smallscale dataset would act as a sort of “universal” historical dataset, which is collected once and for all by the plant manufacturer, and that can then be used to support the scale-up/design exercise on any industrial-scale plant produced by that manufacturer. Notice that the advantage _________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

5

of this approach would be twofold: on the one side, the cost of experimental trials at the industrial site (the customer) would be significantly reduced; on the other side, the commissioning of the industrial plant would be accelerated, even in the case of limited experience of the customer on a particular plant or product. To the best of the authors’ knowledge, so far no attempt has been done to investigate the scale-up/design problem for a wheat roller mill under this perspective. The use of historical data to assist product and process design is not new in the field of industrial chemical processes. Moteki and Arai (1986) used a latent-variable modeling methodology (namely, principal component analysis; PCA) to analyze historical data from a low-density polyethylene process and to infer process conditions for new grades of products. Jaeckle and MacGregor (1998 and 2000) developed a methodology based on the inversion of latent-variable models to estimate the required process conditions for a plant to yield a desired set of properties in the final product. They related the process variable datasets available at different process scales through the common products manufactured at the different scales, and they used latent-variable regression model inversion to estimate the process conditions for a target plant to manufacture a new product of assigned properties. This methodology was then used in several studies. Zhou and Titchener-Hooker (1999) determined operational windows in a bioprocess by plotting pairs of process variables. Chen and Wang (2000) used the same approach with historical data from a fluid catalytic cracking unit, but with multivariate operational regions defined by the scores of a PCA model. The methodology was further refined by Garcìa-Muñoz et al. (2005), which developed the joint-Y projection to latent structures (JYPLS) modeling approach. JY-PLS is a latent-variable regression modeling technique that allows relating datasets of regressors (e.g., raw materials properties and processing conditions) from different sources through the latent space generated by the response variables, which account for the product quality. Liu et al. (2011) used JY-PLS to support the scale-up of a roller _________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

6

compaction process in a pharmaceutical product. Facco et al. (2012) used different latent variables techniques (including JY-PLS) to transfer a process monitoring model between different plants. Tomba et al. (2014) used the JY-PLS model approach to transfer a nanoparticle product between mixers of different geometry. An application of latent-variable modeling to support process scale-up has been reported by Muteki et al. (2011), who combined PLS model design and inversion to estimate the end point for a high-shear wet granulation process in a target plant using data archived from past batches and drug lot properties. This study investigates three alternative process scale-up approaches, all based on latentvariable modeling (Burnham et al., 1999). The approaches differ for the extension of the experimental domain covered by the industrial-scale dataset, which is used to determine the operating conditions required to obtain a product with assigned quality. The relative merits of the three approaches are discussed using data deriving from experimentation on an industrialscale roller mill and a small-scale one.

2. Materials and methods 2.1. Dataset generation The data-driven models discussed later in this section are meant to relate a set of inputs (generally represented by matrix X) to a set of outputs (indicated by matrix Y). The inputs may include process operating parameters (e.g., milling gap, mass flow) and raw materials properties (e.g. wheat variety, moisture content, composition fingerprint). The outputs typically refer to a given quality target for the product under investigation. A key issue in data-driven modeling is appropriate organization of the available data. Details on data organization for the small-scale and industrial-scale mills considered in this study (whose main characteristics are summarized in Table 1) are reported below. _________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

7

Table 1. Main characteristics of the small-scale and industrial-scale roller mills.

Roll length [m] Roll diameter [mm] Feed mass flow (up to) [kg/h] Roll speed [rpm] Roll differential Roll configuration

Small-scale roller mill 0.10 250 300 500 1:2.6 Dull to dull

Industrial-scale roller mill 1.0 250 8350 550 1:2.5 Dull to dull

2.1.1 Small-scale dataset Experiments were designed following a full factorial 23 design of experiments (DoE; Montgomery, 2013) applied to three factors: feed mass flow (MF, kg/h), milling gap between the rolls (MG, µm) and wheat kernel moisture content (MC, %). The ranges for MF, MG and MC were chosen to match typical industrial ranges. Their actual values are not reported here for confidentiality reasons. Each experiment was repeated twice, and average response values were used. In order to take into account the variability of the raw material, the design was repeated for three different wheat varieties (a pure hard, semi-hard and soft variety respectively, all harvested in Switzerland in 2014). The composition fingerprint of each wheat sample feeding the roller mill was characterized using a near infra-red (NIR) spectrometer (TANGO, Bruker Optics, Massachusetts, U.S.A.). Each NIR measurement was repeated three times, and average spectra were used. The product quality was characterized with the mass distribution over a 1120 µm sieve downstream the roller mill. process The input matrix X SS [24×953] for the small-scale (SS) machine was split into an XSS





process NIR NIR process matrix and an XSS matrix; hence XSS  XSS . XSS [24×3] collects the values of 2 , XSS

process parameters (MF and MG) and one wheat property (MC) for 24 different observations. NIR [24×950] collects the absorbance corresponding to 950 wavenumbers that describe the XSS

_________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

8

24 NIR spectra of the feed materials. The wavenumbers range from 11536 cm-1 to 3944 cm-1, evenly spaced every 8 cm-1. The mass fractions wi at the roller mill outlet were collected in YSS [24×2]. Table 2 summarizes the variables collected in the small-scale roller mill. Table 2. Variables collected in the small-scale roller mill.

YSS

XSS no. 1

MF [kg/h]

name

no. 1

2

MC [%]

2

3 from 4 to 953

name

w1 w2

MG [µm] wavenumbers of the NIR spectrum [cm-1]

2.1.2 Industrial-scale dataset Typically, in industrial-scale roller mills the MG cannot be assigned directly as a distance (in [µm]), but is specified by setting a value that can be adjusted directly on the roller mill by means of two rotating knobs (sometimes called “clocks”). In order to counteract the imperfect alignment of the roll pair, slightly different MGs may be set at each end of the roll lengths. These two MGs will be referred to as MGl (left end) and MGr (right end). Experiments were designed in an industrial mill with a 21 DoE on the MG only. The DoE was repeated for three wheat varieties (denoted by Variety 1, Variety 2 and Variety 3). Variety 1 was a mixture of hard wheat coming from North America, Austria and Italy. Variety 2 was a semi-hard mixture, with wheat coming from Slovakia and Italy. On the other hand, Variety 3 is composed of soft wheat cultivated in Hungary and Italy. All wheat was harvested in 2015. Each variety was characterized with a NIR spectrum, collected with a FOSS spectrometer (FOSS, Hillerød, Denmark).

_________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

9

Following the same rationale as in the small-scale plant, the input X IS [12×1052] for the



industrial-scale (IS) machine was split into two matrices: XIS  XIS

process



, XISNIR . X process [12×2] IS

collects the values of MGl and MGr. X ISNIR [12×1050] collects the 1050 wavelengths in nm (from 400 nm to 2498 nm, evenly spaced every 2 nm) describing the wheat type for 12 different experiments. The response matrix YIS [12×2] includes the mass fractions wi of the milled material. The variables collected in the industrial-scale roller mill are listed in Table 3. Table 3. Variables collected in the industrial-scale roller mill.

X IS no.

YIS

1

MGl [-]

name

no. 1

2

MGr [-]

2

from 3 to 1052

name

w1 w2

wavelengths of NIR spectrum [nm]

The twelve available experiments were partitioned into a calibration set and a validation set. Two different situations are analyzed in this study, which differ for the calibration dataset actually used to build the regression models. In Case study 1 (Figure 1), the validation samples (i.e., the new products to be obtained) lie within the experimental domain explored by the calibration samples. Therefore, in this case no model extrapolation is required.

calibration validation

MG

Variety 1

Variety 2

Variety 3

Figure 1. Case study 1: structure of the industrial-scale dataset. The calibration samples (circles) and validation samples (squares) are highlighted. _________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

10

In Case study 2 (Figure 2), a limited industrial-scale dataset was assumed to be available. Namely, only three industrial-scale observations were used (circles in Figure 2), two of which belonging to wheat Variety 1 and one to Variety 2. The validation observations pointed to two Variety 3 observations (squares in Figure 2), meaning that the products to be obtained lie in a subspace of inputs not covered by the available industrial experimentation. On the other hand, it is assumed that a small-scale experimental dataset (shaded dashed line in Figure 2) is available, which is much more extended than the industrial one and covers also the input domain expected for the new products.

calibration validation small-scale experimental domain

MG

Variety 1

Variety 2

Variety 3

Figure 2. Case study 2: sketch of the situation where only a limited dataset is available from the industrial-scale plant. The available industrial-scale data (circles) cover only two wheat varieties and two “large” milling gaps. The available small-scale data (shaded dashed square) cover a wide domain of milling gaps and wheat varieties. The products to be obtained (small squares) lie within a subspace not covered by the industrial-scale dataset.

2.2 Mathematical methods This section describes the latent-variable modeling techniques that have been used in this study. They include principal component analysis (PCA), projection to latent structures (PLS) and joint-Y PLS (JY-PLS). Additionally, the optimization framework utilized to invert the latentvariable models is presented.

_________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

11

2.2.1 Principal component analysis (PCA) Principal component analysis (PCA; Jackson, 1991) is a multivariate statistical technique used to summarize the information embedded in a dataset X of (correlated) variables, by projecting the data onto a new space of latent variables that optimally describe the variability of the data as well as the correlation among them. The latent variables identify the directions of maximum variability in the data and are called principal components (PCs). The dataset X [I×N], where I is the number of observations and N the number of variables, is represented as the product of a score matrix T [I×A] and the transpose (superscript T) of a loading matrix P [N×A]: X  TP T  E .

(1)

Matrix E contains the residuals obtained when the X matrix is built using the A PCs and accounts for the model mismatch. The scores T are the coordinates of the data samples in the new reference system represented by the PCs, and represent how the observations relate to each other. The loadings P of the model are the director cosines between the latent directions and the axes of the original variables space, and can be used to analyze the correlation between the original variables.

The data must be pretreated prior to the analysis. Data pretreatment is essential both to avoid the identification of variability directions due to differences between the variable mean values (Wise et al., 2006) and to make the analysis independent of the measurement units. In this study, auto-

scaling was used, namely the variables were mean-centered and scaled to unit variance. The NIR spectra were pretreated by standard normal variate, although other possibilities exist (van den Berg et al., 2006; Rinnan et al., 2009). The most appropriate number A of PCs to be selected was determined by the SCREE test (Valle et al., 1999).

_________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

12

2.2.2 Projection to latent structures (PLS) Projection to latent structures (PLS; Höskuldsson, 1988) is a multivariate correlative regression method used to relate an X [I×N] matrix of input variables to a Y [I×M] matrix of output responses. The relationship is based on variable projection onto a common space of uncorrelated variables called latent variables (LVs). The LVs explain the major sources of systematic variability of the inputs that are mostly correlated to the variability of the responses. PLS models X and Y as follows: X  TPT  EX ,

(2)

Y  TQT  EY ,

(3)

T  XW ,

(4)

where T [I×A] is the score matrix, P [N×A] and Q [M×A] are the loading matrices and W [N×A] is the weight matrix. EX [I×N] and EY [I×M] are the residual matrices accounting for the model mismatch in the X and Y spaces, respectively. A represents the number of significant LVs chosen to build the model. The same considerations made about PCA scores and loadings are valid also for PLS. In this study, the data were pretreated by auto-scaling, and the SCREE test was used to select the number of LVs.

2.2.3 Joint-Y PLS (JY-PLS) Joint-Y PLS (JY-PLS; Garcìa-Muñoz et al., 2005) is a latent-variable modeling technique that allows relating two or more regressor datasets (e.g., XA [I×NA] and XB [J×NB]) through the joint space formed by their corresponding M response variables datasets (e.g., YA [I×M] and YB [J×M]).

The idea under JY-PLS is that, if similar products are produced in different plants (plant A and plant B), the correlation structure of the response datasets is expected to be the same, and it is identified by a common latent space between YA and YB. This common space of the product _________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

13

quality will be spanned by (part of) the latent space of the regressors, assuming that the regressors are correlated to the response variables (within-plant correlation). In addition, this product quality common space can be used to relate the different datasets and to transfer information between them, as occurs for example when the data come from different plants or different plant scales (between-plant correlation). To find this common latent space, the available datasets are decomposed onto their latent structures in order to maximize at the same time the squared covariance between XA and YA, and between XB and YB, together with the squared joint covariance between them, as identified by the joint response variables datasets YJ  [YAT YBT ]T (Garcìa-Muñoz et al., 2005):

Y  T  YJ   A    A Q TJ  E YJ ,  YB  TB 

(5)

XA  TA PAT  EXA ,

(6)

XB  TBPBT  EX B ,

(7)

TA  XA WA* ,

(8)

TB  XBWB* ,

(9)

where QJ [M×A] represents the joint loading matrix defining the common latent space of YJ (A is the number of LVs used to build the model). The meaning of the other symbols is similar to the one of the PLS models. Note that the JY-PLS modeling framework does not impose any restriction on either the number of columns in XA and XB (i.e., plants A and B may have a different number of measured variables), or on the rows of XA and XB (i.e., the number of observations available in the two plants may be different). Besides sharing the same correlation structure for YA and YB, the only additional restriction is that the number of columns in the response variable matrices be the same. The NIPALS algorithm can be used to compute the JY-

_________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

14

PLS model parameters (Garcìa-Muñoz et al., 2005); in this case, a particular scaling must be applied in addition to auto-scaling:

XA 

XA , IN A

(10)

XB 

XB , JN B

(11)

YA (12) , I Y (13) YB  B , J where I and J are the number of observations in plant A and plant B, respectively, NA is the YA 

number of input variables measured in plant A, and NB the number of input variables of plant B.

2.2.4 Latent-variable model inversion A latent-variable model (LVM) is usually used to predict a set of response variables yˆ NEW [1×M] given a set of regressors x NEW [1×N]. This can be defined as the direct use of the model (Figure 3a). However, if a LVM based on historical data is available and a set of desired response variables yDES [1×M] is assigned, the LVM can be mathematically inverted to estimate the set of regressors xNEW [1×N] that, according to the model, provide the desired response variables (Figure 3b). yˆ 1NEW

x1NEW

... x

NEW N

x

NEW 1

LATENT VARIABLE MODEL

... yˆ MNEW

(a) ... x

NEW N

y 1DES

LATENT VARIABLE MODEL-1

... y DES M

(b) Figure 3. (a) Direct and (b) inverse use of latent-variable models. _________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

15

In this study, the problem of LVM inversion is solved by means of an optimization framework. After data pretreatment, the optimization problem is formulated as follows (Tomba et al., 2012):

 DES DES min  y  yˆ x NEW 



 y T

DES

 yˆ

DES



 A  x NEW w a  G1     a 1  sa 

  

2

    G  SPE  2 x   

subject to : (14)

yˆ DES  x NEW WQ T



SPE x  x NEW  x NEW WP T x

NEW n

c

eq n

 x T

NEW

 x NEW WP T



,

where xNEW is the estimated set of regressors required to obtain yDES, W are the weights of the LVM, Q are the Y loadings of the LVM, G1 and G2 are weighting parameters, wa is the ath column of W, sa is the variance of the ath column of the score matrix T of the LVM, P are the X loadings in the LVM model, and cneq is an equality constraint on element xnNEW of the estimated set of N regressors (a total number   N of x NEW elements may be constrained). Three addenda can be identified in the objective function: 

the first addendum is included to obtain a set of regressors xNEW that provides, according to the model, the desired response variable yDES;



the second addendum represents a soft constraint (i.e., a penalty function) on the Hotelling’s T2 of the solution. The Hotelling’s T2 (Hotelling, 1933) measures the distance of the scores of the solution from the origin of the LV space. This soft constraint keeps the solution close to the average values of the historical dataset X used to build the model. G1 is a weighting parameter, which can be set equal to the inverse of the 95% Hotelling’s T2 limit of the LVM model;



the third addendum identifies a soft constraint on the SPE of the estimated set of regressors. The SPE statistics describes how well an observation is represented by the model. This soft _________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

16

constraint forces the solution to lie close to the model space. G2 is a parameter used to weight the constraint, and can be set equal to the inverse of the 95% SPE limit of the LVM model.

3. Model-based approaches to process design/scale-up This section describes the methodology that is proposed to support process design/scale-up exercises in industrial-scale wheat roller mills. The main objective is estimating the milling gap to be used in an industrial-scale roller mill to manufacture a product with a desired mass distribution over a sieve system from a specific wheat variety, while at the same time minimizing the experiments that need to be carried out in the industrial plant in order to determine the milling gap. To this purpose, three different modeling approaches are investigated: direct PLS modeling, PLS model inversion, and JY-PLS model inversion. The first two approaches use industrial-scale data only, whereas in the latter one the industrial-scale dataset is complemented with a dataset coming from a small-scale roller mill. This latter approach allows assessing whether data integration can accelerate the scale-up exercise.

3.1 Using PCA scores as a proxy of NIR spectra Previous studies (Dal-Pastro et al., 2016) showed that the NIR spectra collected in the smallscale plant are able to discriminate different wheat varieties as well as different MCs in the feed material. The same feature was tested in the industrial-scale plant. Figure 4 shows the score plot of the PCA model built on the NIR spectra of the wheat used in the industrial-scale experimental campaign (the wheat material was completely unrelated to the one used in the small-scale plant). The pattern is similar to that observed in the small-scale plant (Dal-Pastro et al., 2016), but wheat Varieties 1 and 2 are less separated. This is due to the fact that in the industrial plant either variety was obtained as a blend of several constituents, and most of the constituents of _________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

17

Variety 1 and Variety 2 were common. Additionally, the variability observed in the second component is not due to the moisture content as in the small-scale plant (the moisture content was kept constant), but to a variability of the wheat blend. 15 Variety 1 Variety 2 Variety 3

10

PC2 [9.8%]

5 0 -5 -10 -15 -20 -30

-20

-10

0

10

20

30

40

50

PC1 [85.5%]

Figure 4. Score plot of the PCA model built on NIR spectra for the industrial-scale plant.

Two PCs provide a very good representation of the NIR spectra variability, since they overall describe over 95 % of the NIR spectra variance for both the small-scale plant and the industrialscale one. Therefore, to reduce the number of inputs that need to be processed, XNIR (whose column dimension is equal to the number of wavelengths, hence very large) is compressed by PCA modeling, using two PCs at each plant scale. The corresponding score matrices will be wheat denoted as TSS and TISwheat for the small-scale and industrial scale plants, respectively. The

PCA scores are therefore used as a proxy of the whole spectra to describe the wheat variety. process wheat For an assigned new product, the input matrices will be indicated by X'SS  [XSS , TSS ] for

the small-scale plant, and X'IS  [Xprocess , TISwheat ] for the industrial-scale plant. Note that the IS scores are also used to assign the wheat variety in the constrained optimization problem of Eq. (14) (to be discussed in detail later).

_________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

18

3.2 Approach no. 1: direct PLS modeling As sketched in Figure 5, the input variables of the PLS model are the mass distribution over a sieve (collected in YIS ) and the wheat variety (described by TISwheat ). The output variables are the milling gaps, collected in X process . Therefore, according to Approach no. 1 the milling gap IS is estimated directly for the given wheat variety and desired mass distribution of the product. Clearly, only industrial-scale data can be exploited to build the PLS model, which is formulated as:

Y

IS



, TISwheat  TP T  EX

(15)

X process  TQT  E Xprocess . IS

(16)

IS

TISwheat wheat variety

Y IS mass

PLS MODEL

X process IS milling gap

(industrial-scale)

Figure 5. Approach no. 1: model inputs and outputs.

3.3 Approach no. 2: PLS model inversion Also Approach no. 2 relies on PLS modeling using industrial-scale data only. First, a PLS model is built for the industrial-scale plant by relating X'IS (i.e., milling gap and wheat variety) to YIS (mass distribution over a 1120 µm sieve), as shown in Figure 6a. It may be noticed that the way the model is built reflects how the roller mill works (namely, milling gaps and wheat variety are the input variables, and the mass distribution is the response). Once the model is calibrated, it is then inverted to estimate the milling gap to be used to obtain a desired product, for a given wheat variety (Figure 6b). To this purpose, a constrained optimization problem is formulated and solved as in Eqs. (19)-(22).

_________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

19

TISwheat

t wheat

PLS MODEL

wheat variety

X process IS milling gap

wheat variety

Y IS mass

x NEW, process milling gap

(industrial-scale)

(a)

PLS MODEL-1

y DES mass

(industrial-scale)

(b)

Figure 6. Approach no. 2: model inputs and outputs for (a) the direct model and (b) model inversion.

Namely, the PLS model is formulated as follows: X'IS  TP T  EX'

(17)

YIS  TQT  EYIS .

(18)

IS

To accomplish the model inversion, the optimization problem (14) is reformulated as:   DES DES min  y  yˆ x NEW  



 y T

DES

 yˆ

DES



  A  x NEW w 2   a    G2SPE X   G1     a 1  sa      

(19)

subject to :

yˆ DES  x NEW WQ T



SPE X  x NEW  x NEW WP T

(20)

 x T

NEW

 x NEW WP T



t NEW,wheat  t wheat .

(21) (22)

The different addenda in the objective function were described in detail Section 2.2.4. The set of input variables is x NEW  [x NEW,process, t NEW,wheat ] , where x NEW,process is the set of estimated process parameters that should be used to obtain the desired product, and t NEW,wheat are the scores describing the wheat variety for the new product to be obtained. The latter are assigned in (22) through the equality constraint t NEW,wheat  t wheat , where t wheat are the scores obtained by projecting the NIR spectrum onto the PCA model built on the NIR spectra of the calibration samples. Note that, by using constraint (22), only A scores are required to assign the wheat feed _________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

20

variety. Using a whole NIR spectrum instead would have required the assignment of ~1000 wavenumbers, making the optimization problem much more challenging.

3.3 Approach no.3: joint-Y PLS modeling This approach is based on a JY-PLS model, which exploits data from different plant scales: a small-scale roller mill and an industrial-scale roller mill. The inputs to the small-scale equipment are described by X'SS , whereas the outputs are described by YSS . A similar notation is used for the industrial-scale (subscript IS) equipment. Note that the industrial-scale and smallscale input matrices have a different number of columns, due to the fact that different variables are measured in the two plants. Therefore, the two input matrices cannot be simply merged into an augmented matrix to be used as the input to a standard PLS model. A JY-PLS modeling framework is used instead, since it can naturally cope with this data structure. The model is calibrated as shown in Figure 7a, and then inverted as in Figure 7b to estimate the milling gap. Model inversion is carried out according to (14), in which constraints on the PCA model scores are used to assign the input wheat variety.

wheat SS

T

wheat variety process X SS

milling gap

JY-PLS MODEL-1

JY-PLS MODEL

Smallscale

Y

TISwheat wheat variety

X process IS

Smallscale

SS mass

t wheat

Industrialscale (a)

Y IS mass

wheat variety

x NEW, process milling gap

Industrialscale

y DES

(b)

Figure 7. Approach no. 3: model inputs and outputs for (a) the direct model and (b) model inversion.

_________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

mass

21

The model is formulated as†:

Y  T  YJ   SS    A QTJ  EYJ  YIS  TB 

(23)

' XSS  TA PAT  EX'

(24)

X'IS  TB PBT  EX' .

(25)

SS

IS

The optimization problem (14) is formulated as:   DES min y  yˆ DES NEW  x  



 y T

DES

 yˆ

DES



 LV  x NEW w B,lv  G1     lv1  s lv  

   

2

    G SPE  2 X    

(26)

subject to :

yˆ DES  x NEW WBQTJ



SPE X  x NEW  x NEW WB PBT t NEW, wheat  t wheat .

(27)

 x T

NEW

 x NEW WBPBT



(28) (29)

4. Results and discussion In this section, the performances of the three model-based approaches for process scale-up are discussed. It should be recalled that two milling gaps need to be estimated, because in the industrial roller mill two independent knobs are used to set the milling gaps on either side of the roller pair. As discussed in Section 2.2, two alternative occurrences are analyzed, depending on whether an extended dataset (Case study 1) or a limited dataset (Case study 2) are available from the industrial-scale plant.

†

The same symbols for scores (T), loadings (Q) and weights (W) are used for all three approaches, even though each matrix is made of different elements in each approach. _________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

22

4.1 Model design and diagnostics

4.1.1 Case study 1: extended industrial dataset Table 4 summarizes the explained variance for the first two modeling approaches in Case study 1. Two LVs describe almost the entire variability of the responses.

Table 4. Case study 1 (extended industrial dataset): explained variance for the first two modeling approaches (both models are built using 2 LVs).

RX2 (input) [%]

RY2 (output) [%]

75 74

98 98

Approach no. 1 Approach no. 2

The joint-Y PLS model of Approach no. 3 is built using 2 LVs and the variability described by each LV is reported in Table 5. The first LV describes the entire variability of the response variables, whereas the second LV is used to describe a significant amount of the variability of X'IS and of X'SS .

Table 5. Case study 1 (extended industrial dataset): explained variance for each LV of the JY-PLS model in modeling Approach no.3.

RX2 ' [%]

RX2 ' [%]

RY2IS [%]

RY2SS [%]

59 15

20 36

98 0

97 0

IS

LV1 LV2

SS

Figure 8 shows the PA and PB loadings for the first latent variable of the JY-PLS model. For both plants, the main source of variability is the milling gap. Since the first component (that is described by the milling gap) explains almost the entire variability of the response variables, the milling gap is the most influencing process parameter on the responses.

_________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

23

0.6

0.6

0.4

0.4

0.2

PB

PA

0.0

0.2

-0.2 -0.4

0.0

-0.6

-0.2

MF

MC

MG

Score 1 Score 2

-0.8

MGl

MGr

(a)

Score 1

Score 2

(b)

Figure 8. Case study 1: loadings of the JY-PLS model. (a) PA and (b) PB.

4.1.2 Case study 2: limited industrial dataset Table 6 lists the explained data variability for the first two modeling approaches in Case study 2. Two LVs are used to calibrate the models, and they describe almost the entire variability of both regressors and responses.

Table 6. Case study 2: explained variance for the first two modeling approaches (both models are built using 2 LVs).

Approach no. 1 Approach no. 2

RX2 (input) [%]

RY2 (output) [%]

99 99

99 99

The joint-Y PLS model of Approach no. 3 is built using 2 LVs and the variability described by each LV is reported in Table 7. Two LVs are able to explain almost the entire variability of

X'IS , as well as 56 % of the X'SS variability. Interestingly, the first LV explains a very large fraction (90%) of the variability of the industrial-scale input matrix X'IS , but only a limited amount (20%) of the small-scale input matrix X'SS variability; a significant amount of the X'SS variability (36%) is explained by the second LV. On the other side, a single LV describes almost the entire variability of the response matrices in both plants.

_________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

24

Table 7. Case study 2 (limited industrial dataset): explained variance for each LV of the JY-PLS in modeling Approach no.3.

RX2 ' [%]

RX2 ' [%]

RY2IS [%]

RY2SS [%]

90 9

20 36

99 0

98 0

IS

LV1 LV2

SS

The loading analysis is similar to Case study 1 and therefore is not reported.

4.2. Results for Case study 1 The milling gap prediction results for the three approaches are reported in Table 8 and Table 9. The objective was obtaining a product with mass distribution yDES with the wheat variety specified in the column named as “Wheat variety”. Table 8 lists the real milling gap used in the industrial roller mill to obtain that product (column 3) and the milling gap estimations from the three modeling approaches (columns 4, 5 and 6; note that only one milling gap is reported for the sake of conciseness, even though two are estimated). Table 9 reports the relative errors for the estimations shown in Table 8. Several different new products were considered, each one reported in the first row of Table 8 and Table 9. When a given new product was being considered, information about it was obviously not used to build the PLS or JY-PLS models. Therefore, the results in Tables 8 and 9 are validatory for each process scale-up exercise.

Table 8. Case study 1 (extended industrial dataset): validatory milling gap estimations for the three proposed approaches to obtain a desired yDES mass distribution from an assigned wheat variety. For each desired product, the approach providing the best performance is indicated in boldface. yDES [79.4,20.6] [74,26] [80.5,19.5] [74.5,25.5] [80,20] [72.4,27.6]

Wheat variety Variety 1 Variety 1 Variety 2 Variety 2 Variety 3 Variety 3

MGreal 6.25 7.17 6.25 7.17 6.25 7.17

MGApproach 1 6.16 6.89 5.95 6.82 6.23 6.84

MG Approach 2 6.16 6.91 5.95 6.83 6.23 6.83

MG Approach 3 6.08 6.85 5.88 6.77 6.09 6.85

_________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

25

Table 9. Case study 1 (extended industrial dataset): relative milling gap relative estimation errors (RE) for the examples of Table 8. yDES [79.4,20.6] [74,26] [80.5,19.5] [74.5,25.5] [80,20] [72.4,27.6]

Wheat variety Variety 1 Variety 1 Variety 2 Variety 2 Variety 3 Variety 3

REApproach 1 [%] 1.4 3.8 6.7 4.8 0.2 4.5

REApproach 2 [%] 1.3 3.6 4.7 4.6 0.2 4.7

REApproach 3 [%] 2.6 4.4 5.8 5.5 2.5 4.5

All approaches perform very good milling gap estimations, with relative estimation errors well below 7 %. Overall, slightly better estimations are provided by Approach no. 2 (PLS model inversion); in fact, the average relative error is 3.2 % for Approach no.2, whereas it is 3.6 % for Approach no.1 and 4.2 % for Approach no.3. We conclude that, in this case, including the small-scale data in the modeling framework (JY-PLS model) does not provide a major advantage. Nevertheless, using JY-PLS modeling the industrial-scale experimental effort can be somewhat reduced, and the milling gap estimations are still acceptable.

4.3 Results for Case study 2 Results for the three modeling approaches are shown in Table 10 and Table 11 for two different new products. It can be observed that the JY-PLS model inversion approach (Approach no. 3) largely outperforms the other two approaches in estimating the milling gap that is required to obtain the desired product. Notably, Approach no. 2 (PLS model inversion) did not even provide a result: due to the limited domain explored by the industrial-scale experimentation, the PLS model inversion turned out to be unfeasible. However, despite the fact that the industrial-scale data covered a very narrow experimental domain only, complementing them with data from a small-scale plant (Approach no. 3) enabled accurate estimation of the milling gap required to obtain a product lying in a region not explored by the industrial experimentation.

_________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

26

Table 10. Case study 2 (limited industrial dataset): validatory milling gap estimations for the three proposed approaches to obtain a desired yDES mass distribution from an assigned wheat variety. For each desired product, the approach providing the best performance is indicated in boldface. yDES [62.7,37.3] [72.4,27.6]

Wheat variety Variety 3 Variety 3

MGreal 8.00 7.17

MGApproach 1 6.71 6.01

MGApproach 2 constraints not satisfied constraints not satisfied

MGApproach 3 7.83 6.65

Table 11. Case study 2 (limited industrial dataset): relative milling gap relative estimation errors (RE) for the examples of Table 10. yDES [62.7,37.3] [72.4,27.6]

Wheat variety Variety 3 Variety 3

REApproach 1 [%] 16.1 15.8

REApproach 2 [%] unfeasible unfeasible

REApproach 3 [%] 2.1 7.2

Therefore, Approach no. 3 is recommended when only limited experimentation is available from the industrial-scale plant, but a dataset from a similar small-scale plant (or possibly a similar industrial-scale plant located elsewhere) is jointly available that fully covers the input domain wherein the operating conditions required to obtain the desired product are expected to lie.

5. Conclusions In this proof-of-concept study, a mathematical model-based approach to the problem of process design/scale-up in wheat milling has been proposed and validated using an industrial roller mill as a test bed. The objective was to calculate the milling gap to be used in an industrial mill to obtain a product with a desired mass distribution over a sieve from a given wheat variety. The underlying ideas explored in the study were: i) providing a model-driven alternative to standard experience-driven approaches to the design/scale-up of wheat milling operations; ii) assessing the possibility to use near infrared spectra to characterize a given wheat feed; iii) assessing the possibility to move part of the required experimentation from the industrial plant to a smaller-

_________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

27

scale one, where it is quicker to carry out, less expensive and amenable to being shared across different plants. It was shown that the milling gap can effectively be estimated using data-driven models, thus opening the route towards the development of intelligent grain milling systems where the operating conditions required to obtain a given new product can be determined automatically (rather than empirically) on the basis of historical information on products already manufactured. Additionally, NIR spectra were shown to be able to characterize the feed variety and moisture content. Finally, it was shown that data coming from different sources can be merged together to build a single model from which the required operating conditions can be calculated. This means that there indeed exists the potential to move (at least) part of the experimentation from the industrial-scale level (the mill user’s side) to a smaller-scale level (the mill manufacturer’s side), with related savings in time and money. Although more experimentation is required before being able to fully transfer this novel technology to an entire industrial system, we believe that the results presented in this study can boost further research in this direction and stimulate practitioners to consider the use mathematical models to optimize the design and operation of wheat milling processes.

Acknowledgement The authors would like to express their gratitude to Mr. Andrea Quaglia and Mr. Lucio Quaglia for allowing them to carry out the experimentation in the Molino Quaglia industrial plant at Vighizzolo d’Este (Padova; Italy).

References Bühler AG, 2016. Personal communication. http://www.buhlergroup.com. _________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

28

Burnham, A.J., MacGregor, J.F., Viveros, R., 1999. Latent variable multivariate regression modeling. Chemom. Intell. Lab. Syst. 48, 167–180. Campbell, G.M., 2007. Roller milling of wheat, in: Handbook of Powder Technology. Elsevier Ltd, Amsterdam (Netherlands). Chen, F.Z., Wang, X.Z., 2000. Discovery of operational spaces from process data for production of multiple grades of products. Ind. Eng. Chem. Res. 39, 2378–2383. Dal-Pastro, F., P. Facco, F. Bezzo, E. Zamprogna and M. Barolo (2016). Using PLS and NIR spectra to model the firstbreakage step of a grain milling process. In: Computer-Aided Chemical Engineering 38, 26th European Symposium on Computer Aided Process Engineering (Z. Kravanja and M. Bogataj, Eds.), Elsevier, Amsterdam (The Netherlands), p.1171-1176. Facco, P., Tomba, E., Bezzo, F., García-Muñoz, S., Barolo, M., 2012. Transfer of process monitoring models between different plants using latent variable techniques. Ind. Eng. Chem. Res. 51, 7327–7339. Garcìa-Muñoz, S., MacGregor, J.F., Kourti, T., 2005. Product transfer between sites using Joint-Y PLS. Chemom. Intell. Lab. Syst. 79, 101–114. Höskuldsson, A., 1988. PLS regression methods. J. Chemom. 2, 211–228. Hotelling, H., 1933. Analysis of a complex of statistical variables into principal components. J. Educ. Psychol. 24, 417. Jackson, J.E., 1991. A user’s guide to principal analysis. John Wiley & Sons Inc, New York (U.S.A.). Jaeckle, C.M., MacGregor, J.F., 2000. Industrial applications of product design through the inversion of latent variable models. Chemom. Intell. Lab. Syst. 50, 199–210. Jaeckle, C.M., MacGregor, J.F., 1998. Product design through multivariate statistical analysis of process data. AIChE J. 44, 1105–1118. Liu, Z., Bruwer, M.-J., MacGregor, J.F., Rathore, S.S., Reed, D.E., Champagne, M.J., 2011. Scale-up of a pharmaceutical roller compaction process using a joint-Y partial least squares model. Ind. Eng. Chem. Res. 50, 10696–10706. Molino Quaglia, 2016. Personal communication. http://www.molinoquaglia.org/HOME/. Montgomery, D.C., 2013. Design and analysis of experiments. 8th ed., John Wiley & Sons. Moteki, Y., Arai, Y., 1986. Operation planning and quality design of a polymer process. IFAC DYCORD 159–165. Muteki, K., Yamamoto, K., Reid, G.L., Krishnan, M., 2011. De-risking scale-up of a high shear wet granulation process using latent variable modeling and near-infrared spectroscopy. J. Pharm. Innov. 6, 142–156. Rinnan, A.A., van den Berg, F., Engelsen, S.B., 2009. Review of the most common preprocessing techniques for near-infrared spectra. TrAC Trends Anal. Chem. 28, 1201– 1222. _________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)

29

Tomba, E., Barolo, M., García-Muñoz, S., 2012. General framework for latent variable model inversion for the design and manufacturing of new products. Ind. Eng. Chem. Res. 51, 12886–12900. Tomba, E., Meneghetti, N., Facco, P., Zelenková, T., Barresi, A.A., Marchisio, D.L., Bezzo, F., Barolo, M., 2014. Transfer of a nanoparticle product between different mixers using latent variable model inversion. AIChE J. 60, 123–135. Valle, S., Li, W., Qin, S.J., 1999. Selection of the number of principal components: the variance of the reconstruction error criterion with a comparison to other methods. Ind. Eng. Chem. Res. 38, 4389–4401. van den Berg, R.A., Hoefsloot, H.C., Westerhuis, J.A., Smilde, A.K., van der Werf, M.J., 2006. Centering, scaling, and transformations: improving the biological information content of metabolomics data. BMC Genomics 7, 142–157. Wise, B.M., Gallagher, N.B., Bro, R., Shaver, J.M., Windig, W., Koch, R.S., 2006. PLS_Toolbox Version 4.0 for use with MATLABTM. Eigenvector Research Inc., Wenatchee, WA (U.S.A.). Zhou, Y.H., Titchener-Hooker, N.J., 1999. Visualizing integrated bioprocess designs through “windows of operation.” Biotechnol. Bioeng. 65, 550–557.

_________________________________________________________ © 2017 F. Dal-Pastro, P. Facco, E. Zamprogna, F. Bezzo, M. Barolo University of Padova (Italy) and Bühler AG (Switzerland)