Next Article in Journal
Influence of Fluid Properties on Intensity of Hydrodynamic Cavitation and Deactivation of Salmonella typhimurium
Next Article in Special Issue
Effect of Freeze-Drying on Quality and Grinding Process of Food Produce: A Review
Previous Article in Journal
Outdoor Large-Scale Cultivation of the Acidophilic Microalga Coccomyxa onubensis in a Vertical Close Photobioreactor for Lutein Production
Previous Article in Special Issue
Drying Applications during Value-Added Sustainable Processing for Selected Mass-Produced Food Coproducts
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Model-Based Real Time Operation of the Freeze-Drying Process

by
Carlos Vilas
1,*,
Antonio A. Alonso
1,
Eva Balsa-Canto
1,
Estefanía López-Quiroga
2 and
Ioan Cristian Trelea
3
1
(Bio)Process Engineering Group, Instituto de Investigaciones Marinas (CSIC), Eduardo Cabello, 6. 36208 Vigo, Spain
2
School of Chemical Engineering, University of Birmingham, Edgbaston, Birmingham B152TT, UK
3
INRAE, AgroParisTech, UMR SayFood, Université Paris-Saclay, F-78850 Thiverval-Grignon, France
*
Author to whom correspondence should be addressed.
Processes 2020, 8(3), 325; https://doi.org/10.3390/pr8030325
Submission received: 3 January 2020 / Revised: 12 February 2020 / Accepted: 4 March 2020 / Published: 10 March 2020
(This article belongs to the Special Issue Drying Kinetics and Quality Control in Food Processing)

Abstract

:
Background: Freeze-drying or lyophilization is a dehydration process employed in high added-value food and biochemical goods. It helps to maintain product organoleptic and nutritional properties. The proper handling of the product temperature during the operation is critical to preserve quality and to reduce the process duration. Methods: Mathematical models are useful tools that can be used to design optimal policies that minimize production costs while keeping product quality. In this work, we derive an operational mathematical model to describe product quality and stability during the freeze-drying process. Model identification techniques are used to provide the model with predictive capabilities. Then, the model is used to design optimal control policies that minimize process time. Results and conclusion: Experimental measurements suggest splitting the process into two subsystems, product and chamber, to facilitate the calibration task. Both models are successfully validated using experimental data. Optimally designed control profiles are able to reduce the process duration by around 30% as compared with standard policies. The optimization task is introduced into a real time scheme to take into account unexpected process disturbances and model/plant mismatch. The implementation of the real time optimization scheme shows that this approach is able to compensate for such disturbances.

1. Introduction

Freeze-drying is a preservation technique widely used in pharmaceutical and food processing applications [1,2]. In freeze-drying, products are first frozen and then dehydrated by sublimation of the ice formed during freezing (primary drying) and desorption of the unfrozen water (secondary drying) [3,4]. Such particular processing conditions are key to delivering high-quality, shelf-stable dried products. For example, low processing temperatures avoid the degradation of heat-sensitive compounds (e.g., pharmaceutical active principles or nutrients and flavors in foods), while freezing and subsequent sublimation create highly porous microstructures that enhance rehydration/dissolution of food [5,6,7] and drug powders [8]. Characteristic low moisture contents make freeze-drying both time-consuming—cycles are typically over 24 h—and energy intensive [9,10]. The challenge here lies in shortening processing times without compromising the quality of the final product—for example, higher temperatures will increase drying rates, but could cause matrix collapse and/or biodegradation.
In this context, optimization protocols based on mechanistic, first-principles models [11,12,13] can help to achieve the desired time reductions while keeping quality, and also provide a robust tool to deal with process and/or product disturbances during freeze drying. In recent years, the pharmaceutical sector, operating under strict regulations, has adopted a Quality by Design (QbD) approach that looks for product quality in-line, using model-based optimization strategies combined with process analytical tools (PAT)—temperature and pressure sensors providing dynamic information to estimate critical product and process parameters (e.g., heat and mass transfer resistances in both dried product and vials) [14,15,16,17,18,19].
However, such approaches are not yet exploited in freeze-drying of foodstuff, in part due to the lack of operational models suitable for real time applications—i.e., current literature in food systems [20,21,22,23] does not include simple yet physically meaningful dynamic freeze-drying models—and in part because foods are complex biomaterials, with properties strongly dependent on temperature and pressure (e.g., processing conditions) that are difficult to characterize [13]. Previous works on food systems are scarce, and use either lumped inventory models [3,24] that miss temperature/pressure dependences in the system or very detailed, multiscale models that result computationally involved and include a large number of unknown parameters [25].
To fill this gap, we build on previous work and use a physics-based, operational model for freeze-drying [26] (i.e., reduced attending to separation of time scales) to develop a robust two-layer scheme for real time monitoring, control and optimization of freeze-drying in food matrices. This scheme uses input-output data recorded with standard sensing devices (thermocouples and Pirani gauge) to identify model parameters, such thermo physical properties of the food matrix or characteristic transfer coefficients related to the equipment, that vary with operating conditions (e.g., shelf temperature, chamber pressure). The outputs include Pirani [27] and total chamber pressure, as well as product temperature at the (approximated) bottom of the sample. The inputs are both shelf and condenser temperatures, which act as controls for the system. Once both product and equipment are characterized, the scheme can be used to determine sublimation front position and process end point, to update operating conditions under process disturbances or to design optimal operating policies—off-line and in real time (RTO)—that minimize processing times while ensuring final food quality (i.e., reaching target moisture content while avoiding structure collapse). We assume that, at the beginning of the drying process, the product is frozen and its temperature is uniform in the spatial domain. Therefore, focus of this study is on the drying phase of the process.
Overall, the contribution of this work is two-fold:
  • It shows how the identification problem (i.e., characterization of product and equipment) can be simplified by solving a (decoupled) sequence of parameter identification problems, this is, by considering food matrix and freeze-drying chamber as two separate subsystems, which are described by physics-based operational models developed to satisfy identifiability.
  • It provides a robust, digital tool that supports decision-making in real time for complex food structuring processes (such freeze-drying). This has potential for significant reductions in cost and processing times, waste generation (e.g., reducing batch rejection) and energy demand (e.g., off-line and in-line strategies can be used to minimize energy use during processing) [10,28].
To illustrate the capabilities of this approach, we have considered Lactic Acid Bacteria (LAB) as food matrix exemplar. This is a culture starter used in the manufacture of foods (e.g., yogurt, cheese, fermented meats, and vegetables) and probiotic products that is commonly freeze-dried to preserve its functionality [29].

2. Lyophilization Plant and Mathematical Model

2.1. Lyophilization Pilot Plant and Experimental Setup Description

The experiments required to perform the different tasks in this work were carried out on a LyoBeta special freeze-dryer (Telstar, Terrassa, Spain) equipped with three thermocouples, a capacitive manometer and a Pirani gauge.
A schematic representation of the pilot plant is shown in Figure 1. The product is placed over different shelves. Such shelves are connected to a heating system which is used to heat and cool down the product. As the product is heated, ice begins to sublimate starting by the top of the product. At this moment two phases can be distinguished, a dried one at the top of the product and a frozen one at the bottom of the product. As the process evolves the frozen region reduces its size till it disappears. Vapor generated by ice sublimation is removed by the action of a condenser to control the pressure in the chamber.
Lactic acid bacteria were produced by fermentation in controlled conditions of pH and temperature [29]. After concentration, the bacterial cells were re-suspended in a 1:2 cells/protective medium ratio. The protective medium was composed of 200 g L 1 of sucrose and 0.15 M of NaCl. A stainless steel tray was filled with 450 g of bacterial suspension.

2.2. A Physical Model for the Freeze-Drying Process

