Optimal boundary control of a contact thawing process for foodstuff

11th IFAC Symposium on Dynamics and Control of 11th Symposium on and Process including Biosystems 11th IFAC IFACSystems, Symposium on Dynamics Dynamics and Control Control of of 11th IFAC Symposium on Dynamics and Control of Process Systems, including Biosystems Available online at www.sciencedirect.com June 6-8, 2016. NTNU, Trondheim, Norway Process Systems, including Biosystems Process Systems, including Biosystems June 6-8, 2016. NTNU, Trondheim, Norway June June 6-8, 6-8, 2016. 2016. NTNU, NTNU, Trondheim, Trondheim, Norway Norway

ScienceDirect

IFAC-PapersOnLine 49-7 (2016) 183–188

Optimal boundary control of a contact thawing Optimal control of Optimal boundary control of a contact thawing Optimal boundary boundary control of a a contact contact thawing thawing process for foodstuff process for foodstuff process for foodstuff process for foodstuff

Christoph Josef Backi ∗,1 John Leth ∗∗ Jan Tommy Gravdahl ∗∗∗ ∗,1 ∗∗ ∗,1 ∗∗ ∗∗∗ Christoph Josef Backi John Leth Jan Tommy Gravdahl ∗∗∗ Christoph Christoph Josef Josef Backi Backi ∗,1 John John Leth Leth ∗∗ Jan Jan Tommy Tommy Gravdahl Gravdahl ∗∗∗ ∗ Department of Chemical Engineering, Norwegian University of Science and ∗∗ Department of Chemical Engineering, Norwegian University of Science and ∗Technology, DepartmentN-7491 of Chemical Chemical Engineering, Norwegian University of of Science Science and and Trondheim, NorwayNorwegian (e-mail: [email protected]). Department of Engineering, University Technology, N-7491 Trondheim, Norway (e-mail: [email protected]). ∗∗ Department Technology, N-7491 Trondheim, Norway (e-mail: [email protected]). of Electronic Systems, Aalborg University, DK-7220, Aalborg, Technology, N-7491 Trondheim, Norway (e-mail: [email protected]). ∗∗ ∗∗ of Electronic Systems, Aalborg University, DK-7220, Aalborg, ∗∗ Department Department Systems, Aalborg University, Denmark (e-mail: [email protected]) Department of of Electronic Electronic Systems, Aalborg University, DK-7220, DK-7220, Aalborg, Aalborg, Denmark (e-mail: [email protected]) ∗∗∗ Department of Engineering Denmark (e-mail: [email protected]) Cybernetics, Norwegian University of Science Denmark (e-mail: [email protected]) ∗∗∗ ∗∗∗ Department Engineering Cybernetics, Norwegian ∗∗∗ Department of Engineering Cybernetics, Norwegian University of Science Science and of Technology, N-7491 Trondheim, NorwayUniversity (e-mail: of Department of Engineering Cybernetics, Norwegian University of Science and Technology, N-7491 Trondheim, Norway (e-mail: and Technology, N-7491 Trondheim, Norway (e-mail: [email protected]). and Technology, N-7491 Trondheim, Norway (e-mail: [email protected]). [email protected]). [email protected]). Abstract: In this work an approach for thawing blocks of foodstuff, in particular fish, is introduced. The Abstract: this work an approach for thawing blocks of foodstuff, particular is introduced. The Abstract: In this an for thawing blocks in particular fish, is The functional In principle is based on plate freezer technology, which hasin been used infish, industry for decades. Abstract: In this work work an approach approach forfreezer thawing blocks of of foodstuff, foodstuff, inbeen particular fish, is introduced. introduced. The functional principle is based on plate technology, which has used in industry for decades. functional principle is based based on plate plate the freezer technology, which has has beenthawing used in in process industryby formeans decades. The aim ofprinciple this work is to describe temperature dynamics of this of functional is on freezer technology, which been used industry for decades. The aim of this work is to describe the temperature dynamics of this thawing process by means of The aim of this work is to describe the temperature dynamics of this thawing process by means of partial differential equations (PDEs) and control the boundary conditions in an optimal way. The PDE The aim of this work is to describe the temperature dynamics of this thawing process by means of partial differential equations (PDEs) and control the diffusion boundary conditions in an optimal way. The PDE partial differential equations (PDEs) control boundary conditions in an The describing the temperature dynamics is based on the equation with state-dependent parameter partial differential equationsdynamics (PDEs) and and control the diffusion boundaryequation conditions in state-dependent an optimal optimal way. way. parameter The PDE PDE describing the temperature is based on the with describing the temperature temperature dynamics dynamics is is based based on on the the diffusion diffusion equation equation with with state-dependent state-dependent parameter parameter functions. the describing functions. functions. functions. © 2016, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Thawing, Foodstuff, Partial Differential Equations, Optimal Control, MPC Keywords: Thawing, Foodstuff, Partial Differential Equations, Optimal Control, MPC Keywords: Keywords: Thawing, Thawing, Foodstuff, Foodstuff, Partial Partial Differential Differential Equations, Equations, Optimal Optimal Control, Control, MPC MPC capacity leads to increased drip loss, which results in lower 1. INTRODUCTION capacity leads increased drip loss, results in in lower lower 1. 1. INTRODUCTION INTRODUCTION capacity weight. leads to toDrying increased dripsurface loss, which which product of the also results capacity leads to increased drip loss, which results in in lower lower 1. INTRODUCTION weight. Drying of the surface also results weight. Drying of the surface also results in lower product weight, whereas lipid oxidation occurs mostly for fatty Shelf life extension of rapidly spoiling foodstuff is an essential product product weight. Drying of the surface also results in lower product weight, whereas lipid oxidation occurs mostly for fatty Shelf life extension of rapidly spoiling foodstuff is an essential product weight, whereas lipid oxidation occurs mostly for species and can cause rancid tastes and smells. Both, drying and Shelf life extension of rapidly spoiling foodstuff is an essential task for the food industry. This holds especially for fish and product weight, whereas lipidtastes oxidation occursBoth, mostly for fatty fatty Shelf lifethe extension of rapidly spoiling foodstuff is for an essential species and can cause rancid and smells. drying and task for food industry. This holds especially fish and species and can cause rancid tastes and smells. Both, drying and lipid oxidation are most likely to happen in air-based thawing task for the food industry. This holds especially for fish and fish products as they often have to be transported over long species and can cause rancid tastes and smells. Both, drying and task for the food industry. This holds especially forover fish long and lipid oxidation are most likely to happen in air-based thawing fish products as they often have to be transported lipid oxidation are most likely to happen in air-based thawing methods. fish products products as they they often have to in beorder transported over long distances or stored overoften long have periods to deliver them to lipid oxidation are most likely to happen in air-based thawing fish as to be transported over long distances or long periods order deliver them to methods. distances or stored stored over longprocessing periods in instages. order to toBy deliver themthe to methods. the consumer or toover further lowering methods. distances or stored over long periods in order to deliver them to There is a whole range of publications in the scientific comthe consumer or to further processing stages. By lowering the the consumer or to further processing stages. By lowering the temperature, the rate of spoilage is reduced. As a general rule it There aa whole range of publications the scientific comthe consumer or to further processing stages. By lowering the There is of in the communityis different methods forin and thawing temperature, the rate of spoilage is reduced. As a general rule it isintroducing a whole whole range range of publications publications infreezing the scientific scientific comtemperature, thelower rate of of spoilage is reduced. reduced. As aa general general rule it There holds that, thethe thespoilage temperature, the longer the shelfrule life. munity introducing different methods for freezing and thawing temperature, rate is As it munity introducing different methods for freezing and thawing of foodstuff, see e.g. Johnston et al. (1994) and Jason (1974). holds that, the lower the temperature, the longer the shelf life. munity introducing different methods for freezing and thawing holds that, that, the the lower loweristhe the temperature, the longer the shelf shelf life. of foodstuff, see e.g. Johnston et al. (1994) and Jason (1974). Therefore, freezing thetemperature, best methodthe forlonger shelf life extension holds the life. of foodstuff, e.g. al. (1994) and (1974). Moreover, thesee modeling of heatet phenomena Therefore, freezing foodstuff, see e.g. Johnston Johnston etand al. mass (1994)transfer and Jason Jason (1974). Therefore, freezing is is the the best best method method for for shelf shelf life life extension extension of over long periods. Moreover, of heat mass transfer phenomena Therefore, freezing is the best method for shelf life extension Moreover, the modeling of been heat and and masse.g. transfer phenomena with PDEs the hasmodeling extensively studied, in Pham (2006b) over long periods. Moreover, the modeling of heat and mass transfer phenomena over long periods. with PDEs has extensively been studied, e.g. in Pham (2006b) over periods. with PDEs has been studied, e.g. in Pham In Cleland al. (1987) data Afterlong having been transported over long distances or stored over and with PDEs(2006a). has extensively extensively beenet studied, e.g.experimental in Pham Pham (2006b) (2006b) and Pham (2006a). In Cleland et al. (1987) experimental data After having been transported over long distances or stored over and Pham (2006a). In Cleland et al. (1987) experimental data for freezing and thawing of multi-dimensional objects modeled After having been transported over long distances or stored over long periods, the foodstuff might not necessarily be sold in and Pham (2006a). In Cleland et al. (1987) experimental data After having been transported over long distances or stored over for freezing and thawing of multi-dimensional objects modeled long periods, the foodstuff might not necessarily be sold in for freezing and thawing of multi-dimensional objects modeled by finite element techniques are presented. As computational long periods, the foodstuff might not necessarily be sold in frozen state, but processed further. For example, frozen fillets for freezing and thawing of multi-dimensional objects modeled long periods, the foodstuff might not necessarily be sold in element computational frozen further. example, frozen by finite finitehas element techniques are presented. presented. Assimulations computational power growntechniques in recent are years, complexAs can frozen state, but processed further. For example, frozen fillets in blockstate, formbut canprocessed be cut more or lessFor directly to fish sticksfillets after by by finite element techniques are presented. As computational frozen state, but processed further. For example, frozen fillets power has grown in recent years, complex simulations can in block form can be cut more or less directly to fish sticks after power has grown in recent years, complex simulations can provide qualitative and quantitative insights to the system, even blocktempering form can can be be cut more more or less less(1974). directlyOn to the fishother stickshand, after power has grown in recent years, complex simulations ainshort phase, see Jason can in block form cut or directly to fish sticks after provide qualitative and quantitative insights to the system, even afilleting short tempering phase, see Jason (1974). On the other hand, provide qualitative and quantitative insights to the system, even in the absence of analytical solutions. An open-loop approach aa short short tempering tempering phase, see Jason (1974). On the other hand, of headed phase, and gutted fish is not possible in frozen state. provide qualitative and quantitative insights to the system, even see Jason (1974). On the other hand, the absence of analytical solutions. An open-loop approach filleting gutted is in state. in the absence of solutions. An open-loop approach for optimal boundary control of a plate process has filleting of of headed headed and gutted fish fish is not not possible in frozen frozen state. in Therefore, thawingand techniques must be possible used in order to enable in the absence of analytical analytical solutions. Anfreezing open-loop approach filleting of headed and gutted fish is not possible in frozen state. for optimal boundary control of a plate freezing process has Therefore, thawing techniques must be used in order to enable for optimal optimal boundary control of aa plate plate freezing process has been described in Backicontrol and Gravdahl (2013). However, it must Therefore, thawing Thus, techniques mustfreezing be used usedhas in order order to impact enable for further processing. not only a large boundary of freezing process has Therefore, thawing techniques must be in to enable been described in Backi and Gravdahl (2013). However, it must further processing. Thus, not only freezing has a large impact been described in Backi and Gravdahl (2013). However, it must be mentioned that thawing is not simply the reversed freezing further processing. Thus, not only freezing has a large impact on the quality of theThus, final not product, but also has thawing. Thawing been described in Backi and Gravdahl (2013). However, it must further processing. only freezing a large impact mentioned thawing is not in simply reversed freezing on the final but thawing. Thawing be mentioned that thawing is simply the reversed freezing process, as forthat example outlined Olverthe (2014, Section 4.1). on the the quality ofthe thequality final product, product, but also also thawing. Thawing can notquality enhanceof of the frozen product, meaning that be be mentioned that thawing is not not in simply the reversed freezing on the quality of the final product, but also thawing. Thawing process, as for example outlined Olver (2014, Section 4.1). can not enhance the quality of the frozen product, meaning that process, as for example outlined in Olver (2014, Section 4.1). There it is mentioned that the backwards heat equation, where can not enhance the quality of the frozen product, meaning that quality losses and damages caused in earlier stages cannot be process, as for example outlined in Olver (2014, Section 4.1). can not losses enhance thedamages quality of the frozen product, meaning that There it is mentioned that the backwards heat equation, where quality and caused in earlier stages cannot be There it it is past mentioned that the the backwards heat tequation, equation, where future and are changed bybackwards changing from to −t, is where an illquality losses losses andthe damages caused in in earlier earlier stages cannot be There reversed by even most sophisticated way ofstages thawing. Hence is mentioned that heat quality and damages caused cannot be and past in are changed by changing from tt to −t, is an illreversed even the sophisticated way of Hence future and are changed by from −t, an illposed problem sense that small changes initial data reversed by even the most mostto sophisticated wayprocess of thawing. thawing. Hence future it is evenby more important understand the of thawing, future and past past in arethe changed by changing changing fromin t to to −t, is is an can illreversed by even the most sophisticated way of thawing. Hence posed problem the sense that small changes in initial data can it is even more important to understand the process of thawing, posed problem in the sense that small changes in initial data can lead to randomly large changes in the solution arbitrarily close it is even more important to understand the process of thawing, come up with new, gentle methods, derive mathematical models posed problem in the sense that small changes in initial data can it is even more important to understand the process of thawing, to come gentle derive mathematical models lead to randomly randomly large changes changes in in the the solution solution arbitrarily arbitrarily close close to initial time. large comeinup up with new, gentle methods, derive mathematical models lead and thewith endnew, control the methods, process such that an optimal solution lead to randomly large changes in the solution arbitrarily close come up with new, gentle methods, derive mathematical models to initial time. and in the end control the process such that an optimal solution to initial time. andobtained, in the the end end controlthat the process process such that that an optimal optimal solution is meaning quality losses in the thawing solution process to initial time. and in control the such an Applied methods for thawing of fish blocks in industry are is meaning is obtained, obtained, meaning that that quality quality losses losses in in the the thawing thawing process process Applied methods for thawing of fish blocks in industry are can be minimized. is obtained, meaning that quality losses in the thawing process Appliedbased methods forconduction thawing of ofand fishconvection blocks in inbyindustry industry are mostly on heat direct concan be minimized. Applied methods for thawing fish blocks are can be be minimized. minimized. based on heat conduction and convection direct concan mostly based on heat conduction and convection by direct contact with the thawing medium. For this purpose by mostly water As quality is influenced by many factors and parameters, a mostly mostly based on heat conduction and convection by direct conthawing medium. For this purpose mostly water As influenced many factors and tact with the thawing medium. this mostly water or airwith are the used. In water-based techniques the fish blocks are As quality quality is definition influencedforby bythe many factors and parameters, parameters, more preciseis process of thawing is needed.aaa tact tact with the thawing medium. For For this purpose purpose mostly water As quality is influenced by many factors and parameters, or air are used. In water-based techniques the fish blocks are more precise definition for the process of thawing is needed. or air are used. In water-based techniques the fish blocks are put into a water bath. Added air bubbles enhance heat transfer. more precise definition for the process of thawing is needed. The most critical quality parameters influenced by the thawing air area water used. bath. In water-based techniques the fish blocks are more precise definition for the process of thawing is thawing needed. or put into Added air bubbles enhance heat transfer. The most critical quality parameters influenced by the put into a water bath. Added air bubbles enhance heat transfer. In addition the water is sometimes heated in order to speed up The most critical quality parameters influenced by the thawing process are water holding capacity, drying of the surface and put into a water bath. Added air bubbles enhance heat transfer. The most critical quality parameters influenced by the thawing In addition the water method is sometimes sometimes heated intransfer order to tobyspeed speed up process are holding drying of surface and water is order up theaddition process.the Another relies heated on heatin moving process are water waterNot holding capacity, drying the of the the surface and In lipid oxidation. only capacity, do they lower overall quality In addition the water is sometimes heated in order to speed up process are water holding capacity, drying of the surface and the process. Another method relies on heat transfer by moving lipid oxidation. Not only do they lower the overall quality the process. Another method relies on heat transfer by moving air, known as air-blast thawing. Due to the danger of lipid lipid oxidation. Not only do they lower the overall quality directly, but also the product’s yield. Reduced water holding the process. Another method relies on heat transfer by moving lipid oxidation. Not only do they lower the overall quality known as air-blast thawing. to the of lipid directly, but the air, known as air-blast thawing. Due to danger of oxidation of the surfaceDue it should be danger conducted with directly, but also also the product’s product’s yield. yield. Reduced Reduced water water holding holding air, 1directly, air, knownand asdrying air-blast thawing. Due to the the danger of lipid lipid but also corresponding authorthe product’s yield. Reduced water holding oxidation and drying of the surface it should be conducted 1 corresponding author oxidation and and drying drying of of the the surface surface it it should should be be conducted conducted with with 1 oxidation with

