Nonlinear model predictive control of a batch fluidized bed dryer for pharmaceutical particles

Control Engineering Practice 64 (2017) 88–101

Contents lists available at ScienceDirect

Control Engineering Practice journal homepage: www.elsevier.com/locate/conengprac

Nonlinear model predictive control of a batch fluidized bed dryer for pharmaceutical particles

MARK

⁎

Francis Gagnona, André Desbiensa, , Éric Poulina, Pierre-Philippe Lapointe-Garantb, JeanSébastien Simardb a b

LOOP, Department of Electrical and Computer Engineering, Université Laval, Québec City, Québec, Canada G1V 0A6 Process Analytical Science Group, Global Manufacturing Services, Pfizer, Montréal, Québec, Canada H4R 1J6

A R T I C L E I N F O

A BS T RAC T

Keywords: Nonlinear model predictive control Phenomenological model Moving-horizon estimation Batch fluidized bed dryer Near-infrared spectroscopy

The availability of reliable online moisture content measurements exploiting near-infrared (NIR) spectroscopy and chemometric tools allows the application of online control strategies to a wide range of drying processes in the pharmaceutical industry. In this paper, drying of particles with a pilot-scale batch fluidized bed dryer (FBD) is studied using a in-line NIR probe. A consolidated phenomenological state-space model of an FBD based on mass and energy balances is calibrated applying a nonlinear least-square identification to experimental data (grey-box modeling). Then, relying on the calibrated model, a nonlinear model predictive controller and a moving horizon state estimator are designed. The objective is to reach a specific particle moisture content setpoint at the end of the drying batch while decreasing cycle time and limiting particle temperature. A penalty term on the energy consumption can also be added to the usual tracking control cost function. Compared to a typical FBD operation in industry (mostly open-loop), it is shown that the drying time and the energy consumption can be efficiently managed on the pilot-scale process while limiting various operation problems like under drying, over drying, or particles overheating.

1. Introduction Because of the difficulty to get some reliable in-line measurements and the regulatory environment in the pharmaceutical industry, manufacturing in this field is not deploying automatic and real-time optimization systems as fast as other sectors (US Food & Drug Administration, 2004; McKenzie, Kiang, Tom, Rubin, & Futran, 2006). The implementation of in-line process sensors would certainly contribute to a better understanding of the systems. Afterward, model predictive control and real-time optimization would surely improve quality control and performance of the manufacturing processes. The wet granulation process, which involves powder wetting and drying steps to create appropriate granules, is an interesting candidate for those improvements. The fluidization of solid particles is widely used in the pharmaceutical industry for drying wet granules. Because of high thermal diffusivities and good mixing properties, the batch fluidized bed dryer (FBD) is an efficient option in terms of cycle time and drying uniformity (Kunii & Levenspiel, 1991). However, in most industrial applications, they are still operated without feedback using fixed batch times and inlet air conditions (Mujumdar, 2014). This can lead to

⁎

Corresponding author. E-mail address: [email protected] (A. Desbiens).

http://dx.doi.org/10.1016/j.conengprac.2017.04.009 Received 17 August 2016; Received in revised form 13 April 2017; Accepted 15 April 2017 0967-0661/ © 2017 Elsevier Ltd. All rights reserved.

multiple problems like under drying or over drying and particle overheating. During the batch, particle overheating can degrade products that are sensitive to high temperatures. Also, a hot product at the end forces operators to wait for a cooling period before manipulating the dried granules. Besides those problems, it may not be the fastest and the most energy-efficient solution. All those issues raised interests in FBD monitoring with in-line analytical instruments and advanced process control strategies (Briens & Bojorra, 2010). Since most batch processes are highly nonlinear due to their usual lack of steady states and their common unidirectional behavior (Bonvin, Srinivasan, and & Hunkeler, 2006; Nagy & Braatz, 2003), nonlinear models and control strategies are more appropriate. More specifically for batch fluidized bed drying, the way the manipulated variables (inlet air flow rate and temperature) are varied influences the rate of drying but not the final particle moisture content, which is fixed by the inlet air humidity. Also, following positive or negative variations of the manipulated variables and from its initial value, the particle moisture content can only decrease until the final equilibrium point is reached. Finally, FBD exhibits two different drying dynamics through the batch defined as the constant rate and the falling rate periods (Li & Duncan, 2008a). The dynamics of the particle temperature are also very distinct during these two phases.

Control Engineering Practice 64 (2017) 88–101

F. Gagnon et al.

z

Nomenclature Af Ah aw cg cp cw cwv Dc Dg dp db0 dbm db Dg0 ef eh g Hc Hf Hp hs hw Hbc

Hbe Hce

Hmf k Kbc Kbe Kce Lor N Nor Nu pe Pw Pr Rwv Re s T t ts Tw Tamb Tref W xi

elevation, m

Greek symbols

inlet blower speed, % inlet heating power, % specific heat-transfer surface of dryer wall, m−1 specific heat of drying gas, J kg−1 K−1 specific heat of dry particles, J kg−1 K−1 specific heat of liquid water, J kg−1 K−1 specific heat of water vapor, J kg−1 K−1 diameter of bed column, m molecular diffusion coefficient of drying gas, m2 s−1 particle diameter, m bubble minimal diameter, m bubble maximal diameter, m bubble diameter, m molecular diffusion coefficient at Tref, m2 s−1 measured inlet airflow error, m3/h measured inlet temperature error, °C gravitational acceleration, m s−2 predictive controller control horizon in multiple of ts, dimensionless expanded bed height, m predictive controller prediction horizon in multiple of ts, dimensionless heat transfer coefficient between drying gas and solids, J s−1 m−3 K−1 heat transfer coefficient between drying gas and dryer wall, J s−1 m−3 K−1 volumetric heat-transfer coefficient between bubble and cloud-wake regions based on volume of bubbles, J s−1 m−3 K−1 volumetric heat-transfer coefficient between bubbles and emulsion based on volume of bubbles, J s−1 m−3 K−1 volumetric heat-transfer coefficient between cloud-wake region and emulsion based on volume of bubbles, J s−1 m−3 K−1 bed height at minimum fluidizing condition, m discrete time, dimensionless coefficient of gas interchange between bubble and cloudwake regions based on volume of bubbles, s−1 coefficient of gas interchange between bubbles and emulsion based on volume of bubbles, s−1 coefficient of gas interchange between cloud-wake region and emulsion based on volume of bubbles, s−1 center-to-center spacing between adjacent holes on perforated plate distributor, m observer window-size or horizon, dimensionless number of orifices per unit area on perforated distributor, m−1 Nusselt number, dimensionless fluid bed dryer energy index weight in controller criterion, J−1 pressure of saturated water vapor, mmHg Prandtl number, dimensionless water vapor specific gas constant, J kg−1 K−1 Reynolds number for packed beds, dimensionless Laplace variable, s−1 temperature, K time, s control algorithm sampling time, s dryer-wall temperature, K fluidized bed dryer room temperature, K reference temperature, K fluid bed dryer energy consumption, J χpm output disturbance state, %w.b.

αχp

relative tolerance on particle moisture content setpoint for

χ χp*

end-of-batch detection, dimensionless moisture content (dry basis), dimensionless moisture content of drying gas on surface of a particle (dry

χpc δb γ κ μg ν ϕ ψ1 ψ2 ρg ρds ρs ρws ρw σ υ υbr υmf εmf ε0

basis), dimensionless particle critical moisture content (dry basis), dimensionless fraction of fluidized bed consisting of bubbles, dimensionless heat of vaporization at Tref, J kg−1 offset constant in ψ2 function dynamic viscosity of drying gas, kg m−1 s−1 exponent constant in ψ2 function sphericity of particles, dimensionless gas saturated humid curve moisture film correction function based on the absorption characteristics of the particle density of gas, kg m−3 dry bulk density of the bed, kg m−3 dry particle density, kg m−3 wet bulk density of the bed, kg m−3 water density, kg m−3 evaporation coefficient, kg m−2 s−1 gas superficial velocity based on the total cross-sectional area of bed, m s−1 linear velocity of a single bubble, m s−1 gas superficial velocity at minimum fluidizing condition, m s−1 void fraction at minimum fluidizing conditions, dimensionless inlet air relative humidity, %RH