The freeze-drying process mainly consists of the following three stages [30]:
  • Freezing. The product is frozen in controlled conditions to avoid possible damage, by crystal growth, to the food or biological material. Our model does not include this stage. The product temperature is assumed to be uniform at the end of the freezing stage, and it is used as the initial temperature for the primary drying.
  • Primary drying. In this stage, ice is removed from the product by sublimation. Pressure conditions are kept below the triple point and the product is heated from the bottom. An excessive temperature increase in this stage will cause product to collapse so it must be kept below a given value, see Section 2.2.4 for details.
  • Secondary drying. The aim is to remove water bound to the solid matrix by desorption. This stage allows for reaching low moisture contents. The food product is more stable during secondary drying so its temperature may be increased to accelerate the process.
Accurate models describing the primary drying stage are computationally involved as this stage involves complex heat and mass transfer mechanisms, as well as ice sublimation. Besides, such models usually include a large number of unknown parameters which often poses identifiability problems when trying to estimate parameter values from experimental data. Operational models with good predictive capabilities can be obtained by considering the following, widely employed, assumptions [26,31,32]:
  • The frozen region has uniform heat and mass transfer properties.
  • The interface between the frozen and dried layers (sublimation front) is continuous and has infinitesimal thickness.
  • Vapor and ice are at equilibrium at the interface.
  • The matrix pore structure is permeable to the vapor flux and it is not deformable.
Model complexity can be further reduced by focusing on the phenomena of interest from the process control point of view and neglecting mechanisms occurring at time scales very different from such phenomena. The authors in [26] performed a careful dimensionless analysis revealing several time-scales involved during primary drying. Such analysis is used in this work to obtain the final model equations.
The following section presents an operational version of the freeze-drying model. The parameters involved in the model are summarized in Table 1.

2.2.1. The Primary Drying

The phenomenon of interest, from the process control point of view, is the heat transport in the product. Based on the time scale analysis, the authors in [26] concluded that vapor flux contribution to the energy balance can be neglected as it is much faster than heat conduction in the solid matrix. Besides, product dimensions L z and L y are two orders of magnitude larger than L x , see Table 1 and Figure 1. Therefore, the 3D problem can be accurately approximated by a 1D in the L x dimension. Under these conditions, the energy balance leads to the following set of partial differential equations (see [26] for details):
ρ d r c p , d r T d r t = κ d r 2 T d r ξ 2
ρ f r c p , f r T f r t = κ f r 2 T f r ξ 2
T i ( t , ξ ) , with i = d r , f r , represents the food product temperature; ξ are the spatial coordinates (product height); and d r , f r denote dried and frozen regions, respectively.
Note that, during primary drying, sublimation makes the front to move from the top to the bottom of the product. Front velocity w = d x d t , where x ( t ) denotes the front position, is computed using the Stefan condition [32]:
( ρ f r ρ d r ) Δ H s w = κ f r T f r ξ ξ = x + κ d r T d r ξ ξ = x
where Δ H s is the sublimation heat.
At the product top ( ξ = 0 ), radiation is considered as the main heat transport phenomena, hence the heat flux at this boundary reads as:
κ d r T d r ξ ξ = 0 = σ e p f p ( T c h 4 T d r ξ = 0 4 )
where σ is the well-known Stefan–Boltzmann constant, e p is the emissivity, f p is a geometric correction factor and T c h is the chamber temperature. At the product bottom ( ξ = L ), we consider both radiation and conduction via gas-surface collisions [33]:
κ f r T f r ξ ξ = L = h L ( T s h T f r ξ = L )
where T s h is the shelf temperature. Heat transfer coefficient h L is computed as a function of the chamber pressure, P c h , (see [33] for details):
h L = h L , 1 + h L , 2 P c h 1 + P c h h L , 3
At the front ( ξ = x ), continuity of temperature across the front is considered, i.e.,
T d r ξ = x = T f r ξ = x = T f r o n t
At this boundary, vapor and ice are assumed to be in equilibrium so that the front temperature ( T f r o n t ) can be computed from the front pressure using the Clausius–Clapeyron equation:
T f r o n t = 1 1 273.11 K c l a p log P f r o n t 611.72
We use Darcy’s law, which states that the vapor flux per unit area is proportional to the pressure gradient, to compute the front pressure ( P f r o n t ):
P ξ = μ K w ( ρ f r ρ d r ) ρ v
where K corresponds with the permeability of the porous matrix; μ represents the vapor viscosity and ρ v the vapor density. It has been experimentally observed that the mass transfer coefficient in Equation (9) has a linear dependency with the thickness of the dried layer (x). Furthermore, in many freeze-dried products there exists a mass resistance for zero layer thickness [33]. Such resistance, which depends on the chamber pressure, can be attributed to a crust formed on the product or the pore creation process at the sublimation front and should be included in the mass transfer coefficient expression. Taking these issues into account, Equation (9) results into:
P ξ = K c w ( ρ f r ρ d r )
where
K c = k 1 P c h + k 2 x
with k 1 and k 2 being given parameters.
Advanced numerical techniques, like the Arbitrary Lagrangian-Eulerian (ALE) algorithm [32], are usually employed to solve the moving front problem. In general, such techniques require the use of specialized software. We will apply the Landau transform [34,35] to overcome such limitation. As a result, we will obtain an equivalent system with fixed domain in which classical numerical techniques, such as the finite element method, can be applied. See Appendix A for details.

2.2.2. The Secondary Drying

As mentioned before, at this stage of the process there is no frozen water and the bounded water is eliminated by desorption. Since there is only one region (dried region), there is no moving boundary and classical numerical methods can be employed to solve the model equations. Equation (1) with the boundary conditions defined by Equations (4) and (5) , where ξ [ 0 , L ] , is used to describe product temperature evolution and distribution during this stage.

2.2.3. The Condenser Model

Part of the vapor produced as a result of ice sublimation is removed from the plant by the action of a condenser. The evolution of the chamber pressure ( P c h v ) is accurately described by the general gas equation, since working pressures in the chamber are low:
d P c h v d t = R T c h M w V c h d m c h v d t
R is the ideal gas constant, M w denotes the molecular weight of water and V c h is the chamber volume. Vapor accumulated in the chamber can be obtained through a mass balance between the chamber and the condenser:
d m c h v d t = ϕ p v ϕ c v
where ϕ p v and ϕ c v denote, respectively, the flux of vapor from the product to the chamber and from the chamber to the condenser. Flux ϕ p v corresponds with the flux of ice being sublimated, i.e.,:
ϕ p v = ( ρ f r ρ d r ) L z L y d x d t
with L z and L y being, respectively, the length and width of the product sample.
The second term on the RHS of Equation (12) ( ϕ c v ) can be described using the mass transfer theory in binary gas mixtures, including bulk flow and mutual diffusion terms [36]:
ϕ c v = 1 β T c h ln P c h P c v P c h P c h v
where β = l e f f R L z L y M w P c h D v N . D v N and l e f f are the mutual diffusion coefficient and the effective duct length, respectively. The details regarding the derivation of Equation (14) are presented in [36].
Note that Equation (13) can only be used during primary drying where a sublimation front, that evolves with velocity w = d x d t , exists. During secondary drying, flux ϕ p v is assumed to be zero. In the real process, the heat coming trough the edges of the shelf accelerate the sublimation in those parts of the product that are closer to the borders. As a consequence, vapor flux and vapor pressure decrease gradually in the transition between primary and secondary drying. To describe this phenomenon, 2D or 3D models are required [37]. Such models are unsuitable for RTO purposes as they require long simulation times. Alternatively, 1D models, in which the transition between primary and secondary drying is instantaneous, are operational and yet accurate to describe systems where one of the spatial dimensions is much smaller than the other two [26,31,38,39], as it is the case in the present study. Besides, chamber pressure in the transition between primary and secondary drying has little effect on the desorption kinetics [40,41] and, therefore, on the product quality.