1 corresponding corresponding author author Copyright © 2016, 2016 IFAC 183 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © IFAC (International Federation of Automatic Control) Copyright © 2016 IFAC 183 Peer review under responsibility of International Federation of Automatic Copyright © 2016 IFAC 183 Copyright © 2016 IFAC 183Control. 10.1016/j.ifacol.2016.07.243

IFAC DYCOPS-CAB, 2016 184 Christoph Josef Backi et al. / IFAC-PapersOnLine 49-7 (2016) 183–188 June 6-8, 2016. NTNU, Trondheim, Norway

humidified air. The two mentioned methods are only controllable to a small extend as the physical processes underlying heat transfer are complex. The mentioned disadvantages motivate the introduction of another thawing method that is not as widely spread in industry as the two mentioned before. It does not rely on heat transfer by direct contact with the thawing medium, but on principles already known from plate freezer technology. Hereby a refrigerant flows through the plate freezer walls, where it vaporizes and extracts heat from the fish block. The thawing approach investigated here relies on the same principle, but instead of using a refrigerant, a thawing medium with high specific heat capacity can be used to thaw the fish block. The thawing medium, however, does not undergo phase change during the process. To the authors best knowledge no work has been done so far in modeling the temperature dynamics of a contact thawing process. Also the important matter of (optimal) control of such a process has not been covered. Therefore this work presents a novelty with respect to the field of application. The rest of the paper is organized as follows. Section 2 introduces the model and addresses some issues regarding the application of this model to both, freezing and thawing processes. Section 3 describes the problem setting for Model Predictive Control, whereas Section 4 shows numerical simulation examples that highlight the results in the previous section. Finally, Section 5 provides some concluding remarks.