Vectors and matrices

Q R r u xA x ys

observer a priori state estimates (4 × 1) observer estimated augmented states (4 × 1) measured disturbances (1 × 1) augmented model state update function (4 × 1) state-space model state update function (3 × 1) augmented model state output function (3 × 1) state-space model output function (3 × 1) setpoints tracking weights (3 × 3) manipulated input weights (2 × 2) horizon initial observation error estimated covariance (4 × 4) process noise estimated covariance (4 × 4) observation noise estimated covariance (3 × 3) setpoints (3 × 1) manipulated inputs (2 × 1) augmented model states (4 × 1) phenomenological model states (3 × 1) phenomenological model outputs (3 × 1)

y

process measured outputs (3 × 1)

x0 xl d f A(•) f(•) hA(•) h(•) M N P

Subscripts 0 b e m p

89

inlet gas bubble phase emulsion phase: interstitial gas measured emulsion phase: solid particles

Control Engineering Practice 64 (2017) 88–101

F. Gagnon et al.

sp

MPC NGBM NIR NMPC ODE PI PLS RK4 RTD SQP USB

setpoint

Acronyms EKF FBD GMC HMI LOD MHE MISO

extended Kalman filter fluidized bed dryer generic model control human-machine interface loss-on-drying moving-horizon estimator multiple-input single-output

model predictive control nonlinear grey-box model near-infrared nonlinear model predictive control ordinary differential equations proportional-integral partial least squares Runge-Kutta method order 4 resistance temperature detector sequential quadratic programming universal serial bus

desired particle moisture content while limiting particle temperature and weighting energy consumption. A nonlinear moving-horizon estimator (MHE) evaluates the current model states and the unmeasured moisture disturbance (to ensure a zero static error). The approach is tested on both the simulated process and the pilot-scale system.

Many phenomenological models of FBDs have been developed in the past, based on various assumptions. Generally, a more accurate model requires more computation time. If the model is used in a predictive control strategy, it may become necessary to handle this trade-off to ensure an adequate control period. One of the most widespread FBD lumped models is the KuniiLevenspiel (Kunii & Levenspiel, 1991) model. It assumes three distinct phases that exchange mass and energy. The emulsion phase is composed of perfectly mixed solid particles and interstitial gas. The second phase is the bubbles that go up through the emulsion phase. The third phase, the clouds and wakes around the bubbles, is often neglected. One of the drawbacks of this model is the assumption that the moisture content and the temperature for each phase are constant along the height z in the powder bed. This simplification is not realistic for taller fluid beds like the ones used for pharmaceutical tablet production. To counteract this problem, Li & Duncan (2008a) improved the Kunii-Levenspiel two-phase model by including a z-dependency in the balance equations of the bubble phase. Villegas, Duncan, Wang, Yang, and Raghavan (2010) later compared the same model to experimental batches with an instrumented lab-scale fluid bed. The same assumptions are made and experimentally validated on pharmaceutical granules by Wang, Dyakowski, Senior, Raghavan, and Wang (2007). Relying on the modified two-phase model and nonlinear programming algorithms, Villegas et al. (2010) demonstrated that it is possible to reduce the energy consumption on a lab-scale FBD by pre-calculating optimal inlet airflow and temperature profiles for the entire batch. With the same phenomenological model, Li & Duncan (2008b) showed that an extended Kalman filter (EKF) can act as a virtual sensor for the particle moisture content. The virtual measurement combined with a generic model control (GMC) strategy allowed the control of the particle moisture content of a simulated disturbed process. Both of these strategies may lead to static errors in presence of disturbances or modeling mismatches as they do not explicitly incorporate an integral action in the control algorithm. Obregón, Quiñones, and Velázquez (2013) later used a linear model predictive control (MPC) approach to successfully control particle moisture content on a lab-scale FBD. In-line near-infrared (NIR) spectra were obtained with a NIR probe and moisture content measurements were achieved by using a pre-calibrated partial least squares (PLS) model (see Demers, Gosselin, Simard, & Abatzoglou, 2012, for a similar approach). Based on a simple multiple-input singleoutput (MISO) first-order process model, this approach cannot constrain particle temperature (the only output being the moisture content). The energy consumption was not explicitly part of the cost function. It may also be probable that control performances are not optimal because the nonlinear FBD behavior is not accurately represented. In this work, a consolidated phenomenological model of an instrumented pilot-scale FBD is calibrated by nonlinear grey-box model identification. It then serves as the predictive model in a nonlinear model predictive controller (NMPC). A NIR probe provides a particle moisture content measurement. The objective is to reach a

2. Materials and methods To generate wet granulation process data, a raw material that commonly appears in drug formulations is processed in a pilot-scale granulator and fluidized bed dryer. The following subsections describe this methodology in detail. 2.1. Materials For each drying batch, 5.3 kg of calcium carbonate powder is granulated with two excipients in a high-shear wet granulator with 1.19 kg of purified water at 55–60 °C. To perform wet granulation, dry ingredients are first mixed by the impeller of the granulator. This is then followed by water addition and appropriate dispersion. A final wet granulation high shear step is performed with chopper and impeller during at least two minutes. The wet granules are then sieved with a 1.18 mm mesh before proceeding with the drying. The wet batch weight put in the FBD is 5.7 kg . 2.2. Equipment The pilot-scale FBD is an Aeromatic MP-1 (max capacity of 25 L) already equipped with various sensors and modified to insert multiple process analytical instrument probes. A Viavi Solutions (formerly JDSU) MicroNIR Pro near-infrared spectrometer mounted in a custom-made casing is inserted in the fluid bed column to acquire spectra during the drying process. All other measurements and manipulated signals are processed through a National Instrument USB-6343 acquisition card:

• • • • • •

inlet air temperature, with a resistance temperature detector (RTD); inlet air flow rate, with a Pitot tube; inlet air relative humidity, with a capacitive sensor; particle temperature, with a RTD; inlet heating power, with a hot water heat exchanger; and inlet blower speed, with a motor drive.

Off-line measurements of particle moisture content are obtained by a loss-on-drying (LOD) method with a Sartorius moisture analyzer. 2.3. Software A human-machine interface (HMI) and data logger are programed with National Instrument LabVIEW to quickly visualize all relevant FBD signals. The same software allows measurements filtering and implementation of local proportional-integral (PI) control loops (with a 90

Control Engineering Practice 64 (2017) 88–101

F. Gagnon et al.

Table 1 Inlet heater process PI gain-scheduling.

υ0m

50

Gh(s ) =

T0m(s ) Ah (s )

1.03exp( − 0.0 s ) 119.5 s + 1

m3 h

100

m3 h

1.09exp( − 0.5 s ) 103.0 s + 1

200

m3 h

1.19exp( − 8.6 s ) 52.02 s + 1

300

m3 h

1.13exp( − 8.8 s ) 33.89 s + 1

350

m3 h

1.12exp( − 7.9 s ) 28.15 s + 1

Ch(s ) =

Ah (s ) eh(s )

⎛ 1 ⎞ 1.66⎜1 + ⎟ ⎝ 120 s ⎠ ⎛ 1 ⎞ 1.35⎜1 + ⎟ ⎝ 103 s ⎠ ⎛ 1 ⎞ 0.72⎜1 + ⎟ ⎝ 60.0 s ⎠

⎛ 1 ⎞ 0.53⎜1 + ⎟ ⎝ 42.0 s ⎠ ⎛ 1 ⎞ 0.45⎜1 + ⎟ ⎝ 35.0 s ⎠

and the NMPC manipulated variables (PI loop setpoints): Fig. 1. FBD control strategy.

T u = [T0sp υ0sp ]

sampling period of 1 s). The NMPC and the MHE are coded in MATLAB. They both require numerical optimization (MATLAB routine fmincon) with multiple executions of the FBD model. Thus, to accelerate the calculations, the FBD model is compiled into a mex file. At each NMPC sampling period (1 min), the controller and observer are called by LabVIEW HMI with the ”MATLAB script” LabVIEW module.