2.2.4. The Desorption Model

During primary drying, ice sublimation and desorption of bound water are the physical mechanisms for water loss. In the secondary drying stage, no ice remains in the product so desorption is the only mechanism for water loss.
The desorption model will be used to determine final water content in the product, which is related to the final product quality, as well as the glass transition temperature, which is related to product stability during the process. This variable is used to define constraints on the foodstuff temperature that prevents product collapse.
The initially present water is distributed among two “compartments” (compartments A and B), corresponding to different physical states and/or interactions with the solid matrix (bacterial cell, membrane, cryoprotectant, etc.), which have their own temperature-dependent drying kinetics [40]:
d c w , i d t = 1 τ i c w , i c w , i e q ; with i = A , B
with c w , i e q being the equilibrium water content for compartment i at equilibrium. τ i , desorption time functions, are assumed to obey an Arrhenius-like relationship [41]:
τ i = τ i r e f exp E a R 1 T p 1 T r e f
where τ i r e f is the value of the time constant at an arbitrarily fixed reference temperature and E a is the activation energy. Equilibrium water content for compartment i is related to the total equilibrium water content as:
c w , i e q = α i c w e q ; with α A + α B = 1
c w e q is computed using the water activity ( a w ) based on the classical Guggenheim-Anderson-Boer (GAB) equation for the sorption isotherm [29,42]:
c w e q = M g C g K g a w 1 K g a w 1 + K g a w C g 1
where M g , K g and C g are given, product dependent, constants. Water activity is the ratio between the chamber vapor pressure and the equilibrium vapor pressure, i.e., a w = P c h v / P e q v where P e q v is related to product temperature as [43,44]:
P e q v = 3.5941 × 10 12 exp 6144.96 T p
which is valid for ice and:
P e q v = 8.4755 × 10 11 exp 5750.38 T p 1 0.0045 ( T p 273.15 )
whenever water is involved.
The glass transition temperature is related to the water content as follows [29]:
T g ( c w ) = K T g T g , 1 c w + ( 1 c w ) T g , s K T g c w + ( 1 c w )
where K T g , T g , 1 and T g , s are given parameters.
As mentioned above, T g can be used to predict product integrity. In this regard, product temperature should not exceed T g since foodstuff might collapse.

3. Strategies for Model Parameter Identification

The aim of parameter estimation is to find the values of the unknown model parameters that minimize the distance between model simulations and experimental data. To that purpose, the following elements are required:
  • A set of model parameters to be estimated. Usually, these are parameters that are difficult to measure or cannot be found in the literature.
  • Experimental measurements to compare against model results. Measured variables are usually state variables (e.g., temperature, concentration, pressure) or combinations of state variables. The quantities to be measured are known as the observables. Changes on the unknown model parameters should have an effect on the observables, otherwise such parameters cannot be estimated.
  • A measure of the distance between model predictions and experimental data (cost function).
  • An optimization algorithm to find the parameter values that minimize the cost function.
Let us denote by θ the set of unknown model parameters and by y the observable. In this work, the cost function is defined as the Root Mean Square Error (RMSE) value. The parameter estimation problem is formulated as, find:
min θ R M S E ; with R M S E = k = 1 n e i = 1 n s , k y m e a s ( t i , k ) y ¯ m o d ( t i , k ) 2 k = 1 n e n s , k
subject to the model dynamics presented in previous section. In Equation (20), n e is the number of experiments for estimation purposes, n s , k is the number of sampling times in experiment k and y ¯ m o d ( t i , k ) represents the model prediction of the observable at time t i in experiment k.
The optimization method chosen to solve the problem is the sequential hybrid global-local Enhanced Scatter Search (ESS) [45] which is included in the AMIGO2 toolbox [46].
The freeze-drying system was divided in two subsystems to facilitate the estimation task, the condenser and the product, each of them with an associated mathematical model as shown in the previous section. The procedure is schematically represented in Figure 2.
The control variables, represented in Figure 2 as arrows entering the subsystems, correspond to those process variables that can be directly manipulated to actuate into the process. Three controls are available in this pilot plant: condenser pressure ( P c ), chamber pressure ( P c h ) and shelf temperature ( T s h ). These variables are also measured so the real value can be used for parameter estimation tasks. Arrows leaving the subsystems in Figure 2 represent the state variables. The chamber vapor pressure ( P c h v ) is the only state variable computed in the condenser model. State variables in the product subsystem are the front position (x) and velocity (w) as well as product temperature ( T p ), the latter distributed along the product height. Both subsystems, condenser and product, are coupled through two state variables, P c h v and w.
In this system, Pirani pressure ( P p ) and total chamber pressure ( P c h ) are also measured. Since chamber vapor pressure is related to P p and P c h through the formula [36]:
P c h v = P p P c h α 1 ; α = 1.6 ,
we can consider that P c h v is indirectly measured. These indirect measurements can be used, in the product model, so that the condenser model can be neglected in the estimation of the product model parameters. In this regard, the estimation task will proceed in two steps:
  • Estimation of the unknown parameters involved in the product model, using the indirect measurements of P c h v instead of the condenser model. On the one hand, this reduces the number of parameters to be estimated together since the unknown condenser parameter is not taken into account. On the other hand, possible modeling errors on the condenser model are avoided.
  • Estimation of the unknown parameters involved in the condenser model. As it is sketched in Figure 2, this estimation requires the product model and the parameters estimated in the previous step to compute the front velocity. Note that, if a procedure to measure w, such as an on-line sensor of front velocity, were available, product model could be disregarded in this step.
The following freeze-drying experimental protocol was applied is the different tests:
  • Product freezing using −50 ° C as a set point (cooling rate of 0,6 ° C   min 1 ).
  • Primary drying varying from freezing temperature to −10, 0 or 20 ° C depending on the experiment.
  • Secondary drying using 25 ° C as the set point.
  • Total chamber pressure was kept at 20 or 60 Pa depending on the experiment.

3.1. Parameter Estimation in the Product/Chamber Subsystem

The subset of model parameters to be estimated is θ = [ κ d r , k 2 , h L , 2 ] , i.e., the heat conductivity of the dried region; product/chamber mass resistance coefficient—see Equation (11); and heat transfer coefficient—see Equation (6).
The observable in this subsystem is the temperature at the bottom of the product. Because of limitations on the equipment, the precise location of such thermocouple is not known with precision. Therefore, a mean value of the temperatures in the bottom region of the product is used as the observable y = T p , m e a s . Such region is 1 m m long.
In this subsystem, Equation (20) may be rewritten as:
min κ d r , k 2 , h L , 2 R M S E ; with R M S E = k = 1 n e i = 1 n s , k T p , m e a s ( t i , k ) T ¯ p , m o d ( t i , k ) 2 k = 1 n e n s , k
where T ¯ p , m o d ( t i , k ) represents the mean value of the temperatures in the bottom region of the product.
The experimental data employed in the parameter estimation were obtained from five experiments. Each experiment differs from the others in the control variables, i.e., shelf temperature ( T s h ) and chamber pressure ( P c h ). Besides, although the set point for the freezing temperature was set to −50 ° C , the refrigerant group was not able to reach such temperatures. Experimental measurements showed that the actual values turn out to range between −37 and −43 ° C , depending on the experiment. Therefore, initial conditions for primary drying also differ from one experiment to another.
A cross validation procedure was applied [47] to robustly estimate the parameters. In this procedure, part of the experimental data is saved for validation purposes, i.e., such data is not used in the estimation task. In this case, we use four of the five experiments for the parameter estimation task, leaving the remaining one for validation. To increase robustness of the results, we have performed five estimations, each one excluding a different experiment.
The product model, together with the measured values of the temperature at the bottom of the product ( T p , m e a s ) in the experiments considered for estimation, were implemented in the AMIGO2 toolbox. RMSE values were obtained by comparing experimentally measured values of the product temperature with the simulation results.
The estimation results, which include optimal parameter values as well as RMSE values for both estimation and validation tests, are summarized in Table 2.
The values found for the parameters are in accordance with the expected physical range. Besides, variability in the parameter estimates of the cross validation procedure for h L , 2 , k 2 and κ d r is 4.4%, 14.7% and 23.7%, respectively. While the variability of κ d r is much larger than the variability of the other two parameters, the results are considered to be robust.
Figure 3 illustrates the predictive capabilities of the model.
Four subplots are shown, each one corresponding to the validation experiment in the cross validation procedure. Blue line represents the shelf temperature (control variable) evolution in each experiment whereas dashed lines correspond to model predictions of the product temperature evolution at different spatial points. The figure shows that model predictions at the bottom of the product (continuous black line) reproduce the experimental behavior (black dots). In particular, the end of the primary drying stage as well as the product temperature during such stage are predicted by the model with a reasonable degree of accuracy. Temperature during secondary drying is also accurately described by the model. However, some mismatch between predicted and measured temperatures can be observed in the transition from primary to secondary drying (highlighted in Figure 3 by two vertical black lines). In fact, the model seems to predict a faster temperature increase than the one measured by the thermocouple.
So far, no clear explanation has been found to such disagreement although it most probably can be attributed to an interphase structure which is more complex than the sharp transition between a frozen and a dried region as it is considered by the model. Nevertheless, as mentioned on Section 2.2.3, this does not have a significant impact on the product quality. Besides, product temperature predicted by the model in the transition is usually larger than the experimental measurements. Therefore, as it will be described in Section 4.1, operation policies that satisfy safety constraints according to model predictions will also satisfy such constraints in the experiments.