2.1 System structure In Figure 1 we propose a control structure for the thawing process consisting of a cold water feed, a hot water buffer tank, two pumps and three (controlled) valves. The idea is to use an algorithm based upon Model Predictive Control (MPC) controlling the unit consisting of the two pumps and the valve mixing the hot and the cold water such that the desired boundary temperature in the thawer can be adjusted. The return water from the thawer can be split up to be either fed back to the hot water buffer tank or be disposed via the outlet (parts of it could also be fed back to the cold water feed, for example). In the present work only the temperature dynamics of the thawer will be considered. In addition, we assume that the required temperature of the thawing medium can be exactly provided by the valve combining the hot and cold water feed. It is explicitly not intended to provide an underlying control structure for the unit providing the desired water (and thus boundary) temperature to the thawer. This is a task for future work and can result in energy-efficient operation.   







2. PROBLEM FORMULATION

 

The process of thawing foodstuff is often performed based on thermodynamic steady state calculations (prior to execution), meaning that the temperature dynamics are neglected. In order to obtain better monitoring and in the end control the quality during thawing, the temperature dynamics play a vital role. However, these temperature dynamics and the temperature distribution throughout the foodstuff are hard to model in standard ways of thawing, such as immersion or air-blast processes. Thus we propose a contact thawing process in order to monitor and control the foodstuff’s interior temperature field, which directly influences the quality. The choice of contact thawing is justified by the fact that the related process of freezing in plate freezers is well-understood and its results in modeling can be directly applied to the contact thawing process. Existing monitoringalgorithms to estimate the temperature field inside the foodstuff can be adapted to the contact thawing application, see e.g. Backi et al. (2015). The overall aim of this work is to provide an optimal boundary control formulation for a contact thawing process such that determined temperature margins can be satisfied during the whole operation. This is beneficial for the final product quality. It must be mentioned here that the time it takes to thaw foodstuff is significantly larger than to freeze it. The reason for this is that ice is a better conductor of heat compared to water. This means that for freezing the heat flow from the inner domain towards the boundary layer is accelerated over time due to ice layer formation at the boundary propagating towards the interior. For thawing, however, an initial liquid layer builds up at the boundary and propagates towards the interior causing a decrease in heat flow from the boundary to the inside of the food over time. 184

 