The NMPC also exploits the availability of the inlet relative humidity (%RH) as a measured disturbance:

d = [ φ0 ]T

(3)

The NMPC setpoint vector is denoted:

⎡ ⎤T r = ⎢* χpsp *⎥ ⎣ ⎦

2.4. NIR probe

(4)

Only the particle moisture content setpoint χpsp is considered in this paper since the other measurements only need constraining (zero weights are assigned to the other setpoints).

In-line NIR spectra are read by selecting the auto-saving ”timed acquisition mode” of the Viavi MicroNIR software. Hence, at each second, a new spectrum is automatically saved as a comma-separated values (csv) file in a specific folder. Meanwhile, from LabVIEW, the most recent csv file is found at each second, opened and logged. If moisture content measurement is desired, a statistical model is applied to the spectrum which is first filtered with the Savitzky-Golay algorithm, with a window-size of 15 and a polynomial order of 2 (Demers et al., 2012). For a reliable moisture content measurement, the statistical model, which is a partial-least square (PLS) model, must be pre-calibrated with calcium carbonate particles (Demers et al., 2012). The PLS model is calculated with Umetrics SIMCA using the nonlinear iterative partial least square (NIPALS) algorithm (Wold, 1975) from five FBD batches of the same raw material, granulated the same way. At every two minutes, granules are sampled from the FBD with a powder thief and moisture content is measured with the LOD analyzer (the PLS Y -matrix). As data pretreatment, the same Savitzky-Golay algorithm is applied to the full NIR spectrum and bundled into the PLS X -matrix. The resulting regression coefficients and Y-intercept (unscaled values) are exported from SIMCA and used in the LabVIEW HMI for moisture content measurement. Ten or even more batches are preferable but the model was precise enough for our purposes since both the root mean square error of estimation and the root mean square error of cross validation were around 0.62% in wet basis.

2.6. PI control loops Inlet blower and heater models are obtained by system identification using experimental data collected on the pilot FBD. Models with different structures are estimated. The models are cross-validated (i.e. with data not used for identification purposes) with techniques such as the analysis of the covariance matrix of the estimated parameters, the autocorrelation and cross-correlation of the residuals with the input, the comparison of the predicted model output to the measured output, etc. (Ljung, 1999). First or second order models are selected since they provide a good trade-off between simplicity and accuracy. The sampling period for their associated PI controller is fixed at 1 s. 2.6.1. Inlet blower identification and control The identified second-order model is:

Gf (s ) =

υ0m(s ) 5.4 = A f (s ) (12.3 s + 1)(11.6 s + 1)

(5)

Based on the model, a controller is designed to obtain a 60° phase margin and a 10 − 90% rise-time of 20 s. The resulting PI compensator is:

Cf (s ) = 2.5. Overview of the whole control strategy Fig. 1 depicts the control strategy. There are two local PI loops. They respectively manipulate the inlet heating power Ah (%) and the inlet blower speed Af (%) to control the inlet temperature measurement T0m (°C) and airflow measurement υ0m (m3/h). The MHE provides the required state estimates xl to the NMPC from the inlet temperature (°C), particle moisture content (% wet basis) and particle temperature (°C) measurements:

A f (s ) ef (s )

⎛ 1 ⎞ = 0.222⎜1 + ⎟ ⎝ 16.2 s ⎠

(6)

2.6.2. Inlet heater identification and control The inlet heater dynamics significantly depend on the air flow rate. Therefore, five different first-order plus dead-time models are identified on the pilot FBD, as shown in the first two columns of Table 1. For each model, the third column gives the corresponding PI controller that leads to a closed-loop time constant of approximately 70 s. The gainscheduled PI is implemented by linearly interpolating the proportional gain and integral time from the inlet airflow measurement υ0 m as the scheduling variable.

T

y = [T0m χpm Tpm ]

(2)

(1) 91

Control Engineering Practice 64 (2017) 88–101

F. Gagnon et al.

Fig. 2. Local PI closed-loop step responses.

and the solid particles. 3. The emulsion phase remains in minimum fluidization conditions. The excess airflow above minimum fluidization velocity passes through the bed as bubbles in the form of a plug flow. 4. The bubble phase contains no particles and the clouds surrounding them have zero thicknesses. Therefore, the bubbles exchange mass and energy only with the interstitial gas. The gas inside the bubble is perfectly mixed. Bubbles moisture content, temperature and size depend only on their position in z along the fluid bed. 5. All changes in physical properties of the gas and particles due to the change of temperature are neglected (except water diffusion coefficients, see Wang et al., 2007). Particle size remains constant during the drying process so it can be assumed that the total bulk volume of the bed is the one occupied with the particle at the critical moisture content χpc.

2.6.3. PI control results Inlet heater and blower PI loops are tested with several step responses, as shown in Fig. 2. Vertical dotted lines indicate bumps on the airflow setpoint υ0sp , hence large variations in the scheduling variable υ0m are about to occur. The loops behave as expected. Since the closed-loop model of the heater is required in the next section, it is identified from the data plotted in the first graph of Fig. 2:

T0m(s ) 1 = T0sp(s ) 71.2 s + 1

(7)

The dead time is neglected since it is considerably much shorter than the time constant. The time response of the second PI loop being much shorter than all other dynamics in the system, the model is simply approximated as follows:

υ0m(s ) ≈1 υ0sp(s )

The model mass and energy balance equations are presented in the following sub-sections. The remaining equations of the model are given in Appendix A and symbol definitions are listed in the Nomenclature ∼ section. In all equations, f is defined as the spatial average of the function f(z) along the bed height:

(8)

3. FBD model description 3.1. Assumptions

1 ∼ f ≜ Hf

The batch FBD model from Li & Duncan (2008a) is consolidated relying on an earlier work from Lai, Chen, and Fan (1986) describing this same model applied on the continuous dryer. Modeling elements from Wang et al. (2007) related to particle and gas physical properties are also added to include more versatility on simulation conditions and to simplify calculations. The lumped capacitance model assumes two distinct phases (Kunii & Levenspiel, 1991): a bubble phase and a perfect emulsion phase of the solid particles and interstitial gas. The main assumptions are:

∫0

Hf

f (z )dz (9)

3.2. Balance equations for the bubble phase From Li & Duncan (2008a) and Lai et al. (1986), the bubble phase mass balance analytical solution with the proper boundary conditions is:

⎛ K δ ⎞ χb (z ) = χe − (χe − χ0 )exp⎜ − be b z⎟ υb ⎠ ⎝

1. All particles in the fluid bed are identical. They all have the same physical properties, moisture content and temperature at a given time during the drying. The moisture content and temperature inside each particle are uniform. 2. The emulsion phase is perfectly homogeneous, i.e. the solid particles are perfectly mixed with the interstitial gas. The interstitial gas mediates mass and thermal exchanges between the bubble phase

(10)

From the simplification suggested by Li & Duncan (2008a) about the enthalpy of the water vapor in the emulsion gas and with the approximation that χb (z ) ≈ χ͠ b , the bubble phase energy balance analytical solution with the proper boundary conditions becomes:

