CET vol 99 DOI: 10.3303/CET2399111 Paper Received: 14 December 2022; Revised: 30 March 2023; Accepted: 22 April 2023 Please cite this article as: Balogh L., Egedy A., Barkanyi A., 2023, Model Predictive Control of CMSMPR Crystalliser, Chemical Engineering Transactions, 99, 661-666 DOI:10.3303/CET2399111 CHEMICAL ENGINEERING TRANSACTIONS VOL. 99, 2023 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Sauro Pierucci, Flavio Manenti Copyright © 2023, AIDIC Servizi S.r.l. ISBN 978-88-95608-98-3; ISSN 2283-9216 Model Predictive Control of CMSMPR Crystalliser László Balogh*, Attila Egedy, Ágnes Bárkányi University of Pannonia, Faculty of Engineering, Department of Process Engineering, Egyetem street 10, Veszprém, 8200 Hungary balogh.laszlo@mk.uni-pannon.hu One of the most critical components of the chemical industry in terms of crystallisation is the pharmaceutical sector. Most medicine components are expensive and require complex processes for their production, so producing waste is highly inefficient. Another concern is the high-quality standards for most pharmaceutical products. Therefore, optimising the crystallisation process is critical from a quality perspective, with the main concerns being the product's crystal structure and particle diameter distribution. Regardless efficient control in batch processes such as crystallisation is a difficult task due to the inherently nonlinear behaviour of the system. Using a priori model of the system as the basis for nonlinear model predictive control could provide a useful tool for handling the crystallisation process, mitigating the effects of disturbance and noise and ensuring appropriate product quality. In this work, we wish to showcase the possibility of controlling a crystallisation process using model predictive control to enable the production of crystal products with desired particle diameter distribution and crystalline product average size. The method is shown using citric acid as a model substance in a case study of a continuous crystallisation procedure in a stirred tank reactor. The crystalliser model includes an energy balance, so the system's behaviour depends on the cooling rate and residence time. Accordingly, the control problem can be formulated as multiple inputs and multiple outputs (MIMO) system. Moreover, the two controlled (average particle size and crystal size dispersion) variables are not easily detached from each other. So, the traditional controlling strategies, for example, the decoupling controller, is challenging to apply. The MPC (model predictive control) as an advanced control algorithm can be a solution to this. 1. Introduction Control of crystallisers is challenging because the process is influenced by many variables, for example, the residence time, temperature, supersaturation, or cooling rate. However, using an MPC, the cross effects between the variables and the controlled properties become manageable. The classical control strategies base problem, as the intervention effect will apply later, can be eliminated. This is especially true for the crystallization processes, where the manipulated variables (residence time, cooling medium temperature and flow rate has a complex relationship with the controlled valriables: particle size, and distribution etc. Yang and Nagy (2015) examined a two-stage cascade, continuous MSMPR (mixed suspension mixed product removal) crystalliser and the application of NMPC (nonlinear MPC) for controlling the average particle size and the yield of the crystalliser system. In their case, the operational variables are the antisolvent flow rates and the temperatures in the two equipment. Szilágyi et al. (2018) investigate an experimental implementation of PBM (population balance model) based nonlinear MPC to l-ascorbic acid crystallisation. The aim of that study is quality-based control (QbC). Shi et al. (2005) were engaged with the particle sizes distributions predictive control. In their implementation, the MPC was supplemented with Lyapunov-based bounded control; thus, a supervisor algorithm keeps the system in a stable region. Moldoványi et al. (2004) studied an MPC of continuous MSMPR crystallisers. They examined more operation variables based on the system's RGA (relative gain array), for example, inlet concentration, flow rate, and seeding. In this study, the base of the MPC is the linearised steady-state model of the crystalliser. Based on the literary review it can be concluded that while many new advancements have been made on the use of MPC the utilization of robust nonlinear first principle models as the basis for optimized control is still a relevant topic with many opportunities for innovation and study. 661 mailto:balogh.laszlo@mk.uni-pannon.hu In our study, we first explain the mathematical model of the CMSMPR (continuous MSMPR) crystalliser and the assumptions we use. An essential difference between our work and previous studies, we implemented an energy balance to the crystalliser model, and one of the manipulating variables acts through it. Thus, the control dynamics and the possibility of intervention appear more realistic. After the presentation of the model, we examined the system’s (citric acid crystalliser) dynamic behaviour with repeating sequence stair test functions. Then, we present the result of the closed-loop MPC. 2. Mathematical model for CMSMPR crystalliser The particle size distribution change can be described with a population balance equation similar to the general balance equations. The local change depends on the change along coordinates and the source terms. (Ramkrishna, 2000). That can be formalised the follows. ∂𝛹 𝜕𝑡 + ∇𝑟(𝛹 �̇�) + 𝐺 = 0 (1) It is assumed that the following conditions are satisfied: • The working volume of the crystallisation is constant; we assume ideal level control. • The inlet does not contain crystal seeds or any impurities. • The crystal growth (𝑣𝑟) is independent of the crystal size, so the growth rate (𝑘𝑔) is linear (𝛼 = 0) and can be formalised with a phenomenology power law equation. The 𝑘𝑔 and 𝑔 are greater than zero and constants in Eq (3). The supersaturation formalised the following where the 𝑐𝑓 is the concentration of the inside of the crystalliser and the 𝑐𝑓 ∗ is the solubility (Eq 2). Δ𝑐 = 𝑐𝑓 − 𝑐𝑓 ∗(𝑇) (2) The crystal growth can be formalised in the following (Eq 3). 𝑣𝑟 = 𝑘𝑔𝛥𝑐𝑔(1 + 𝛼𝑟 ) (3) • Nucleation (𝐵) consists of primary (𝐵𝑝) and secondary (𝐵𝑠) nucleation and both phenomena can be formalised with a phenomenology power law equation. In the case of primary nucleation, the nucleation rate depends only on supersaturation (Δ𝑐). In contrast, in the case of secondary nucleation, the volume of the solid crystals (𝑉𝑠) influences the rate of nucleation. The 𝑘𝑝, 𝑝, 𝑘𝑠 and 𝑠 are greater than zero and constants. The crystal grown can be formalised in the following (Eq 4). 𝐵 = 𝐵𝑝 + 𝐵𝑠 = 𝑘𝑝 Δ𝑐𝑝+ 𝑘𝑠Δ𝑐𝑠𝑉𝑆 𝜔 (4) • The suspension inside the crystalliser is perfectly mixed. Therefore, the outlet stream has the same composition as the one inside the equipment. • For easy manageability, we introduce the crystal volume fraction, which denotes the epsilon (Eq 5). 𝜖 = 𝑉𝑠 𝑉𝑠 + 𝑉𝑓 (5) Considering these assumptions, Eq (1) takes the following form (Eq 6) 𝜕𝛹(𝑟, 𝑡) 𝜕𝑡 + 𝑣𝑟 𝜕𝛹(𝑟, 𝑡) 𝜕𝑟 = 𝑞 𝛹(𝑟, 𝑡) 𝑉 + 𝐵 [1 − 𝜖(𝑡)]𝛿(𝑟 − 𝐿0) (6) We want to design a model predictive control for the equipment, so we need a relatively fast solver method; therefore, we chose the moments method. To describe the present system, we need the first 3 moments. The definition of the mth moments can be formalised the follows (Eq 7). 𝜇𝑚 = ∫ 𝑟𝑚𝛹(𝑟, 𝑡) 𝑑𝑟 ∞ 0 where 𝑚 = 0, 1, … , 3 (7) After the moments transformation, we get an ordinary differential equation system, the general form of which is as follows (Eq. 8). 𝑑𝜇𝑚 𝑑𝑡 = 𝑚 𝑣𝑟 𝜇𝑚−1 − 𝑞 𝑉 𝜇𝑚 + 𝐵 𝜖 𝐿0 𝑚 where 𝑚 = 0, 1, … , 3 (8) The size of 𝐿0 is difficult to measure, and a small quantity and the powers of that are even smaller in the case of 𝑚 > 0 is negligible. The change in the macroscopic properties (𝑐𝑓 , 𝑐𝑤 , 𝑉𝑓 , 𝑉𝑠, 𝑇𝑟 , 𝑇𝑗) of the crystalliser can be formalised in the following equation system (Eq 9, 10, 11, 12, 13, 14). 662 𝑑(𝑐𝑓𝑉𝑓) 𝑑𝑡 = 𝑞𝑐𝑓,𝑖𝑛 − (1 − 𝜖)𝑞𝑐𝑓 − 𝑣𝑟𝜙3𝜇2𝑉𝜌𝑠 (9) 𝑑(𝑐𝑤𝑉𝑓) 𝑑𝑡 = 𝑞𝑐𝑤,𝑖𝑛 − (1 − 𝜖)𝑞𝑐𝑤 (10) 𝑑𝑉𝑓 𝑑𝑡 = 𝑞 − (1 − 𝜖)𝑞 − 𝑣𝑟𝜙3𝜇2𝑉 (11) 𝑑𝑉𝑠 𝑑𝑡 = −𝜖𝑞 + 𝑣𝑟𝜙3𝜇2𝑉 (12) 𝑑𝑇𝑟(𝑉𝑠𝜌𝑠 + 𝑉𝑓𝑐𝑓𝑐𝑝𝑐 + 𝑉𝑓𝑐𝑤𝑐𝑝,𝑤) 𝑑𝑡 = 𝑞𝑇𝑟,𝑖𝑛(𝑐𝑤,𝑖𝑛𝑐𝑝,𝑤 + 𝑐𝑓,𝑖𝑛𝑐𝑝𝑐) − 𝑞(1 − 𝜖)𝑇𝑟(𝑐𝑤𝑐𝑝,𝑤 + 𝑐𝑓𝑐𝑝𝑐) −𝑞𝜖𝑇𝑟𝜌𝑠𝑐𝑝,𝑠 − 𝑈𝐹(𝑇𝑟 − 𝑇𝑗) + Δ𝐻𝑘𝑟𝑖𝑠𝑡𝑅𝑉𝑠 𝑉 (13) 𝑑𝑇𝑗𝑉𝑗𝜌𝑗𝑐𝑝,𝑗 𝑑𝑡 = 𝑞𝑗𝜌𝑗𝑐𝑝,𝑗𝑇𝑗,𝑖𝑛 − 𝑞𝜌𝑗𝑐𝑝,𝑗𝑇𝑗 + 𝑈𝐹(𝑇𝑟 − 𝑇𝑗) (14) The temperature dependence of thermodynamic properties and solubility was calculated based on Apelblat 2014, with fitted polynoms. The heat of crystallisation was neglected. The operating parameters of the equipment were chosen arbitrarily, except for the 𝑐𝑓,𝑖𝑛 because this was calculated from the solubility curve. The parameters of the model can be found in Table 1. Table 1: Kinetic and operating parameters of the case study of the crystalliser model Parameter Value Unit Reference Parameter Value Unit 𝑘𝑝 2.869 ⋅ 103 [ # 𝑚3𝑠 ] Nemdili et al. (2016) 𝑉 1.0 [𝑚3] 𝑝 1.585 [−] 𝑉𝑗 0.235 [𝑚3] 𝑘𝑠 1.72 ⋅ 108 [ # 𝑚3𝑠 ] Févotte and Févotte (2009) 𝐹 4.2 [𝑚2] 𝑠 0.47 [−] 𝑈 1500 [ 𝑊 𝑚2𝐾 ] 𝜔 1.14 [−] 𝑇𝑡,𝑖𝑛 50 [°𝐶] 𝑘𝑔 7.18 ⋅ 10−6 [ 𝑚 𝑠 ] Févotte and Févotte (2009) 𝑇𝑗,𝑖𝑛 -10 [°𝐶] 𝑔 1.58 [−] 𝑐𝑓,𝑖𝑛 774.066 [ 𝑘𝑔 𝑚3] 𝜙 𝜋 6 [−] Caillet et al. 2007 3. Model predictive control design In the MPC structure, the control algorithm is a model of the controlled object supplemented with an optimiser (Szilágyi et al. 2018). The MPC can be formulated as a multivariable, conditional optimisation problem. In our case, the controller model and the controlled object are the same a priori model. Another approximation compared to the real case is that the controller model gets all variable from the controlled model and use them as an initial condition. The structure of MPC can be seen in Figure 1. Figure 1: Structure of MPC For the control, the flow rate to the inside of the crystalliser (𝑞) and to the jacket (𝑞𝑗) was chosen as an manipulated variable, and the controlled variables are the average particle size (𝐸), the variance of the crystals 663 (𝜎) and the volume of the crystalline product (𝑉𝑠) because these properties are the most important industrially (Moldoványi 2005). So, the control (𝑦) and the manipulated variables (𝑢) are formalised as follows (Eq 15). 𝑦1 = 𝐸 = 𝜇1 𝜇0 𝑦2 = 𝜎 = 𝜇2 𝜇0 − ( 𝜇1 𝜇0 ) 2 𝑦3 = 𝑉𝑠 𝑢1,2 = [𝑞; 𝑞𝑗] (15) The optimisation problem is formalised as the following (Eq 16): min 𝑢𝑖(𝑡𝑀,0),𝑢𝑖(𝑡𝑀,0+𝜒),… ,𝑢𝑖(𝑡𝑀,0+𝑡𝑀) ∑ 𝜁𝑗 ∫ (𝑦𝑗 − 𝑦𝑗 ref) 2 𝑡𝑃 𝑡𝑃,0 𝑑𝑡, 𝑛 𝑗 𝑖 = 1,2 𝑗 = 1,2,3 (16) The model length of the horizon is 𝑡𝑀,0 to 𝑡𝑀 and the length of time steps is 𝜒, while the prediction horizon is 𝑡𝑃,0 to 𝑡𝑃, and the 𝜁 is a weight factor. The optimisation problem is solved online parallel with the model at certain intervals, and the length of this interval (𝜏) is also a parameter in the MPC algorithm. The solution is made in MATLAB/Simulink environment with fmincon function. The fmincon function was chosen because the method can be parallelised, fast, and easily parametrised. The constraints of the minimum search for both variables is the same, the minimum is 0.001 and the maximum is 0.1 m3/s. Table 2 shows the parameters of the optimisation problem. Table 2: Parameters of the MPC Parameter Value Parameter Value Parameter Value 𝑡𝑀,0 0 s 𝑡𝑃,0 500 s 𝜏 100 s 𝑡𝑀 3000 s 𝑡𝑃 3000 s 𝜒 0.001 s 𝜁1 0.5 𝜁2 0 𝜁3 0.5 The control of average crystal size and the size dispersion simultaneously by the supersaturation and the residence time would be challenging (Moldoványi 2005). So, in this study, we consider only the crystalline product's average crystal size and volume. Another manipulating variable should be examined in further work, for example, seeding. Along the tests, at first, the dynamic behaviour of the system was examined 3.1 System dynamic analysis For the dynamic analysis, the response of the system was examined with different flowrates (𝑞) with permanent coolant flow rates (𝑞𝑗), and inversely. The results can be seen in Figure 2. a) b) Figure 2: Dynamic behaviour of the crystallisation In Figure 2a can be seen that if the flow rate is decreasing, the examined parameters are increased, but the rate of change is different in the case of solid volume and expected value of crystal size. This can be explained by the resident time of the crystals being higher if the flow rate is lower because the crystals have more time to grow. Figure 2b shows the opposite trends for the intervention; this is caused by the higher coolant flow rate cooling it down more the inside space and pushing away the supersaturation. In comparing the two cases, the system's dynamic behaviour is faster for volume flow change because the inlet flow rate directly affects the 664 change of moments, while the coolant flow rate is indirectly through the energy balance. The behaviour of the model is as expected, which shows the same trends as in Yang and Nagy (2015). 3.2 Nonlinear model predictive control evaluation For the tests, we examined objects with closed-loop MPC. The set points change every 1000 s (to allow the system to reach a stationary state), and the system response was investigated. In Figure 3. we indicated the set points, the predicted output (at the end of the prediction horizon) and the response of the system, as well as the two manipulated variable values. Figure 3 shows that the two controlled variables follow the set points; the predicted output and the system response coincide with a good approximation. The first time, the system starts from an initial state, which causes the deviation from the set point, but the predicted output in the present points also shows the same values as the set point. The swings at 3000 to 4000 seconds (and in approximately 6000 and 7000 s) are likely the extreme value search algorithm’s mistakes. Figure 4 shows the temperature changes during the tests, which also follow the flow rate change. Figure 3: Operation of the MPC Figure 4: The temperature changes during the tests. 665 4. Conclusions The results show that the crystalline product volume and the expected volume of the crystal are independently adjustable quantities. However, it should be noted that only in a narrow environment. In our method, the energy balance was implemented, and the one manipulating variables acts through the energy balance, and the algorithm can adjust the controlled system to the setpoint. During further investigations, we will develop the appropriate manipulating variables based on the system's RGA (relative gain array), and we will develop state estimation procedures; thus, based on the controlled object measurable variables, we can operate the MPC. Moreover finally, we will examine the effects of tuning parameters. Nomenclature t – time, s tM0, tM – times of model horizon, s tP0, tP – times of prediction horizon, s r – crystal size, m vr – grow rate, m/s �̇� – Vector of the rate of change along internal crystal properties G – general source term RVs – volume source term, m3/s cf – concentration of fluid phase, kg/kg solvent cf,in – inlet concentration of fluid phase, kg/kg solvent cf * – solubility concentration, kg/kg solvent Δc – driving force, kg/kg solvent cw – solvent concentration cw,in – inlet solvent concentration cp,j, cp,w, cp,s – specific heat capacity, J/kgK ρj, ρw, ρs – density, kg/m3 Bp – primary nucleation, #/m3s Bs – secondary nucleation, #/m3s kg, g, α – parameters of growing function kp, p – parameters of primary nucleation ks, s, ω – parameters of secondary nucleation 𝜖 – crystal fraction factor Vf – volume of fluid phase, m3 Vs – volume of solid phase, m3 V – volume of crystalliser, m3 F – heat transfer surface, m2 U – Overall heat transfer coefficient, W/m2K m – order of the moments q – inlet volume flow rate, m3/s qj – coolant volume flow rate, m3/s Tr – temperature of inside, °C Tj – temperature of jacket, °C ΔHr – heat of crystallization u – manipulating variable y – controlled variables E – expected value of crystals size, m 𝜎 – variance of crystals size, m2 L0 – smallest crystal size, m 𝜁 – weight factor 𝜒 – simulation step time, s 𝜏 – MPC call time frequency, s Φ – crystal shape factor μm – moments, #/m4-m-1 δ – Dirac delta function Ψ – size distribution density function, #/m4 References Apelblat A, 2014, Properties of Citric Acid and Its Solutions, Chapter in: Citric Acid, Springer International Publishing, Switzerland, 13–130. Caillet A., Sheibat-Othman N., Fevotte G., 2007, Crystallization of Monohydrate Citric Acid. 2. Modeling through Population Balance Equations, Crystal Growth & Design, 7, 10, 2088-2095 Févotte F., Févotte G., 2009, A new approach for the modelling of crystallisation processes in impure media using Population Balance Equations (PBE), IFAC Proceedings Volumes, 42, 11, 52-61. Moldoványi N., Lakatos B. G., Szeifert F., 2014, Model predictive control of MSMPR crystallisers, Journal of Crystal Growth, 275, 1-2, 1349-1354. Nemdili L., Koutchoukali O., Bouhelassa M., Seidel J., Mameri F., Ulrich J., 2016, Crystallisation kinetics of citric acid anhydrate, Journal of Crystal Growth, 451, 88-94. Ramkrishna D., 2000, The Framework of Population Balance, Chapter in: Population Balance - Theory and Applications to Particulate Systems in Engineering, Academic Press: Cambridge, MA, USA Shi D., El-Farra N. H., Li M., Mhaskar P., Christofides P. D., 2015, Predictive control of particle size distribution in particulate processes, Chemical Engineering Science, 61, 1, 268-281. Szilágyi B., Borsos Á., Pal K., Nagy Z. K., 2018, Experimental implementation of a Quality-by-Control (QbC) framework using a mechanistic PBM-based nonlinear model predictive control involving chord length distribution measurement for the batch cooling crystallisation of l-ascorbic acid, Chemical Engineering Science, 195, 335-346. Yang Y. and Nagy Z. K., 2015, Application of nonlinear model predictive control in continuous crystallisation systems, American Control Conference (ACC), 4282-4287. 666 393balogh.pdf Model Predictive Control of CMSMPR Crystalliser