Fig. 1. Functional principle of the contact thawing process

2.2 Temperature dynamics model in the thawer To model the temperature dynamics of the food in the thawer we use the diffusion equation with state-dependent parameters in the variable T = T (t, x), which can be expressed in the following standard form: ρ (T ) c (T ) Tt = [λ (T ) Tx ]x subject to (controlled) Dirichlet boundary conditions

(1)

T (t, 0) = T (t, L) = u, (2) where ρ (T ) > 0 denotes the density, c (T ) > 0 indicates the specific heat capacity and λ (T ) > 0 describes the thermal conductivity of the good inside the thawer. Note that ρ(T ), c(T ) and λ (T ) all depend on the temperature T . For completeness we point out that u is the control variable (input), whereas its derivative u˙ = v is the optimization variable, see Section 3.2. The diffusion equation (1) can be rewritten as ρ (T ) c (T ) Tt = λT (T ) Tx2 + λ (T ) Txx .

(3)

due to λx (T ) = λT (T )Tx . To keep notation simple, two new parameters can be introduced as

IFAC DYCOPS-CAB, 2016 June 6-8, 2016. NTNU, Trondheim, Norway Christoph Josef Backi et al. / IFAC-PapersOnLine 49-7 (2016) 183–188

λ (T ) , ρ (T ) c (T ) λT (T ) κ (T ) = , ρ (T ) c (T ) leading to a rewritten form of (3):