⎛ ⎞ Hbeδb Tb(z ) = Te − (Te − T0 )exp⎜⎜ − z⎟⎟ ρ ( c + χ c ) υ ͠ ⎝ g g b wv b ⎠ 92

(11)

Control Engineering Practice 64 (2017) 88–101

F. Gagnon et al.

υ0 = 3.3. Balance equations for the emulsion phase: interstitial gas From Li & Duncan (2008a) and Lai et al. (1986) and by adding the particle sphericity ϕ as proposed by Wang et al. (2007), the interstitial gas mass and energy balance equations are respectively:

⎛ 6 ⎞ ⎞ ⎛ ⎟⎟(1−εmf )(1−δb )σ ⎜χp* −χe ⎟ ρg (χe −χ0 ) = ρg δbKbe(χ͠ b −χe ) + ⎜⎜ ⎝ Hf ⎠ ⎝ dpϕ ⎠

υ0sp 900πDc2

(19)

3.7. State-space representation The state vector of the FBD model described by equations (10)–(19) is:

υmf

T

(12)

x(t ) = [T0(t ) χp (t ) Tp(t )]

and

The system of ordinary differential equations (ODE) given by the time derivative of eq. (20) are solved with the fourth order Runge-Kutta method (RK4) with a sampling period of 1 min. The resulting discrete nonlinear state-space equations are:

⎛ 6 ⎞ υmf ⎟⎟(Tp − Te ) ρg (cg + cwvχ0 )(Te − T0 ) = (1 − εmf )(1 − δb )⎜⎜ Hf ⎝ dpϕ ⎠ ⎡ ⎤ × ⎢σ ( χp* − χe )cwv + hs ⎥ + δbHbe(T͠ b − Te ) ⎣ ⎦ + awh w(Tw − Te )

(20)

(13)

x(k + 1) = f(x(k ), u(k ), d(k ))

(21a)

ys(k ) = h(x(k ))

(21b)

with f(x , u, d) given by the RK4 solver and h(x) given by:

⎡ T0 − 273.15 ⎤ ⎢ ⎥ ⎞−1⎥ ⎢ ⎛ h( x ) = ⎢100⎜χp−1 + 1⎟ ⎥ ⎠ ⎥ ⎢ ⎝ ⎢⎣ Tp − 273.15 ⎥⎦

3.4. Balance equations for the emulsion phase: solid particles Again, from Li & Duncan (2008a) and Lai et al. (1986), and by taking into account the particle sphericity (Wang et al., 2007), the solid particle mass and energy balances are respectively:

ρds

dχp dt

=−

6 σ ( χp* − χe ) ϕdp

As shown in eq. (21b) and (22), the model outputs ys are the inlet temperature in °C , the particle moisture content in % wet basis and the particle temperature in °C . The units are the same as the measured outputs y .

(14)

and

ρds (cp + χp cw )

dTp dt

=

(22)

6 ⎡ ⎢hs(Te − Tp ) − σ (χp* − χe ) ϕdp ⎣

3.8. Energy consumption

⎤ × [cwv(Te − Tref ) − cw(Tp − Tref ) + γ ]⎥ ⎦

As proposed by Villegas et al. (2010), who also worked with a twophase model, the FBD energy consumption is:

(15)

To simplify the model, the reference temperature is fixed at room temperature:

Tref = Tamb

W = cgρg

πDc2 4

∫0

tf

υ0(T0 − Tamb )dt

(23)

The energy consumed by the blower is neglected and a perfectly efficient inlet air heater is assumed.

(16)

3.5. Inlet air relative humidity

4. FBD model calibration

In equations (10)–(15), the inlet gas moisture content χ0 is the absolute humidity in dry basis. Since the inlet gas moisture content measurement is a relative humidity percentage, a conversion between absolute humidity χ0 and relative humidity φ0 is required. Knowing that the wet basis moisture content is (χ −1 + 1)−1 and from the ideal gas law, if the measure is taken at temperature Tamb, the relation is:

⎛ ⎞−1 ρg R wvTamb ⎜ ⎟ χ0 = ⎜ φ0 − 1⎟ ⎜ 133.32Pw(Tamb ) ⎟ ⎝ ⎠ 100

4.1. Calibration method In order to use the batch FBD model described in Section 3 for the design of the NMPC and MHE, physical, chemical and empirical Table 2 Phenomenological model calibrated parameters.

(17)

3.6. Inlet air heater and blower closed-loop dynamics Two of the inputs of the FBD model described in sections 3.2 to 3.5 are T0 (in K) and v0 (in m s−1). Because of the cascade hierarchy, the links between them and the real manipulated variable are eq. (7) and eq. (8). Equation (7) corresponds to the following expression when converting the setpoint from °C to K:

(T0sp + 273.15) − T0 dT0 = dt 71.2

(18) 3

When converting the inlet air velocity setpoint units from m /h to m s−1, eq. (8) becomes: 93

(a) Fixed parameters

(b) Free parameters

par.

value

par.

value ± std. dev.

Dc Hmf Lor cg cw cwv ρg ρs ρw kg g Rwv Tamb κ

0.30 m 0.20 m 0.005 m 1.0/(kg K) 4.2/(kg K) 1.9/(kg K) 1.1 kg/m2 2800 kg/m2 1000 kg/m2 29 mW/(m K) 9.8 m/s2 462 J/(kg K) 293 K 0.45

cp dp ϕ ν χpc μg Dg0 Tw

(3.9 ± 0.04) kJ/(kg K) (1.10 ± 0.01) μm (0.653 ± 0.0006) (1.11 ± 0.03) (0.027 ± 0.0005) (71 ± 0.8) mg/(m s) (21 ± 0.3) μm2/s (297 ± 0.2) K

Control Engineering Practice 64 (2017) 88–101

F. Gagnon et al.

Fig. 3. Particle sphericity sensitivity analysis, ϕ = 0.6 ± 0.12 .

parameters must be calibrated. An empirical calibration method based on nonlinear optimization is selected. This approach allows to optimize the fitting of the predicted model outputs, which is highly desired for control purposes. More precisely, nonlinear grey-box model identification (NGBM) functions from the MATLAB identification toolbox are employed. The calibration relies on a least-square algorithm method available in the MATLAB optimization toolbox. NGBM tools include features like parameter and state constraints, fixed parameters and initial state values, measured output weighting, variable-step ordinary differential equation (ODE) solvers and model validation tools. 4.2. Selection of optimized parameters Many parameters of the model do not need to be estimated by NGBM identification as they are either known or can easily be measured. A sensitivity analysis is made for all remaining parameters. If the model behavior is not sensitive to a parameter, a reasonable expected value is simply assigned. The parameters with a fixed value are given in Table 2, and their corresponding description are listed in the Nomenclature section. The parameter sensitivity analysis is similar to what was described by Li and Duncan (2008a). At first, all the simulator parameters are set near their theoretical or expected value. A drying batch simulation of 40 min with constant manipulated inputs and measured disturbance is executed as a reference. Then, for each parameter, two other cases are simulated: by increasing and decreasing the parameter value by 20%. The resulting particle moisture content and temperature are then compared in a graph. An example of a high-sensitivity parameter (the particles sphericity) is given in Fig. 3. 4.3. Measurements, initial state and output weights For both model calibration and validation data, T0m, Tpm, υ0m and φ0 are measured at each second with the probes already available on the pilot fluid bed. The particle moisture content measurement is obtained by one of the two following methods: 1. For χpm > 5% w.b., LOD scale is used by sampling powder from the fluid bed every two minutes approximately. Spline interpolation is then used to generate a moisture content measurement at each second. 2. For χpm ≤ 5% w.b., the NIR probe with the PLS model provides a measurement at each second. With the actual calibrated PLS model, particle moisture predictions are more precise for values smaller than 5%, being very similar to LOD measurements, which are the reference to calibrate the PLS model.