3.2. Parameter Estimation for the Condenser Model

The unknown model parameter in the condenser model is coefficient θ = β , see Equation (14). In this plant, it is possible to experimentally measure the total chamber and Pirani pressures so that the chamber vapor pressure ( P c h v ) can be obtained, see Equation (21). y = P c h v is used as the observable in the estimation task.
In this case, the optimization problem defined by Equation (20) reads as:
min β R M S E ; with R M S E = k = 1 n e i = 1 n s , k P c h , m e a s v ( t i , k ) P c h , m o d v ( t i , k ) 2 k = 1 n e n s , k
The same experiments used in the previous section are employed here to estimate the condenser parameter β . In this case we use both the product model, with the parameter values computed in previous section, and the condenser model.
As in previous section, cross validation is used to evaluate the performance of the model. The estimation results, including optimal value of β and the RMSE values corresponding to both the estimation and validation experiments, are presented in Table 3.
The estimated values of β are within the expected physical range and close to the values presented in the literature [36]. Besides, the variability in the parameter estimates for the different cases in the cross validation is lower than 4 % which evidences the estimation robustness.
Model predictive capabilities, in terms of chamber pressure prediction, are illustrated in Figure 4. As in the previous case, the four subplots correspond to the validation experiment in the cross validation tests.
Only primary drying is represented in the figures since during secondary drying vapor pressure is almost zero. As mentioned in Section 2.2.4, the transition between primary and secondary drying is assumed to be instantaneous. Simulation results for the chamber vapor pressure are represented by continuous black lines whereas black circles correspond to experimental measurements. As shown in the figure, the model is able to reproduce the experimental behavior with a high degree of accuracy.

4. Real time Optimization of the Freeze-Drying Process

4.1. Off-Line Dynamic Optimization

The objective of the dynamic optimization problem is to find the shelf temperature time profile that minimizes the process duration while fulfilling product quality requirements at final time and product temperature constraints during the whole process to avoid product collapse. The other control variable on this system, i.e., chamber pressure, as well as initial conditions for the state variables are taken from the first experiment employed in the parameter estimation task. Mathematically, the optimization problem is stated as:
min T s h t f
subjet to 50 ° C T s h 25 ° C 1 60 ° C s 1 d T s h d t 1 60 ° C s 1
and to model dynamics, see Section 2.2.
If no other constraints are considered, it seems clear that the optimal profile would be that one corresponding to the maximum shelf temperature allowed, i.e., 25 ° C . However, this would result into product collapse and the rejection of the batch. In this way, constraints on the allowed product temperature are defined as function of the glass transition temperature, see Equation (19), to avoid product collapse:
T p ( ξ , t ) T g ( ξ , t ) 0
In Section 3.1, validation results showed that, in some experiments, the model might overstimate product temperature in the transition between primary and secondary drying. Mathematically, this can be formulated as: T p , m e a s ( ξ , t t ) < T p , m o d ( ξ , t t ) , where t t represents the transition period time. If the product temperature predicted by the model satisfies constraint (25), then product temperature in the experiment also satifies it:
T p , m e a s ( ξ , t ) T g ( ξ , t ) < T p , m o d ( ξ , t ) T g ( ξ , t ) 0
We handle constraints by minimizing the square of the difference, T p , m o d ( ξ , t ) T g ( ξ , t ) , in Equation (26) as suggested in [48]. This ensures that condition in Equation (26) is always satisfied.
Product quality is related to the water content at the end of the process. The lower the water content, the better the quality. The following constraint is used to ensure product quality:
ξ c w ( ξ , t f ) d ξ c w , m a x
The value c w , m a x = 9.15 × 10 3 k g water k g total was taken from the final water content in the first experiment.
There are several alternatives for the solution of dynamic optimization problems [11]. In this case, we have chosen the control vector parameterization (CVP) approach [49], because of its capacity to handle large-scale dynamic optimization problems without solving very large non-linear programming (NLP) problems and without dealing with extra junction constraints [50]. Basically, the control vector parameterization proceeds dividing the process duration into a number of elements and approximating the control functions typically using low order polynomials. The polynomial coefficients become the new decision variables and the solution of the resulting NLP problem (outer iteration) involves the system dynamics simulation (inner iteration).
NLP problems arising from the application of the control vector parameterization method are frequently multimodal (i.e., presenting multiple local optima), due to the highly nonlinear nature of the dynamics [51]. In this work, we will use a method that combine global stochastic methods with local gradient based methods which has proven to perform well for solving NLP problems. The method selected, eSS, is based on the scatter search method, see [45] for details.
In this problem, we have chosen piecewise linear functions, with seven discretization points, to approximate the control in the CVP approach. Figure 5 represents the evolution and spatial distribution of the product temperature (gray lines on the top figures) and water content (bottom figures) as well as the evolution of the shelf temperature (black lines on the top figures).
As expected, the shelf temperature profile on the optimally designed experiment (see right column subplots) starts at a low value, to avoid product damage, and increases as the product becomes more stable. Note that this profile is more aggressive than the one used in experiment 1 and, as a result, the final time is reduced by a 30% (13 h ). Despite the resulting sharp shelf temperature profile, constraints defined in Equations (25)–(27) are fulfilled, i.e., the product integrity is ensured and quality, in terms of final water content, are the same in both experiment 1 and the optimally designed one.

4.2. Real Time Optimization

The optimal profile for shelf temperature obtained in previous section was computed assuming that the model represents perfectly the plant behavior and the controller is able to follow precisely the set points. However, in real situations, plant/model mismatch is unavoidable and the PI controller implemented in the freeze-drying plant may be not able to exactly follow the piecewise linear optimal profile. Besides, unexpected plant disturbances, which cannot be taken into account in the off-line profile, may emerge during process operation.
To take into account all these issues, a real time optimization (RTO) scheme is implemented in the plant. The main idea behind the RTO scheme is to periodically use process information to recompute on-line the optimal profiles. The steps of the RTO procedure are the following:
  • An optimal profile for the shelf temperature is computed as in previous section.
  • Such profile is sent to the freeze-drying plant as the set point for the PI control.
  • As the process is running, measurements of relevant variables are being recorded. In this case, we measure, shelf temperature, chamber temperature, Pirani pressure and product temperature.
  • After a given period of time, 1 h in this case, measured information is introduced in the model. Then, a new optimal profile is computed by solving the optimization problem defined in Equation (24) and taking into account the new available information.
  • Steps 2–4 are repeated till the end of the process.