−6

k (T ) =

2 −1

(4)

k(T) [m s ]

1

x 10

0.8 0.6 0.4 0.2

Tt =κ (T ) Tx2 + k (T ) Txx .

0 235

(5)

240

2 −1

2 −1

k(T) [m s ]

2 1

271 272 273 Temperature [K]

8.6 8.4 8.2 8 250

274

255 260 265 270 Temperature [K]

Figure 4 presents the parameter function κ(T ), which is only nonzero inside TF ±ΔT for the piecewise continuous definitions of λ (T ) and c(T ). For the arctan-approximation the values are nonzero also outside TF ± ΔT , however, very small. −9

0 −2 −1

x 10

−1 −2 −3 −4 240

245

250

255 260 265 Temperature [K]

270

275

280

285

271.5 272 272.5 Temperature [K]

273

−10

x 10

1

x 10

0

−4.6

−2 −1

k(T) [m s ]

−2 −1

k(T) [m s ]

−4.5

−4.7 −4.8 −4.9 −5 271

−1 −2 −3 −4

271.5 272 272.5 Temperature [K]

273

−5 271

Fig. 4. κ(T ), based on piecewise constant definitions of c(T ) and λ (T ) (blue) and approximated by arctan-expressions (red), the two bottom plots represents zooms into a specific temperature region

ΔT TF

285

TF = 272 K and ΔT = 0.5 K causes a significant decrease of the thermal diffusivity k(T ) in that region.

λl

cl

280

8.8

−9

ci

275

x 10

Fig. 3. k(T ), based on piecewise constant definitions of c(T ) and λ (T ) (blue) and approximated by arctan-expressions (red), the two bottom plots represents zooms into a specific temperature region

−5 235

T

270 −7

9

0 270

λs

ΔT

255 260 265 Temperature [K]

k(T) [m s ]

3

k(T) [m s ]

The model (5) does not explicitly model the physical behaviour at and around the freezing point TF . The presence of latent heat of fusion causes so called thermal arrest, which means that the temperature around the freezing point remains constant until the specific amount of latent heat of fusion is removed from (in the freezing case) or added to (in the thawing case) the object. One approach to model this phenomenon is to overestimate c(T ) in a region around the freezing point IΔT = TF ± ΔT , such that heat conduction is lowered in this region. This method is called apparent heat capacity method and has been described generally in e.g. Muhieddine et al. (2008). Figure 2 displays a qualitative sketch of the parameters λ (T ) and c(T ) over T . We remark that the temperature dependence of ρ(T ) is insignificant compared to λ (T ) and c(T ) and therefore approximated as ρ(T ) = 950 kg m−3 = const.

cs

250

x 10

2.3 Parameters

c

245

−7

For freezing, (5) has extensively been studied, for example in Backi et al. (2015), where a stability investigation as well as an observer design have been conducted, and in Backi et al. (2014), where in addition similarities to the (potential) Burgers’ equation were presented.

λ

185

T

Fig. 2. Qualitative sketch of parameters λ and c over T Parameter functions based on real experimental data defining c(T ) and λ (T ) as thermodynamic alloys of different substances such as water, fat, carbohydrates, etc. are described in Harðarson (1996) resulting in [λs λl ] = [1.8 0.5] W m−1 K−1 and [cs ci cl ] = [2200 283000 3800] J kg−1 K−1 . In Figures 3 and 4 the blue plots represent the parameter functions k(T ) and κ(T ) calculated with (4) and the definitions in Figure 2 whereas in red the approximations of k(T ) and κ(T ) based on arctan-expressions are shown. These arctan-functions are displayed in Appendix A and are used in the optimization software package ACADO (see Section 4). As can be seen from Figure 3, the overestimation of c(T ) in the region TF ±ΔT , with 185