Fig. 4. Model validation.

with constraints at ±0.5 %w.b. of the measured value. The inlet and particle temperature measurements are accurate and thus the corresponding initial state values are selected equal to the measurements. For the calibration procedure, a relative weight must be assigned to each output of the model in the minimized loss function (weighted sum of squared errors). The inlet temperature weight is set to zero since this part of the model has already been identified at Section 2.6. Weights of 1 and 0.7 are given to χpm and Tpm respectively. The highest priority is given to the particle moisture content since it is the controlled output.

Initial values have to be given to the NGBM identification and validation algorithms. Since interpolation is used for large values of χpm, which may be less reliable, the particle initial moisture content is set free 94

Control Engineering Practice 64 (2017) 88–101

F. Gagnon et al.

values for the discrete time i estimated at the discrete time j. It is worth noticing that, even if all states are measured, a state estimator is still recommended since:

4.4. Calibrated model validation The model calibration is performed with one drying batch with various bumps on the inlet temperature setpoint T0sp and airflow setpoint υ0sp . Fig. 4a compares validation data (i.e. data from a second batch, not used for calibration) and the outputs predicted with the model using the signals shown in Fig. 4b as inputs. Data from more batches is recommended but this was sufficient for this proof-ofconcept as the predictions are very good with fits around 90% for all three measured outputs (normalized root-mean-square error). The optimized value of the free parameters and their associated standard deviations are given in the Table 2. Parameter descriptions are provided in the Nomenclature section. At first, the selected optimization variables were constrained to an expected range of realistic values. However, the resulting model was less accurate than the model identified without constraining the parameters. Since the predicting capability of the model is by far the most important criterion when selecting a model for MPC design, it was decided to select the second approach even if some calibrated parameters (for instance particle diameter and viscosity of the gas) converged to unrealistic values. As stated by Ljung (1999), the identification procedure provides us a particular model for the selected structure model and ”the crucial question is whether the identified model is good enough for the intended purpose”, which is here the adequate control of the FBD. In this project, the model itself is not a goal but a tool. The standard deviations of the estimated parameters being relatively small, the convergence of some parameters outside their expected range of values may be caused by all the limiting assumptions of the model. For instance, by assuming a constant particle size dp, particle attrition and fine generation phenomena are neglected, which are known to be significant in most fluid bed applications (Hao, Zhao, Ye, & Lius, 2016). Moreover, appreciable changes in the dryer wall temperature was noted during the batches implying that a constant wall temperature assumption is not realistic.

• • •

the unmeasured moisture disturbance estimation incorporates an integral action into the controller; it filters out stochastic noises by relatively weighting the model predictions and the measurements (with the respective covariance matrices); and it can act as a soft-sensor in case of a sensor failure. For instance, the NIR moisture measurement is reliable only if the probe window is clear of powder accumulation and if the particle fluidization is good enough. In the worst case, there is no measurement at all. The state observer is then used as a backup virtual sensor.

Also, we can take advantage of the filtering effect for an end-ofbatch detection. More precisely, drying is terminated when the estimated moisture content value is close enough to the desired setpoint (see Section 5.4). The selected observer algorithm is a moving-horizon state estimator, which is an optimization based estimator described by Rao and Rawling (2002). The main advantages of this observer are the handling of nonlinear dynamics (without linearization) and estimated state constraints. To avoid the optimization problem to grow indefinitely with time, the observation window is truncated at N data points in the past. However, the calculations are much heavier than an extended Kalman filter. At each time-step k, the following MHE optimization problem is solved:

min

m xlk (k − N ), W

Jmhe

(26)

where the criterion is:

Jmhe = [xlk(k − N ) − xk]T P−1[xlk(k − N ) − xk] m TQ−1W m+V l TR−1V l +W

5. Control algorithms

(27)

based on the disturbed model:

5.1. Augmented states model In real-life control problems, modeling mismatches and disturbances are unavoidable. Thus, it is necessary to add an integral action in the predictive control strategy to eliminate steady-state errors. One way of doing that is by estimating the output disturbances, represented as a random walk (Maciejowski, 2000). In the present case, a disturbance xi is added only to the particle moisture content since it is the unique controlled variable. The resulting augmented state vector is: (24)

The corresponding augmented model equations are:

⎡ f(x(k ), u(k ), d(k ))⎤ xA(k + 1) = ⎢ ⎥ xi(k ) ⎣ ⎦ = f A( xA(k ), u(k ), d(k ))

(25a)

⎡ 0 ⎤ yA(k ) = h(x(k )) + ⎢ xi(k )⎥ ⎢ ⎥ ⎣ 0 ⎦ = hA(xA(k ))

(25b)

(28a)

y(k ) = hA( xl(k ) ) + vl(k )

(28b)