The objective is to compute the shelf temperature optimal profile that minimizes process time. Process is assumed to be finished when the product mean water content is below 0.0142 k g Water k g Total . Chamber pressure set point is fixed to 10 Pa . In this case, bounds on the shelf temperature and its slopes are:
50 ° C T s h 40 ° C ; 1 120 ° C   s 1 d T s h d t 1 120 ° C   s 1
It must be pointed out that computing optimal policies takes a few minutes. During this time, although the process continues evolving, plant measurements are not taken into account and the RTO scheme assumes that the plant follows the predicted profiles. Unexpected disturbances during this period can make the computed profile non optimal. Such disturbances are, nonetheless, taken into account when the next profile is computed.
Figure 6 illustrates an RTO implementation on the freeze-drying process. Figures at the left hand side, i.e., Figure 6a,c,e, represent shelf temperature and chamber pressure. Figure at the right hand size, i.e., Figure 6b,d,f represent product mean water content. Figure 6a,b correspond with the off-line optimal results, i.e., before the primary drying begins. Therefore, no plant measurements were used in the computation of this profile.
As process evolves, plant measurements are being recorded and new optimal profiles are being recomputed every hour using such information and following the procedure previously indicated. Subplots Figure 6c,d show the situation after 16 h, i.e., around half of the process duration. Blue and red parts of the lines correspond with the measured and recorded data. Black and gray parts of the lines represent the optimal solution from that time. As shown in the figures, the optimal profile has deeply changed. The most important factors for this are:
  • Product freezing temperature differs from the expected one, i.e., from the one used in Figure 6a,c. This was because the refrigerant group was not able to reach −50 ° C .
  • Initial chamber pressure was around 60 Pa and it was supposed to be 10 Pa . Besides, chamber pressure controller did not perform as expected and most of the time it remains below the set point (10 Pa ).
  • Shelf temperature controller is not able to exactly follow the optimal profiles.
Once these unexpected disturbances are measured, they are introduced in the RTO scheme which computes a new optimal shelf temperature profile. Such profile allows to reach the quality constraint (water content at the end of the process) without significant increase in process time.
Figure 6e,f show the real (measured) values at the end of the process. Optimal profile has also changed with respect to the situation at 16 h to compensate for plant unexpected disturbances. However, the on-line procedure allowed to reach the desired final product water content.
Finally, we use a simulation experiment to illustrate one of the advantages of the RTO scheme as compared with the off-line optimal control approach (see Figure 7).
To that purpose, we use the experimental measurements shown in Figure 6. After the process has been running for 1 h , i.e., the first time the control profile is updated, we consider two situations:
  • We follow the RTO procedure, as in the previous case, and use experimental measurements—shelf temperature (blue segment in Figure 7a) and chamber pressure (red segment in Figure 7a)—to compute a new optimal profile (continuous black line in Figure 7a).
  • We do not use available plant measurements and force the system to follow the off-line profile previously computed (dashed black line in Figure 7a).
As shown in the figure, the RTO scheme results into 1.3 h reduction in process duration as compared with the off-line profile. Besides, the RTO scheme would allow to face other two issues [52,53,54,55]:
  • Unknown disturbances cause lower product temperatures than the predicted ones in the off-line scheme. In this case, final product water content might be larger than the targeted one, i.e., product quality would be lower than required. The RTO scheme can either increase process time or increase shelf temperature to reach the required quality.
  • Unknown disturbances cause higher product temperatures than the predicted ones in the off-line scheme. In this case, safety constraints might be not fulfilled which can cause product collapse. The RTO scheme can reduce shelf temperature to avoid product collapse.

5. Conclusions and Future Work

In this work, we derived and experimentally validated an operational model describing the behavior of the freeze-drying process. The model focuses on the time scale of the state variables related to product stability and quality, i.e., product temperature and water content. It also describes the evolution of chamber vapor pressure.
Validation results show that the model satisfactorily predicts the evolution of the relevant variables during primary and secondary drying stages. In the transition between these stages minor deviations are found. However, such deviations do not impact the prediction of the product quality. Besides, the model is conservative in the sense that temperature predictions in the transition are larger than the experimental measurements. As a result, control policies derived with the model also satisfy safety constraints in the experiments.
The validated model was subsequently used to derive optimal operation policies aiming at reducing process duration while achieving quality and stability targets. A 30% reduced process duration (13 h ) was achieved as compared to standard processes. Finally, we implemented a real-time optimization scheme in which control policies were recomputed on-line using the data measured during the process. This allowed to compensate for unexpected disturbances and model/process mismatch.
The validated operational model and the proposed RTO scheme offer new possibilities to optimize the performance of the freeze-drying process with the only requirement of recalibrating the model for a specific chamber/product.
Future research efforts will be focused on including low level controller (PID) dynamics in the computation of the optimal profiles. In this work, we assume that the low level regulator is able to perfectly follow the optimal profiles so the time required to reach the set point or possible offsets are not taken into account.

Author Contributions

Conceptualization, C.V. and A.A.A.; methodology, C.V., A.A.A., E.B.-C. and I.C.T.; software, C.V. and E.L.-Q.; validation, C.V. and I.C.T.; formal analysis, C.V., A.A.A. and I.C.T.; investigation, C.V., A.A.A., I.C.T., E.B.-C. and E.L.-Q.; resources, A.A.A.; writing–original draft preparation, C.V.; writing–review and editing, C.V., A.A.A., E.B.-C., E.L.-Q. and I.C.T.; visualization, C.V.; supervision, C.V. and A.A.A.; project administration, A.A.A.; funding acquisition, A.A.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by E.U. H2020 research and innovation programme (CoPro project, No 723575), E.U. 7th Framework Programme (CAFE Project, KBBE-2007-2-3-01) and Spanish Ministry of Science, Innovation and Universities (InCo4In project, PGC2018-099312-B-C31, MCIU/AEI/FEDER, UE). E. Lopez-Quiroga acknowledges funding received from EPSRC grant no. EP/S023070/1.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. The Landau Transform

The well-known Landau transform [34,35] is applied to the Freeze-Drying model presented in Section 2.2 to obtain an equivalent system with a fixed domain. In this section we will present the model equations after the application of the Landau transform. The following transforms are defined for the dried and frozen regions:
Q d r : ξ z ; z R | z = ξ x
Q f r : ξ y ; y R | y = ξ x L x
Note that z , y [ 0 , 1 ] . The time coordinate is not transformed although it will be renamed to avoid confusion θ R | θ = t whereas the temperature will be denoted by T , this is: T d r ( t , ξ ) T d r ( θ , z ) and T f r ( t , ξ ) T f r ( θ , y ) .
In the new system of coordinates ( z , y ), Equations (1) and (2) read as:
T d r θ = α d r x 2 2 T d r z 2 + z w x T d r z
T f r θ = α f r ( L x ) 2 2 T f r y 2 + w 1 y L x T d r y
where α d r = κ d r / ( ρ d r c p , d r ) and α f r = κ f r / ( ρ f r c p , f r ) are, respectively, the thermal diffusivity of the dried and frozen regions.
Note that a negative convection term “appears” in the model as a consequence of the application of the transformation. Moreover, two phases are required to apply Stefan condition, see Equation (3), i.e., 0 < x ( t ) < 1 .
Boundary conditions at the top and the bottom of the sample read as follows in the new coordinate system:
κ d r x T d r z z = 0 = σ e p f p ( T c h 4 T d r z = 0 4 )
κ f r L x T f r y y = 1 = h L ( T s h T f r y = 1 )
The Stefan Equation (3) reads now:
( ρ f r ρ d r ) Δ H s w = κ f r L x T f r y y = 0 κ d r x T d r z z = 1