3. MODEL PREDICTIVE CONTROL We are going to optimally control the temperature distribution inside the block of foodstuff during thawing by using an MPC algorithm. Therefore it will be necessary to discretize the PDE. Both, discretization in space and time are possible, however, as we are going to design a continuous MPC scheme, the discretization will only be performed for the spatial domain. 3.1 Discretization As outlined above, the model (5) is discretized in the spatial domain only, a so-called semi-discretization approach. This

IFAC DYCOPS-CAB, 2016 186 Christoph Josef Backi et al. / IFAC-PapersOnLine 49-7 (2016) 183–188 June 6-8, 2016. NTNU, Trondheim, Norway

leads to a set of continuous, coupled (interdependent) ODEs. Forward, backward and central difference schemes are used in order to discretize the term Tx , whereas only central differences are applied for the term Txx . The overall scheme holds for n = 1, 2, . . . , N, with N odd, and can be seen in the following ⎧ Tn+1 − Tn ⎪ ⎪ , if n ∈ RL = {1, 2, . . . , N−1 ⎪ 2 } ⎪ Δx ⎪ ⎨ Tn+1 − Tn−1 Tx = , if n = RM = N+1 2 ⎪ 2Δx ⎪ ⎪ (6) ⎪ T − T n n−1 ⎪ N+3 ⎩ , . . . , N}, , if n ∈ RR = { 2 , N+5 2 Δx Tn+1 − 2Tn + Tn−1 , Δx2 The fictional states (ghost cells) at T0 and TN+1 represent the boundary conditions (inputs; control variables) and enter the equations exclusively in the discretization of the second order partial derivative with respect to the state variable. Txx =

Discretizing (5) using (6) ultimately leads to a set of interdependent, nonlinear ODEs: ⎧ (Tn+1 − Tn )2 Tn+1 − 2Tn + Tn−1 ⎪ ⎪ + k(Tn ) , n ∈ RL ⎪κ(Tn ) ⎪ 2 ⎪ Δx Δx2 ⎨ 2 − T ) − 2T + T (T T n n−1 n+1 n−1 T˙n = κ(Tn ) n+1 + k(Tn ) , n = RM 2 2 ⎪ 4Δx Δx ⎪ ⎪ 2 ⎪ ⎪ ⎩κ(Tn ) (Tn−1 − Tn ) + k(Tn ) Tn+1 − 2Tn + Tn−1 , n ∈ RR , Δx2 Δx2 (7) It is clear that for large numbers of discretization steps the results will converge towards the “real” (analytical) solution of (5). In this work, however, the spatial discretization is chosen as N = 25 discretization steps. This is due to reduced computational complexity as the MPC algorithm is supposed to deliver its results in real-time. Nevertheless, we point to Backi (2015, Appendix A.2), where a comparative investigation between different discretization resolutions using the scheme in (6) and R the commercial MATLAB solver pdepe is presented. This investigation leads to the conclusion that a resolution of N = 25 covers the dynamics sufficiently well for the present application. 3.2 Algorithm The MPC consists of a cost function incorporating reference tracking for states, the rate of temperature change at the boundaries (optimization variable) and the temperature at the boundaries (control variable; input), respectively. In addition, we are going to introduce constraints to the system. We require the temperature at the boundaries to be in the band 283 K − 293 K, as temperatures below this band cause a too small driving force for heat exchange and temperatures above this band will lead to quality degradation (caused by protein denaturation). Furthermore, the temperature rate of change v is limited to 1 K per 100 s to account for the dynamics of the supplying components displayed in Figure 1. Hence, the cost functions looks as follows min J = v

�τ 0

(T − Tre f )T Q(T − Tre f ) + R(u − ure f )2 + S(v − vre f )2 dt

subject to 186

T˙ − f (T, u) = 0, u˙ = v T (0) = Tinit , u(0) = uinit , T ≤ T, v ≤ v ≤ v, u ≤ u ≤ u, where T˙ − f (T, u) = 0 denotes the set of interdependent ODEs (7) and u˙ = v indicates that the optimization variable v is the rate of change of the boundary temperature u. The remaining five equations indicate initial values for the temperatures (Tinit ) and the boundary temperatures (uinit ) as well as an upper bound for the temperatures (T ) and lower and upper bounds for the rate of change of the boundary temperature (v and v) and the boundary temperature itself (u and u). The positive definite matrix Q has only nonzero entries on the diagonal, whereas R and S are scalars. All numerical values are collected in the table in the Appendix. The nature of the thawer dictates that at first the outer layers of the block are heated up before the warm front propagates to the center of the spatial domain. Thus the reference temperature Tre f in the objective function gets chosen such that the reference values are increasing from the middle towards the outer layers. In the very center the reference temperature is chosen to be 272.5 K, which is just outside the latent zone, meaning that the block is fully thawed. The reference for the boundary temperature u is defined as 283 K, which is assumed to be the thawing medium’s initial temperature and the temperature of the cold water feed. This shall account for saving energy as less water will have to be drained from the hot water buffer tank. The reference for the rate of change of the boundary temperature, v, is set to 0 K s−1 as a steady boundary temperature is desired. 3.3 Stability of the MPC scheme Stability of the MPC algorithm is an important property and can be enforced for example by introducing so called terminal region constraints and/or terminal penalty terms (also for nonlinear MPC). This can be achieved by quasi-infinite horizon nonlinear MPC schemes, where the terminal region and the terminal penalty term are calculated off-line. The aim is to find an upper approximation for the infinite horizon cost functional, which itself guarantees stability, as described for example in Findeisen and Allgöwer (2002). These stability results hold for systems, whose steady state is the origin. The steady state for model (5), however, is dictated by the boundary conditions, as it is known that the steady state will converge to a straight line between the boundary conditions (Backi et al. (2015)). This means that with every change of the boundary temperature, different steady state values will be defined. To implement this into an MPC scheme seems impractical for the present case and therefore the system will not be transformed to have the origin as equilibrium point. However, it must be mentioned that there exist methods to handle stability enforcing MPC schemes for families of setpoints, as introduced in Findeisen and Allgöwer (2000). These methods are based upon pseudolinearization of the nonlinear system dynamics, which can be done off-line. This, however, has the disadvantage of being computationally excessive.

IFAC DYCOPS-CAB, 2016 June 6-8, 2016. NTNU, Trondheim, Norway Christoph Josef Backi et al. / IFAC-PapersOnLine 49-7 (2016) 183–188

The fact that the present case is a batch process, where a tracking rather than a classical stabilization problem has to be solved (the reference values explicitly differ from the steady state values), and that, in addition, the system dynamics are stable (bounded) for bounded inputs, makes a stability investigation unnecessary. 4. SOLVING STRATEGY AND SIMULATION RESULTS

187

and rendering the initially symmetric problem non-symmetric. The initial conditions are typically estimated; Observer designs related to these classes of problems are introduced in e.g. Backi et al. (2015) and Backi (2015). The simulation is terminated when all states are larger than 272.5 K, meaning that all states are above the latent zone IΔT . 4.2 Simulation results

Simulations have been conducted with the Optimization package ACADO for Matlab, see e.g. Houska et al. (2011) and Ariens et al. (2010). In order to use ACADO, the parameter functions k(T ) and κ(T ) have to be defined as continuous functions. The approximations based on arctan-expressions and the simulation parameters can be found in the Appendix.

Figure 5 shows the optimization variable (rate of change of the boundary temperature) together with its constraints in red, whereas the control variable (boundary temperature) itself and the respective constraints in red can be seen in Figure 6. It must be pointed out that actuation is only needed in roughly the first 4000 s with the chosen simulation parameters.

4.1 Solution strategy for the optimization problem 0.01

The time horizon is set to τ = 500 s, where 5 equidistant timesteps are chosen, and hence Δτ = 100 s. The first time sequence of 100 s is used as an input to the plant. This represents a trade-off between numerical accuracy and computational speed. After each MPC iteration the states obtained by simulation get perturbed by white Gaussian noise with power = 0.01 dBW before being fed back as new initial conditions to the optimizer. This shall model for both measurement and modeling errors, leading to an application, which is closer to the real world 187

Temperature change [Ks -1 ]

0.006 0.004 0.002 0 -0.002 -0.004 -0.006 -0.008 -0.01 0

0.5

1

1.5

2

2.5

3

Time [s]

×10

4

Fig. 5. Rate of change v of the boundary temperature over time

294

292

Temperature [K]

Discretizing the optimization and respectively the control variable (input) in the aforementioned fashion and insertion into the objective value and constraints defines a Nonlinear Programming (NLP) problem. An NLP can be solved by different methods (e.g. Interior point methods), however, Sequential Quadratic Programming (SQP), an iterative method, is often the most effective for shooting methods as they generate small and dense NLPs. This is done by defining quadratic (QP) optimization subproblems, which are solved by an active-set strategy relying on linearized dynamics and constraints. Here, line search algorithms are used in order to find a search direction leading to a decrease in the objective function. Thus, one SQP iteration delivers an optimal search direction, however, this may have been achieved by several iterations solving the QP. The integration method used to solve the QP is a Runge Kutta method of order 7/8 with relative tolerance 10−6 . The built-in QP solver in ACADO is qpOASES, an open-source C++ implementation of an active set strategy. The Hessian is chosen to be an approximation of Gauss-Newton type and the tolerance level for the Karush-Kuhn-Tucker conditions is set to 10−6 .

0.008

290

288

286

284

282 0

0.5

1

1.5

2

2.5

3

Time [s]

×10

4

Fig. 6. Boundary temperature u over time In Figure 7 specific temperatures over time are shown. Note that the temperatures for e.g. n = 1 and n = 25 are not identical due to the perturbation with white Gaussian noise. Block temperatures

285 280 275 270

Temperature [K]

The optimization problem is solved by a direct single shooting method. This represents a sequential approach, for which the time horizon τ gets divided into a fixed time grid 0 = t0 < t1 < · · · < tG = τ. This time grid does not necessarily have to be equidistant, however, due to convenience it is chosen as such. Nevertheless, it should be pointed out that some problems might be solved more efficiently and faster if the time grid is defined non-equidistant. Between each of these time instances a piecewise constant input function, the optimization variable, gets chosen and implemented into the ODEs. However, piecewise constant input functions are often not practical due to their stepwise definition. Therefore, as mentioned before, the rate of change of the boundary temperature v is chosen as the optimization variable, which, after integration, will result in a continuous boundary temperature u, which is inserted into (7).

265 n=1 n=3 n=5 n=7 n=9 n=11 n=13 n=15 n=17 n=19 n=21 n=23 n=25

260 255 250 245 240 235 230 0

0.5

1

1.5

2

Time [s]

Fig. 7. Specific temperatures over time

2.5

3 ×10

4

IFAC DYCOPS-CAB, 2016 188 Christoph Josef Backi et al. / IFAC-PapersOnLine 49-7 (2016) 183–188 June 6-8, 2016. NTNU, Trondheim, Norway

In Figure 8, temperatures at the end of the simulation time are visualized in blue together with the reference temperatures Tre f in red. 282 281

Temperature [K]

280 279 278 277 276 275 274 273 272 0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

Length [m]

Fig. 8. Temperatures at the end of simulation (blue) and Tre f (red) 5. CONCLUSION In this paper a method for thawing blocks of foodstuff in a contact thawing device was introduced. In addition, a functional principle for such a contact thawing device was proposed. Furthermore, a mathematical model for thawing processes based on the diffusion equation with state-dependent parameters was defined. In a simulation study, Model Predictive Control (MPC) was used in order to find an optimal solution for controlling the temperature field inside the block of foodstuff in the thawer. The simulation examples presented in Section 4 illustrate that the control objectives could be achieved. The constraints were satisfied even in the presence of added white Gaussian noise, which represents measurement noise and modeling errors. REFERENCES Ariens, D., Ferreau, H.J., Houska, B., and Logist, F. (2010). ACADO for Matlab User’s Manual. Optimization in Engineering Center (OPTEC) and Department of Electrical Engineering, K. U. Leuven, Leuven, Belgium, Version 1.0beta v2022 edition. Backi, C.J., Bendtsen, J.D., Leth, J., and Gravdahl, J.T. (2014). The nonlinear heat equation with state-dependent parameters and its connection to the Burgers’ and the potential Burgers’ equation. In Proceedings of the 19th IFAC World Congress. Cape Town, South Africa. Backi, C.J., Bendtsen, J.D., Leth, J., and Gravdahl, J.T. (2015). A heat equation for freezing processes with phase change: Stability analysis and applications. International Journal of Control, 89(4), 833–849. Backi, C.J. and Gravdahl, J.T. (2013). Optimal boundary control for the heat equation with application to freezing with phase change. In Proceedings of the 3rd Australian Control Conference. Perth, Australia. Backi, C.J. (2015). Modeling, Estimation and Control of Freezing and Thawing Processes - Theory and Applications. Ph.D. thesis, Norwegian University of Science and Technology, Trondheim, Norway. Cleland, D.J., Cleland, A.C., Earle, R.L., and Byrne, S.J. (1987). Experimental data for freezing and thawing of multidimensional objects. International Journal of Refrigeration, 10, 22–31. 188

Findeisen, R. and Allgöwer, F. (2000). A Nonlinear Model Predictive Control Scheme for the Stabilization of Setpoint Families. Journal A, Benelux Quarterly Journal on Automatic Control, 1(41), 37–45. Findeisen, R. and Allgöwer, F. (2002). An Introduction to Nonlinear Model Predictive Control. In Proceedings of the 21st Benelux Meeting on Systems and Control, 119–141. Veldhoven, The Netherlands. Harðarson, V. (1996). Matvarens termofysiske egenskaper og deres betydning ved dimensjonering av frysetunneler. Ph.D. thesis, Universitetet i Trondheim – Norges Tekniske Høgskole. Houska, B., Ferreau, H.J., and Diehl, M. (2011). ACADO Toolkit - an open-source framework for automatic control and dynamic optimization. Optimal Control Applications and Methods, 32, 298–312. Jason, A.C. (1974). Thawing Frozen Fish. Torry Advisory Note. Ministry of Agriculture, Fisheries and Food, Torry Research Station. Johnston, W.A., Nicholson, F.J., Roger, A., and Stroud, G.D. (1994). Freezing and refrigerated storage in fisheries. Technical Report 340, FAO Fisheries. Muhieddine, M., Canot, É., and March, R. (2008). Various approaches for solving problems in heat conduction with phase change. International Journal of Finite Volume Method, 6(1). Olver, P.J. (2014). Introduction to Partial Differential Equations. Undergraduate Texts in Mathematics. Springer. Pham, Q.T. (2006a). Mathematical modeling of freezing processes. In Handbook of Frozen Food Processing and Packaging, chapter 7. Taylor & Francis Group, LLC. Pham, Q.T. (2006b). Modelling heat and mass transfer in frozen foods: a review. International Journal of Refrigeration, 29, 876–888. Appendix A. SIMULATION PARAMETERS AND APPROXIMATIONS FOR k(T ) AND κ(T ) Tre f ure f vre f Tinit uinit u u v v T Q

[278.5 278 277.5 277 276.5 276 275.5 275 274.5 274 273.5 273 272.5 273 273.5 274 274.5 275 275.5 276 276.5 277 277.5 278 278.5]T K u 0 Ks-1 235 · 1(25×1) u 283 K 293 K −0.01 K s-1 0.01 K s-1 285 · 1(25×1) I(25×25)

R S

10−2 10−2

1(n×1) denotes a vector of ones with n elements.      T π k (T ) = 7 · 10−8 + 2.737 · 10−7 − arctan − 1 20000π + TF − ΔT 2     T −8 + 4.35 · 10 − 1 50000π arctan TF + ΔT     T κ (T ) = 1.543 · 10−9 − arctan − 1 50000π TF − ΔT     T −9 + 1.543 · 10 − 1 50000π arctan TF + ΔT