l k ) and v( l k ) are stochastic processes with their In (28a) and (28b), w( associated covariance matrix Q and R respectively. In (27), the weighting matrices are:

T

xA(k ) = [T0(k ) χp (k ) Tp(k ) xi(k )]

l(k ) xl(k + 1) = f A( xl(k ) , u(k ) , d(k ) ) + w

⎡Q ⎢ 0 Q=⎢ ⎢⋮ ⎢⎣ 0

0 Q ⋮ 0

… … ⋱ …

0⎤ ⎥ 0⎥ ⋮⎥ Q ⎥⎦ 4N ×4N

(29)

⎡R ⎢ R = ⎢0 ⎢⋮ ⎣0

0 R ⋮ 0

… … ⋱ …

0⎤ 0⎥ ⎥ ⋮⎥ R ⎦3N ×3N

(30)

The estimated process noise over the observation horizon is:

⎡w l(k ⎢ l(k w m ⎢ W= ⎢ ⎢⎣

− N + 0)⎤ ⎥ − N + 1)⎥ ⋮ ⎥ ⎥⎦ l(k ) w

(31)

and the estimated observation noise over the observation horizon is:

with f(•) and h(•) being defined by equations (21a) and (21b).

⎡ y(k − N + 0) − hA( xlk(k − N + 0) )⎤ ⎢ ⎥ l = ⎢ y(k − N + 1) − hA( xlk(k − N + 1) )⎥ V ⎢ ⎥ ⋮ ⎢ ⎥ y(k ) − hA( xlk(k ) ) ⎣ ⎦

5.2. State observer algorithm In what follows, the notation xlj(i ) refers to the augmented state 95

(32)

Control Engineering Practice 64 (2017) 88–101

F. Gagnon et al.

⎡ r(k )⎤ ⎢ ⎥ l = ⎢ r(k )⎥ R ⎢ ⋮ ⎥ ⎢⎣ r(k )⎥⎦

At a beginning of a batch and as long as k < N , since there are not enough points to fill the complete observation window, all matrices are correctly truncated to correspond to the appropriate dimensions. In the criterion (27), P is the covariance matrix of the state estimation error between x(k − N ) and the horizon initial value xk , which is:

⎧x0 if k ≤ N xk = ⎨ ⎩ xlk − N −1(k − N ) if k > N

(33)

∀ i ∈ [ − 1 .. min{k , N}]

(34)

⎡N ⎢ N = ⎢0 ⎢⋮ ⎣0

0 N ⋮ 0

… … ⋱ …

… … ⋱ …

0⎤ 0⎥ ⎥ ⋮⎥ M ⎦3H ×3H p

p

0⎤ 0⎥ ⎥ ⋮⎥ N ⎦2H ×2H c

c

(35)

(43)

cgρg ts 3600

Hp −1

∑

(Tlk −1(k + i ) − Tamb )υ0sp(k + i )

i =0

(44)

where Tlk−1(k + i ) is the first element of xlk−1(k + i ) and υ0sp(k + i ) is the second element of u(k + i ). Predictive controller criterion. At each time-step k, the complete NMPC optimization problem is thus:

min Jtrk + pe Jw

(45)

ΔU

where Jtrk is given by (36) and Jw by (44) while pe weights the energy index. The minimization can be calculated with respect to the following constraints:

Δumin < Δu(k + i ) < Δumax

(36)

∀ i ∈ [0 .. Hc[

(46)

u min < u(k + i ) < u max

∀ i ∈ [0 .. Hc[

(47)

ylmin < ylk(k + i ) < ylmax

∀ i ∈ [1 .. Hp]

(48)

Notice that it is therefore possible to constrain the particle temperature. According to the receding horizon principle, only the first optimal input increment Δu(k ) is applied to the FBD, and the whole problem is solved again at the following control period. (37) 5.4. Stopping condition Drying is assumed to be completed when the current particle moisture content estimate calculated by MHE is smaller than the moisture content setpoint. Hence, the stopping condition is:

(38)

[ 0 1 0 ] ylk(k ) ≤ [ 0 αχp 0 ]r(k )

In the cost function (36), the manipulated input increments over the control horizon Hc are:

⎡ Δu(k + 0) ⎤ ⎢ ⎥ Δu(k + 1) ⎥ ΔU = ⎢ ⋮ ⎢ ⎥ ⎢⎣ Δu(k + Hc − 1)⎥⎦

(42b)

Jw =

where the weighting matrices are:

0 M ⋮ 0

ylk(k + i ) = hA(xlk−1(k + i ))

Energy index. The energy consumption is defined by (23). The numerical integration is obtained with a rectangular approximation (Euler's method). The integration is calculated from time k to time k + Hp − 1. In (23), T0 is replaced by its prediction calculated with (42a). Substituting (19), the energy consumption performance index becomes:

In this section, the notation ylj(i ) represents the augmented model output values predicted for the discrete time i at the discrete time j, and thus based on the state estimations available at the discrete time j. Relying on the receding horizon principle, the proposed NMPC leads to the optimal manipulated variables by minimizing at each sampling time k a weighted sum of two criteria. The first one is the usual setpoint tracking cost function. The second criterion is an energy consumption index. Setpoint tracking. The usual model predictive controller criterion is:

⎡M ⎢ M=⎢0 ⎢⋮ ⎣0

(42a)

l (k + i ) = d(k ) d k

5.3. Predictive control algorithm

l−Y l)T M(R l−Y l) + (ΔU)T N(ΔU) Jtrk = (R

l (k + i ) ) xlk −1(k + i + 1) = fA(xlk −1(k + i ) , u(k + i ) , d k

In (42a), all future measured disturbances are simply assumed equal to the current value measurement:

m values are estimated by When the optimal xlk(k − N ) and W minimizing the cost function, the next state estimation xlk(k + 1) can be deduced using (34) with i=0. This estimate is kept in memory and becomes the required state estimate by the NMPC at the next control period. This allows at each control period to first solve the NMPC problem and then the observer problem. The optimal manipulated variables are hence obtained more rapidly than the other way around. Because the MHE observer is optimization-based, constraints on state estimations can be added: xl min < xlk(k − i ) < xl max

(41)

In (41), the output predictions are recursively calculated using the augmented state-space model defined by (25a) and (25b). The initial state vector xlk−1(k ) was estimated by the state observer at the previous control iteration. The prediction equations are:

where x 0 is an a priori state estimate for discrete-time 0. In (32), state estimates are obtained by recursively applying the augmented disturbed model equation:

l(k − i ) xlk(k − i + 1) = f A( xlk(k − i ) , u(k − i ) , d(k − i ) ) + w

and

⎡ yl (k + 1) ⎤ ⎢ k ⎥ l = ⎢ ylk(k + 2) ⎥ Y ⎢ ⎥ ⋮ ⎢ ⎥ l k H y ( + ) p ⎦ ⎣ k

(49)

with αχp fixed at 1.05, giving a relative tolerance. It is worth mentioning that, in this end-of-batch detection, the particle moisture content measurement is indirectly considered through the unmeasured moisture disturbance estimation xli .

(39) 6. Results and discussion

where:

Δu(i ) = u(i ) − u(i − 1)

The control algorithm described in Section 5 is first tested on the calibrated simulator. Multiple batches are simulated to illustrate the effect of the weight pe on the batch time and the energy consumption. Then, three selected closed-loop scenarios are compared: (1) NMPC

(40)

The output setpoints are considered constant over the prediction horizon Hp and the predicted outputs are calculated for k + 1 to k + Hp : 96

Control Engineering Practice 64 (2017) 88–101

F. Gagnon et al.

with the objective of reducing the drying time, and NMPC with a (2) low and (3) high weight on the energy consumption. Finally, a usual open-loop and the second closed-loop scenarios are evaluated with the pilot FBD. The selected NMPC parameters are:

• • • • • • • • • • •

control period: ts = 60 s; Hp = Hc = 25, thus covering the whole batch and giving full flexibility to the manipulated variables (without a weight on the energy index, the drying time should be reduced); r(k ) = [ 0 1 0 ]T ; M = diag{0, 1, 0}; N = 0.001 × diag{1, 0.4}; Δumin = [−20 −50 ]T ; Δumax = [+20 +50 ]T ; u min = [ 40 75]T ; u max = [ 90 200 ]T ; ylmin = [−∞ −∞ −∞]T ; and ylmax = [+∞ +∞ +45]T .

Thus, the desired particle moisture content is 1%w. b.. To smooth the manipulated variable moves, N , Δumin and Δumin are not zero (knowing that this may lightly increase the drying time). Inlet temperature constraints and the maximal air flow rate constraint correspond to FBD physical limits. Minimal air flow rate is limited to ensure proper fluidization. Maximal particle temperature is also constrained. The MHE parameters are:

• • • • • •

Fig. 6. FBD simulator in closed-loop: priority on batch time (no weight on energy consumption).

Wtb =

cgρg ts 3600

kend

∑ (T0m(k ) − Tamb )υ0m(k ) k =0

(50)

6.1. Simulation results

N = 10; R = diag{0.436, 0.0625, 0.109}, from the magnitude of the sensor noises on the pilot FBD; Q = 10−4 × diag{100.0, 0.0025, 1.0, 39.0}; P = 0.01Q ; xl min = [ 0 0 0 −∞]T ; and xl max = [+∞ +∞ +∞ +∞]T .

The inlet relative humidity is kept constant throughout the simulation:

d(k ) = 30 The initial conditions are:

x(0) = [ 323.15 0.1765 298.15]T No disturbances nor modeling mismatches are included in the simulations (even if they have been successfully tested) to ease the comparison of the batch durations and energy consumptions.

Q and P are not easy to select. They were simply tuned from several simulations by using direct state measurements for the a priori state estimate x 0 . A more refined method to select P based on model linearization is described in Rao and Rawling (2002). Both NMPC and MHE problems are solved with the MATLAB function fmincon. The combination of sequential quadratic programming (SQP) method and central finite difference calculation lead to the most reliable and coherent results for both the controller and state estimator. The stopping condition is defined in Section 5.4. The final discrete time is kend. The total energy consumption of a terminated batch is evaluated with:

6.1.1. Closed-loop: energy index weight impact on simulations The effect of the energy index weight pe on the total drying time and energy consumption is investigated on the FBD simulator. Relying on (49) for a final time value tend = tskend and (50) for the energy consumption, multiple closed-loop batches are simulated with different pe values. All other control algorithm parameters are kept constant. Notice that since the stopping condition (49) is a discrete-time threshold, the resulting tend values are discontinuous. As presented in Fig. 5, by increasing pe, the batch duration tend increases and the total energy consumption Wtb decreases. Both

Fig. 5. Simulations batch time and energy consumption for different weights on energy index pe (all other parameters being equal).

97

Control Engineering Practice 64 (2017) 88–101

F. Gagnon et al.

Table 3 LOD moisture content in wet basis (final target: 1%). Open-loop operation

Closed-loop operation

Sample

Init (%)

Final (%)

Init (%)

Final (%)

A B C

13.60 13.33 13.28

0.50 0.43 0.60

13.25 13.74 13.15

1.15 1.77 1.75

mean

13.40

0.51

13.38

1.56

6.1.3. Closed-loop: low weight on energy consumption The weight on the energy index is:

pe = 0.9 × 10−4 As shown in Fig. 7, in this case, the airflow setpoint is slowly decreased at the sixth minute. The resulting batch time is 20 min and the total energy consumption is 2.55 MJ . Fig. 7. FBD simulator in closed-loop: low weight on energy consumption.

6.1.4. Closed-loop: high weight on energy consumption The weight on the energy index is:

benchmarks change drastically for smaller values of pe and tend to plateau for larger values of pe > 18 × 10−4 . This saturation is due to the manipulated input minima u min . To keep energy consumption as low as possible with high pe values, the control algorithm forces both manipulated inputs on their respective lower bounds for the entire batch.

pe = 9.0 × 10−4 For high pe values, Fig. 8 shows that the airflow setpoint is promptly decreased at the beginning of the batch. With even larger pe values, both manipulated inputs u would be forced to their respective minimum u min for the entire batch. The resulting batch time is 31 min and the total energy consumption is 1.52 MJ .

6.1.2. Closed-loop: priority on batch time To prioritize cycle time, energy consumption is ignored in the NMPC criterion (45):

6.2. Pilot unit results

pe = 0

Because of a lack of equipment availability, the same raw material is granulated with a different granulator than the one used for the PLS and FBD model calibrations. This represents a real-world disturbance that should be absorbed by the integral action and illustrates the robustness of the proposed control law. The NIR probe provides the particle moisture content measurement. To verify the accuracy of the measurement, some powder samples are taken before, during and after the batch and measured

Fig. 6 shows that the airflow setpoint is maintained at its maximal value for the full batch while the inlet temperature setpoint is slowly decreased at the eleventh minute to have the particle moisture content smoothly reaching its desired value at the end of the batch and limit their temperature. Since the airflow has a very little impact on drying in the falling-rate period (Villegas et al., 2010), the airflow setpoint is kept high at the end. The batch time is 15 min and the total energy consumption is 3.40 MJ .

Fig. 9. Pilot FBD in closed-loop: weight on energy consumption (integral action xli is Fig. 8. FBD simulator in closed-loop: high weight on energy consumption.

shown with error bars).

98

Control Engineering Practice 64 (2017) 88–101

F. Gagnon et al.

7. Conclusion

with the LOD analyzer. To improve the LOD measurement of the moisture content before and after the batch, three samples at different locations in the bed were taken and their measurements averaged.

Advanced process control of particle moisture content in a fluidized bed drying process has proven to be more advantageous than the usual open-loop operation. By explicitly including energy consumption in a nonlinear predictive controller cost function, it is possible to either aim at a reduced cycle time or to take into account the energy consumption. An in-line NIR probe was successfully pre-calibrated on particle moisture content, but with identified improvement opportunities. Relying on the NIR measurements, a specific moisture content setpoint is reached while constraining particle temperatures to limit issues like under drying, over drying and particle overheating. The proposed strategy was successfully tested on a pilot FBD and exhibited robustness to modeling mismatches (for instance, different powder granulations). When solving the NMPC and MHE problems with an optimization routine, it is capital to include exception handling and an alternate solution (from the preceding sampling period for instance) should be planned and always ready to be employed. This problem did not occur during the tests presented in this paper. The nonlinear optimizations require heavy calculations. It would certainly be recommended to perform a complete compilation of the control law instead of programming with an interpreted language like MATLAB. Also, a commercial programmable logic controller would be more reliable for data acquisition, data logging and human-machine interface than a USB data acquisition card with LabVIEW. Further investigations include the assessment of a control solution based on a simpler model (linear model, family of linear models, etc.). To reduce the process instrumentation costs, the MHE should be evaluated as a potential particle moisture content virtual sensor. Finally, since the FBD model was first developed for a continuous process (see Lai et al., 1986), it would be interesting to see if a similar control approach could be adapted for this process.

6.2.1. Usual open-loop operation A basic open-loop control strategy is applied to the pilot FBD. As it is often the case in industry, the drying time, the inlet temperature setpoint and the inlet air flow rate setpoint are selected before the batch. Their respective values are 30 min, 65 °C and 150 m3/h. As shown in Table 3, the particles are slightly over dried by 0.49% at the end and their temperature is about 53.3 °C, which may be inappropriate. The calculated total energy consumption with tend = 30 min is 3.00 MJ . 6.2.2. Closed-loop: weight on energy consumption The NMPC control approach with a low weight on energy consumption, as described in Section 6.1.3, is applied to the pilot FBD. The impact of pe weight value over simulations, as depicted in Section 6.1.1, were following the same tendencies on the pilot FBD. Fig. 9 shows that the measured y and estimated ylk moisture contents (in blue) drift apart. The situation would be even worse without the integral action (the estimate xli ). In that case, the estimated moisture content ylk would pass by the top of the error bars. The correction of the integral action is guaranteed to be perfect only in steady-state. Nonetheless, the calculated manipulated inputs bring the particle moisture content near the desired 1% within 13 min. The calculated total energy consumption based on the real measurements is 2.15 MJ . The particle temperature at the end of the batch is very acceptable without the corresponding constraint of 45 °C being hit during the drying. An important variation detected from the use of the NIR probe for several batches is that the NIR moisture content measurement error of prediction varies depending on the air flow rate, since the PLS model was calibrated for non-static powders at higher flow rates. The same phenomenon was also noticed by Obregón et al. (2013). With this closed-loop test, Table 3 shows that the batch is in fact (according to LOD measurements) slightly under dried by 0.56%. When trying to reduce the energy consumption, the air flow rate decreases at the end of the batch making this problem more apparent. A better PLS model by varying the flow rate at calibration, or a family of PLS models depending on the flow rate could be a potential improvement. Also, the particle moisture content distribution at the end of the batch seems a bit more spread than in the open-loop case with a larger inlet flow rate (Table 3). This may be explained by a decrease in the mass diffusivity during the drying at lower flow rates. It would probably be worth defining a better stopping condition.

Acknowledgements The authors would like to thank Fonds de recherche du Québec Nature et technologie (FRQNT) and National Sciences and Engineering Council of Canada (NSERC) for the funding of this research, and to Pfizer Montreal and Viavi Solutions for providing the necessary raw materials and laboratory equipment respectively. Special thanks to Mr. Steve Hammond (Pfizer) for his time and support, and to Prof. Haigang Wang (The Institute of Engineering Thermophysics) for sharing some of his works on FBD simulations.

Appendix A. FBD model: additional equations All the other FBD model parameter calculations are listed in this appendix with their definitions and units given in the Nomenclature section. They are gathered into subsections, which follow the correct calculation sequence to respect dependency. A.1. Fluidization regime parameters From Wang et al. (2007):

υmf =

ηdp2

3 2 εmf ϕ

150μg 1 − εmf

with η = ( ρws − ρg )g

(A.1)

⎛ μ2 ⎞0.029⎛ ρ ⎞0.021 g ⎜⎜ g ⎟⎟ εmf = 0.586ϕ−0.72⎜⎜ 3 ⎟⎟ ηd ⎝ ρws ⎠ ⎝ p⎠

(A.2)

From Li & Duncan (2008a):

99

Control Engineering Practice 64 (2017) 88–101

F. Gagnon et al.

ρs

ρds =

1+

ρs χ ρw pc

and ρws = (1 + χp )ρds (A.3)

A.2. Bubble parameters From Kunii and Levenspiel (1991):

⎛ −0.3z ⎞ db(z ) = dbm − (dbm − db0 )exp ⎜ ⎟ ⎝ Dc ⎠

(A.4)

⎞ ⎛ πD 2 dbm = 26 × 10−6⎜⎜ c (υ0 − υmf )⎟⎟ ⎠ ⎝ 4

(A.5)

0.4

For perforated plate distributor: 0.4 1.3 ⎛ υ0 − υmf ⎞ ⎜ ⎟ g 0.2 ⎝ Nor ⎠

db0 =

(A.6)

where:

⎧ 2 L −2 for triangular orifice pattern ⎪ or Nor = ⎨ 3 ⎪ L −2 for square orifice pattern ⎩ or

(A.7)

For porous plate distributor:

db0 = 0.376(υ0 − υmf )2

(A.8)

From Li & Duncan (2008a):

υb = υ0 − υmf + υbr δb =

υbr = 0.711 gd͠ b

with

(A.9)

υ0 − υmf υb

(A.10)

From Lai et al. (1986):

Hf =

Hmf 1 − δb

(A.11)

A.3. Heat and mass transfer coefficients From Li & Duncan (2008a) or Kunii and Levenspiel (1991):

Kbe =

1 1/ Kbc + 1/ Kce

Kbc = 4.5

(A.12)

υmf Dg1/2g1/4 + 5.85 5/4 d͠ b d͠ b

(A.13)

⎛ε Dυ ⎞ mf g b ⎟ Kce = 6.78⎜⎜ 3 ⎟ ⎝ d͠ b ⎠

(A.14)

1 1/ Hbc + 1/ Hce

(A.15)

1/2

Hbe =

⎞1/2 ⎛ ⎜kgρg cg⎟ υmf ρg cg ⎠ ⎝ Hbc = 4.5 + 5.85 g 5/4 5/4 d͠ b d͠

(A.16)

⎞1/2 ⎛ εmf υb ⎞ ⎛ Hce = 6.78⎜ρg cgkg⎟ ⎜⎜ 3 ⎟⎟ ⎠ ⎝ d͠ ⎠ ⎝

(A.17)

b

1/2

b

From Wang et al. (2007):

⎛ Tp ⎞3/2 ⎟⎟ Dg = Dg0⎜⎜ ⎝ Tref ⎠

(A.18)

100

Control Engineering Practice 64 (2017) 88–101

F. Gagnon et al.

γ = 103(3168 − 2.4365Tref )

(A.19)

From Li & Duncan (2008a) and Lai et al. (1986):

σ=

hsρg Dg kg

Nup = 2 +

Rep =

and hs =

Nupkg dp

(A.20)

1.8Rep1/2Pr1/3

dpU0ρg

and

μg (1 − εmf )

(A.21)

Pr =

c pμg kg

(A.22)

From Li and Finlayson (1977), for cylindrical bed packing, we have: 0.93 ⎛ kg ⎞⎛ dpU0ρg ⎞ ⎟ h w = 0.16⎜⎜ ⎟⎟⎜⎜ ⎟ ⎝ dp ⎠⎝ μg ⎠

aw =

πDcHf (1/4)πDc2Hf

=

(A.23)

4 Dc

(A.24)

A.4. Moisture content on the particle surface From Li & Duncan (2008a):

χp* = ψ1(Tp )ψ2(χp )

ψ1(Tp ) = 0.622

(A.25)

Pw(Tp ) 760 − Pw(Tp )

(A.26)

T −273.15

Pw(T ) = 100.622+7.5 T −35.15

(A.27)

⎧1 if χp ≥ χpc ⎪ ⎪ ν⎛ ν ⎞ ψ2(χp ) = ⎨ χp ⎜⎝χpc + κ ⎟⎠ if χp < χpc ⎪ ⎛ ⎞ ν ⎜χ ν + κ ⎟ ⎪ χpc ⎩ ⎝p ⎠

(A.28)

Maciejowski, J. (2000). Predictive control: With constraints (1st ed.). Prentice Hall. McKenzie, P., Kiang, S., Tom, J., Rubin, E., & Futran, M. (2006). Can pharmaceutical process development become high tech? AIChE Journal, 52(12), 3990–3994. Mujumdar, A. S. (2014). Handbook of industrial drying (4th ed.). CRC Press. Nagy, Z. K., & Braatz, R. D. (2003). Robust nonlinear model predictive control of batch processes. AIChE Journal, 49, 1776–1786. Obregón, L., Quiñones, L., & Velázquez, C. (2013). Model predictive control of a fluidized bed dryer with an inline nir as moisture content sensor. Control Engineering Practice, 21, 509–517. Rao, C. V., & Rawling, J. B. (2002). Constrained process monitoring: Moving-horizon approach. AIChE Journal, 48, 97–109. U.S. Food and Drug Administration (2004). Pharmaceutical CGMPS for the 21st century – a risk-based approach: final report. Technical report. Department of Health and Human Services, September. Villegas, J., Duncan, S., Wang, H., Yang, W., & Raghavan, R. (2010). Optimal operating conditions for a batch fluidised bed dryer. IET Control Theory Appl, 4, 294–302. Wang, H., Dyakowski, T., Senior, P., Raghavan, R., & Wang, W. (2007). Modelling of a batch fluidised bed drying of pharmaceutical granules. Chemical Engineering Science, 62, 1524–1535. Wold, H. (1975). Quantitative Sociology: International perspectives on mathematical and statistical model building. , in: Blalock, H. M., Aganbegian, A., Borodkin, F., Boudon, R., & Cappecchi, V. (Eds.). (1975). Path models with latent variables: The NIPALS approach . . Academic Press, 307–357.

References Bonvin, D., Srinivasan, B., & Hunkeler, D. (2006). Control and optimization of batch processes. IEEE Control Systems Magazine, 26(6), 34–45. Briens, L., & Bojorra, M. (2010). Monitoring fluidized bed drying of pharmaceutical granules. AAPS PharmSciTech, 11(4), 1612–1618. Demers, A.-M., Gosselin, R., Simard, J.-S., & Abatzoglou, N. (2012). In-line near infrared spectroscopy monitoring of pharmaceutical powder moisture in a fluidised bed dryer: an efficient methodology for chemometric model development. The Canadian Journal of Chemical Engineering, 90, 299–303. Hao, J., Zhao, Y., Ye, M., & Lius, Z. (2016). Influence of temperature on fluidized-bed catalyst attrition behavior. Chemical Engineering Technology, 39(5), 927–934. Kunii, D., & Levenspiel, O. (1991). Fluidization engineering (2nd ed.). ButterworthHeinemann. Lai, F. S., Chen, Y., & Fan, L.-T. (1986). Modelling and simulation of a continuous fluidized-bed dryer. Chemical Engineering Science, 41(9), 2419–2430. Li, C.-H., & Finlayson, B. A. (1977). Heat transfer in packed beds – a reevaluation. Chemical Engineering Science, 32, 1055–1066. Li, M., & Duncan, S. (2008). Dynamic model analysis of batch fluidized bed dryers. Particle and Particle Systems Characterization, 25, 328–344. Li, M., & Duncan, S. (2008). Model-based nonlinear control of batch fluidized bed dryers. Particle and Particle Systems Characterization, 25, 345–359. Ljung, L. (1999). System identification – theory for the user (2nd ed.). Prentice Hall.

101