References

  1. Andrieu, J.; Vessot, S. A review on experimental determination and optimization of physical quality factors during pharmaceuticals freeze-drying cycles. Dry. Technol. 2018, 36, 129–145. [Google Scholar] [CrossRef]
  2. Harguindeguy, M.; Fissore, D. On the effects of freeze-drying processes on the nutritional properties of foodstuff: A review. Dry. Technol. 2019. [Google Scholar] [CrossRef]
  3. Trelea, I.C.; Passot, S.; Fonseca, F.; Marin, M. An interactive tool for the optimization of freeze-drying cycles based on quality criteria. Dry. Technol. 2007, 25, 741–751. [Google Scholar] [CrossRef]
  4. Ratti, C. Freeze-Drying Process Design. In Handbook of Food Process Design; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2012; Chapter 22; pp. 621–647. [Google Scholar] [CrossRef]
  5. Malik, N.; Gouseti, O.; Bakalis, S. Effect of freezing on microstructure and reconstitution of freeze-dried high solid hydrocolloid-based systems. Food Hydrocoll. 2018, 83, 473–484. [Google Scholar] [CrossRef] [Green Version]
  6. Lopez-Quiroga, E.; Prosapio, V.; Fryer, P.J.; Norton, I.T.; Bakalis, S. Model discrimination for drying and rehydration kinetics of freeze-dried tomatoes. J. Food Process Eng. 2019, e13192. [Google Scholar] [CrossRef] [Green Version]
  7. Wang, R.; Gouseti, O.; Lopez-Quiroga, E.; Fryer, P.J.; Bakalis, S. Water Crystallization in Highly Concentrated Carbohydrate-Based Systems. Cryst. Growth Des. 2019, 19, 2081–2088. [Google Scholar] [CrossRef]
  8. Pisano, R.; Arsiccio, A.; Nakagawa, K.; Barresi, A. Tuning, measurement and prediction of the impact of freezing on product morphology: A step toward improved design of freeze-drying cycles. Dry. Technol. 2019, 37, 579–599. [Google Scholar] [CrossRef]
  9. Prosapio, V.; Norton, I.; De Marco, I. Optimization of freeze-drying using a Life Cycle Assessment approach: Strawberries’ case study. J. Clean. Prod. 2017, 168, 1171–1179. [Google Scholar] [CrossRef]
  10. Ladha-Sabur, A.; Bakalis, S.; Fryer, P.; Lopez-Quiroga, E. Mapping energy consumption in food manufacturing. Trends Food Sci. Technol. 2019, 86, 270–280. [Google Scholar] [CrossRef]
  11. Vilas, C.; Balsa-Canto, E.; García, M.S.G.; Banga, J.R.; Alonso, A.A. Dynamic optimization of distributed biological systems using robust and efficient numerical techniques. BMC Syst. Biol. 2012, 6, 79. [Google Scholar] [CrossRef] [Green Version]
  12. Balsa-Canto, E.; Alonso, A.; Arias-Méndez, A.; García, M.; López-Núñez, A.; Mosquera-Fernández, M.; Vázquez, C.; Vilas, C. Modeling and optimization techniques with applications in food processes, bio-processes and bio-systems. SEMA SIMAI Springer Ser. 2016, 9, 187–216. [Google Scholar] [CrossRef]
  13. Datta, A. Toward computer-aided food engineering: Mechanistic frameworks for evolution of product, quality and safety during processing. J. Food Eng. 2016, 176, 9–27. [Google Scholar] [CrossRef]
  14. Fissore, D.; Pisano, R.; Barresi, A.A. A model-based framework to optimize pharmaceuticals freeze drying. Dry. Technol. 2012, 30, 946–958. [Google Scholar] [CrossRef]
  15. Pisano, R.; Fissore, D.; Barresi, A.A. In-line and off-line optimization of freeze-drying cycles for pharmaceutical products. Dry. Technol. 2013, 31, 905–919. [Google Scholar] [CrossRef] [Green Version]
  16. Bosca, S.; Barresi, A.A.; Fissore, D. Fast freeze-drying cycle design and optimization using a PAT based on the measurement of product temperature. Eur. J. Pharm. Biopharm. 2013, 85, 253–262. [Google Scholar] [CrossRef] [Green Version]
  17. Bosca, S.; Barresi, A.A.; Fissore, D. On the use of model-based tools to optimize in-line a pharmaceuticals freeze-drying process. Dry. Technol. 2016, 34, 1831–1842. [Google Scholar] [CrossRef]
  18. Fissore, D. Model-Based PAT for Quality Management in Pharmaceuticals Freeze-Drying: State of the Art. Front. Bioeng. Biotechnol. 2017, 5, 5. [Google Scholar] [CrossRef] [Green Version]
  19. Shivkumar, G.; Kazarin, P.; Strongrich, A.; Alexeenko, A. LyoPRONTO: An Open-Source Lyophilization Process Optimization Tool. AAPS PharmSciTech 2019, 20. [Google Scholar] [CrossRef] [Green Version]
  20. Nakagawa, K.; Ochiai, T. A mathematical model of multi-dimensional freeze-drying for food products. J. Food Eng. 2015, 161, 55–67. [Google Scholar] [CrossRef]
  21. Warning, A.; Arquiza, J.; Datta, A. A multiphase porous medium transport model with distributed sublimation front to simulate vacuum freeze drying. Food Bioprod. Process. 2015, 94, 637–648. [Google Scholar] [CrossRef]
  22. Aydin, E.; Yucel, O.; Sadikoglu, H. Modelling and simulation of a moving interface problem: Freeze drying of black tea extract. Heat Mass Transf./Waerme- Stoffuebertragung 2017, 53, 2143–2154. [Google Scholar] [CrossRef]
  23. Tarafdar, A.; Shahi, N. Application and comparison of genetic and mathematical optimizers for freeze-drying of mushrooms. J. Food Sci. Technol. 2018, 55, 2945–2954. [Google Scholar] [CrossRef] [PubMed]
  24. Antelo, L.T.; Passot, S.; Fonseca, F.; Trelea, I.C.; Alonso, A.A. Toward Optimal Operation Conditions of Freeze-Drying Processes via a Multilevel Approach. Dry. Technol. 2012, 30, 1432–1448. [Google Scholar] [CrossRef] [Green Version]
  25. Fissore, D.; Pisano, R.; Barresi, A.A. Applying quality-by-design to develop a coffee freeze-drying process. J. Food Eng. 2014, 123, 179–187. [Google Scholar] [CrossRef] [Green Version]
  26. Lopez-Quiroga, E.; Antelo, L.T.; Alonso, A.A. Time-scale modeling and optimal control of freeze-drying. J. Food Eng. 2012, 111, 655–666. [Google Scholar] [CrossRef] [Green Version]
  27. Jennings, T.A. Lyophilization: Introduction and Basic Principles; Interpharm/CRC: Boca Raton, FL, USA, 1999. [Google Scholar]
  28. Saguy, I.S. Challenges and opportunities in food engineering: Modeling, virtualization, open innovation and social responsibility. J. Food Eng. 2016, 176, 2–8. [Google Scholar] [CrossRef]
  29. Passot, S.; Cenard, S.; Douania, I.; Trelea, I.C.; Fonseca, F. Critical water activity and amorphous state for optimal preservation of lyophilized lactic acid bacteria. Food Chem. 2012, 132, 1699–1705. [Google Scholar] [CrossRef]
  30. Song, C.; Nam, J.H.; Kim, C.J. A finite volume analysis of vacuum freeze-drying processes of skim milk solution in trays and vials. Dry. Technol. 2002, 20, 283–305. [Google Scholar] [CrossRef]
  31. Litchfield, R.J.; Liapis, A.I. An adsorption-sublimation model for a freeze-dryer. Chem. Eng. Sci. 1979, 34, 1085–1090. [Google Scholar] [CrossRef]
  32. Mascarenhas, W.J.; Akay, H.U.; Pikal, M.J. A computational model for finite element analysis of the freeze-drying process. Comput. Methods Appl. Mech. Eng. 1997, 148, 105–124. [Google Scholar] [CrossRef]
  33. Pikal, M.J. Transport Processes in Pharmaceutical Systems. Drugs and The Pharmaceutical Sciences. In Heat and Mass Transfer in Low Pressure Gases: Applications to Freeze Drying; Marcel Dekker: New York, NY, USA, 2000; Volume 102, Chapter 16; pp. 611–686. [Google Scholar]
  34. Landau, H.G. Heat conduction in a melting solid. Q. Appl. Math. 1950, 8, 81–94. [Google Scholar] [CrossRef] [Green Version]
  35. Illingworth, T.C.; Golosnoy, I.O. Numerical solutions of diffusion–controlled moving boundary problems which conserve solute. J. Comput. Phys. 2005, 209, 207–225. [Google Scholar] [CrossRef] [Green Version]
  36. Trelea, I.C.; Fonseca, F.; Passot, S.; Flick, D. A Binary Gas Transport Model Improves the Prediction of Mass Transfer in Freeze Drying. Dry. Technol. 2015, 33, 1849–1858. [Google Scholar] [CrossRef]
  37. Nam, J.H.; Song, C.S. Numerical simulation of conjugate heat and mass transfer during multi-dimensional freeze drying of slab-shaped food products. Int. J. Heat Mass Transf. 2007, 50, 4891–4900. [Google Scholar] [CrossRef]
  38. Millman, M.J.; Liapis, A.I.; Marchello, J.M. Guidelines for the desirable operation of batch freeze dryers during removal of free water. J. Food Technol. 1984, 19, 725–738. [Google Scholar] [CrossRef]
  39. Boss, E.A.; Filho, R.M.; de Toledo, E.C.V. Freeze-drying process: Real-time model and optimization. Chem. Eng. Process. 2004, 43, 1475–1485. [Google Scholar] [CrossRef]
  40. Pikal, M.J.; Shah, S.; Roy, M.L.; Putman, R. The secondary drying stage of freeze drying: Drying kinetics as a function of temperature and chamber pressure. Int. J. Pharm. 1990, 60, 203–217. [Google Scholar] [CrossRef]
  41. Pisano, R.; Fissore, D.; Barresi, A.A. Quality by Design in the Secondary Drying Step of a Freeze-Drying Process. Dry. Technol. 2012, 30, 1307–1316. [Google Scholar] [CrossRef]
  42. Trelea, I.C.; Fonseca, F.; Passot, S. Dynamic modeling of the secondary drying stage of freeze drying reveals distinct desorption kinetics for bound water. Dry. Technol. 2016, 34, 335–345. [Google Scholar] [CrossRef]
  43. Pikal, M.J.; Shah, S. The collapse temperature in freeze drying: Dependence on measurement methodology and rate of water removal from the glassy phase. Int. J. Pharm. 1990, 62, 165–186. [Google Scholar] [CrossRef]
  44. Wilhelm, L.R. Numerical calculation of psychrometric properties in SI units. Trans. ASAE 1976, 19, 318–325. [Google Scholar] [CrossRef]
  45. Egea, J.A.; Vazquez, E.; Banga, J.R.; Marti, R. Improved scatter search for the global optimization of computationally expensive dynamic models. J. Glob. Opt. 2009, 43, 175–190. [Google Scholar] [CrossRef] [Green Version]
  46. Balsa-Canto, E.; Banga, J.R. AMIGO, a toolbox for advanced model identification in systems biology using global optimization. Bioinformatics 2011, 27, 2311–2313. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  47. Ljung, L. System Identification: Theory for the User, 2nd ed.; PTR Prentice Hall: Prentice-Hall, NJ, USA, 1999. [Google Scholar]
  48. Vassiliadis, V.S.; Sargent, R.W.H.; Pantelides, C.C. Solution of a Class of Multistage Dynamic Optimization Problems .2. Problems with Path Constraints. Ind. Eng. Chem. Res. 1994, 33, 2123–2133. [Google Scholar] [CrossRef]
  49. Vassiliadis, V.S.; Sargent, R.W.H.; Pantelides, C.C. Solution of a Class of Multistage Dynamic Optimization Problems. 1. Problems Without Path Constraints. Ind. Eng. Chem. Res. 1994, 33, 2111–2122. [Google Scholar] [CrossRef]
  50. Balsa-Canto, E.; Banga, J.R.; Alonso, A.A.; Vassiliadis, V. Dynamic optimization of distributed parameter systems using second-order directional derivatives. Ind. Eng. Chem. Res. 2004, 43, 6756–6765. [Google Scholar] [CrossRef]
  51. Banga, J.R.; Seider, W.D. Global optimization of chemical processes using stochastic algorithms. State of the Art in Global Optimization—Computational Methods and Applications. In Nonconvex Optimization and Its Applications; Floudas, C.A., Pardalos, P.M., Eds.; Kluwer Academic Publ.: Dordrecht, The Netherlands, 1996; Volume 7, pp. 563–583. [Google Scholar]
  52. Akterian, S.G. On-line control strategy for compensating for arbitrary deviations in heating-medium temperature during batch thermal sterilization processes. J. Food Eng. 1999, 39, 1–7. [Google Scholar] [CrossRef]
  53. Simpson, R.; Figeroa, I.; Teixeira, A. Simple, practical, and efficient on-line correction of process deviations in batch retort through simulation. Food Control 2007, 18, 458–465. [Google Scholar] [CrossRef]
  54. Alonso, A.A.; Arias-Méndez, A.; Balsa-Canto, E.; García, M.R.; Molina, J.I.; Vilas, C.; Villafín, M. Real time optimization for quality control of batch thermal sterilization of prepackaged foods. Food Control 2013, 32, 392–403. [Google Scholar] [CrossRef] [Green Version]
  55. Marchetti, A.G.; François, G.; Faulwasser, T.; Bonvin, D. Modifier Adaptation for Real-Time Optimization-Methods and Applications. Processes 2016, 4, 55. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Schematic representation of the lyophilization pilot plant.
Figure 1. Schematic representation of the lyophilization pilot plant.
Processes 08 00325 g001
Figure 2. Schematic representation of the freeze drying process. Two subsystems are considered, the condenser and the product. The arrows entering the subsystems represent the process control variables whereas output arrows represent the state variables.
Figure 2. Schematic representation of the freeze drying process. Two subsystems are considered, the condenser and the product. The arrows entering the subsystems represent the process control variables whereas output arrows represent the state variables.
Processes 08 00325 g002
Figure 3. Predictive capabilities of the freeze-drying model from the product temperature point of view. Blue line corresponds to the shelf temperature. Continuous black line and black dots represent, respectively, the model predictions and experimental measurements at the bottom of the product. Dashed black lines correspond to the model prediction at other spatial points in the product. Vertical lines represent the end of primary drying and the beginning of secondary drying.
Figure 3. Predictive capabilities of the freeze-drying model from the product temperature point of view. Blue line corresponds to the shelf temperature. Continuous black line and black dots represent, respectively, the model predictions and experimental measurements at the bottom of the product. Dashed black lines correspond to the model prediction at other spatial points in the product. Vertical lines represent the end of primary drying and the beginning of secondary drying.
Processes 08 00325 g003
Figure 4. Predictive capabilities of the freeze-drying model for the chamber vapor pressure. Black circles represent the chamber vapor pressure obtained from measurements of Pirani and total chamber pressure whereas continuous black lines correspond to the model predictions.
Figure 4. Predictive capabilities of the freeze-drying model for the chamber vapor pressure. Black circles represent the chamber vapor pressure obtained from measurements of Pirani and total chamber pressure whereas continuous black lines correspond to the model predictions.
Processes 08 00325 g004
Figure 5. Evolution of the product temperature and shelf temperature (top figures) as well as water content (bottom figures). Figures on the left correspond with experiment 1 whereas the optimal design results are represented on the figures on the right. Horizontal line on the bottom figures represent the constraint on the final water content.
Figure 5. Evolution of the product temperature and shelf temperature (top figures) as well as water content (bottom figures). Figures on the left correspond with experiment 1 whereas the optimal design results are represented on the figures on the right. Horizontal line on the bottom figures represent the constraint on the final water content.
Processes 08 00325 g005
Figure 6. Offline and real time optimization (RTO) scheme results. Subfigures (a,c,e) represent the evolution of shelf temperature (decision variable) and chamber pressure. Black and gray parts of the profiles correspond with the optimal solution. Blue and red parts correspond with the plant measurements. Subfigures (b,d,f) represent the evolution of product mean water content. Dashed line indicates the quality objective to reach.
Figure 6. Offline and real time optimization (RTO) scheme results. Subfigures (a,c,e) represent the evolution of shelf temperature (decision variable) and chamber pressure. Black and gray parts of the profiles correspond with the optimal solution. Blue and red parts correspond with the plant measurements. Subfigures (b,d,f) represent the evolution of product mean water content. Dashed line indicates the quality objective to reach.
Processes 08 00325 g006
Figure 7. Comparison, in terms of process duration, between the RTO approach (continuous black lines) and the off-line control profile (dashed black lines). (a) Shelf temperature and chamber pressure. (b) Mean water content evolution. Blue and red segments in subplot (a) correspond with the real plant measurements for shelf temperature and chamber pressure, respectively.
Figure 7. Comparison, in terms of process duration, between the RTO approach (continuous black lines) and the off-line control profile (dashed black lines). (a) Shelf temperature and chamber pressure. (b) Mean water content evolution. Blue and red segments in subplot (a) correspond with the real plant measurements for shelf temperature and chamber pressure, respectively.
Processes 08 00325 g007
Table 1. Parameters involved in the freeze-drying model. Parameters to be estimated are indicated in the “Value” column as t.b.e.
Table 1. Parameters involved in the freeze-drying model. Parameters to be estimated are indicated in the “Value” column as t.b.e.
ParameterValueUnitsDescription
ρ d r 200.31 k g   m 3 Dried region density
ρ f r 1001.6 k g   m 3 Frozen region density
c p , d r 1254 J   k g 1   K 1 Dried region heat capacity
c p , f r 1818.8 J   k g 1   K 1 Frozen region heat capacity
κ d r t.b.e. W   m 1   K 1 Dried region heat conductivity
κ f r 2.4 W   m 1   K 1 Frozen region heat conductivity
σ 5.6704 × 10−8 W   m 2   K 4 Stefan–Boltzmann constant
e p 0.78-Thermal emissivity at the product top
f p 0.99-Geometrical correction factor
K c l a p 1.6548 × 10−4 K 1 Constant in the Clapeyron equation
Δ H s 2791.2 × 10−3 J   k g 1 Sublimation heat
R8314 Pa m 3   K 1   kmol 1 Ideal gas constant
L x 5.75 × 10−3 m Food product height
L z 0.242 m Food product length
L y 0.307 m Food product width
M w 18 k g   koml 1 Water molecular mass
h L , 1 3.3 W   m 2   K 1 Heat transfer coefficient constant
h L , 2 t.b.e. Pa 1 Heat transfer coefficient constant
h L , 3 34.4 Pa Heat transfer coefficient constant
k 1 430.0 s   m   k g 1 Mass transfer coefficient constant in Darcy’s equation
k 2 t.b.e. s   Pa   k g 1 Mass transfer coefficient constant in Darcy’s equation
β t.b.e. k g   s 1   K 1 Mass transfer coefficient constant in
chamber/condenser flux
T c h 293.15 K Chamber temperature
V c h 0.202 m 3 Chamber volume
τ A r e f 2.689 × 104 s Compartment A reference time constant
τ B r e f 6.493 × 105 s Compartment B reference time constant
E a 4.271 × 104 k J   k g 1 Activation energy in desorption model
T r e f 273.15 K Reference temperature in desorption model
α A 0.669-Compartment A ratio between equilibrium water content
α B 0.331-Compartment B ratio between equilibrium water content
M g 0.0434 k g -water   k g 1 -totalConstant of the GAB equation
C g 7.4789-Constant of the GAB equation
K g 0.9827-Constant of the GAB equation
K T , g 8.2-Constant of the glass transition temperature
T g , l −135°CConstant of the glass transition temperature
T g , s 75.58°CConstant of the glass transition temperature
α 1.6-Ratio of molecular heat conductivities of vapor
and nitrogen
Table 2. Cross validation results for the product model. Five estimations of the parameters were performed, each one excluding a different experiment. R M S E e and R M S E v denote, respectively, the Root Mean Square Error (RMSE) values obtained for the estimation and validation experiments.
Table 2. Cross validation results for the product model. Five estimations of the parameters were performed, each one excluding a different experiment. R M S E e and R M S E v denote, respectively, the Root Mean Square Error (RMSE) values obtained for the estimation and validation experiments.
ParameterExcluded Experiment
12345
h L , 2 1.961.861.931.822.03
k 2 8.60 × 1077.86 × 1076.52 × 1079.41 × 1079.51 × 107
κ d r 2.45 × 10−44.13 × 10−42.41 × 10−42.85 × 10−42.85 × 10−4
R M S E e 2.151.642.112.102.00
R M S E v 1.623.102.012.492.54
Table 3. Cross validation results for the product/condenser model. Five estimations of parameter β were carried out, each one excluding a different experiment. R M S E e and R M S E v denote, respectively, the RMSE values obtained for the estimation and validation experiments.
Table 3. Cross validation results for the product/condenser model. Five estimations of parameter β were carried out, each one excluding a different experiment. R M S E e and R M S E v denote, respectively, the RMSE values obtained for the estimation and validation experiments.
ParameterExcluded Experiment
12345
β 2.31 × 1032.17 × 1032.31 × 1032.37 × 1032.37 × 103
R M S E e 1.61.72.22.01.8
R M S E v 2.82.10.60.62.3

Share and Cite

MDPI and ACS Style

Vilas, C.; A. Alonso, A.; Balsa-Canto, E.; López-Quiroga, E.; Trelea, I.C. Model-Based Real Time Operation of the Freeze-Drying Process. Processes 2020, 8, 325. https://doi.org/10.3390/pr8030325

AMA Style

Vilas C, A. Alonso A, Balsa-Canto E, López-Quiroga E, Trelea IC. Model-Based Real Time Operation of the Freeze-Drying Process. Processes. 2020; 8(3):325. https://doi.org/10.3390/pr8030325

Chicago/Turabian Style

Vilas, Carlos, Antonio A. Alonso, Eva Balsa-Canto, Estefanía López-Quiroga, and Ioan Cristian Trelea. 2020. "Model-Based Real Time Operation of the Freeze-Drying Process" Processes 8, no. 3: 325. https://doi.org/10.3390/pr8030325

APA Style

Vilas, C., A. Alonso, A., Balsa-Canto, E., López-Quiroga, E., & Trelea, I. C. (2020). Model-Based Real Time Operation of the Freeze-Drying Process. Processes, 8(3), 325. https://doi.org/10.3390/pr8030325

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop