Applying machine learning (ML) techniques is a complex task when the data quality is poor . Integrating first-principle models and ML techniques, namely hybrid modelling significantly supports this task . This paper introduces a novel approach to develop ing a hybrid model for dynamic chemical systems . The case in analysis employs one first -principle structure and t wo ML-based predictors . Two training approaches ( serial and parallel), two optimi sers ( particle swarm optimi sation and differential evolution) and two ML functions ( multivariate rational function and polynomial) are tested. The polynomial function trained with the differential evolution showed the most accurate and robust results . The training approach does not significantly affect the hybrid model accuracy . However , the main effect of the training approach is on the robustnes s of the parameter predictions . The coefficients of determination (R2) on the test batch es are above 0.95. In addition, it showed satisfactory extrapolation capabilitie s on different production scales with R2>0.9. Keywords : hybrid modelling, dynamic system, solvent switch, optimization, statistical modelling, data value 1. Introduction In recent year s, the process industry is gaining several benefits in the application of big data and machine learning (ML) techniques to model their processes (Bogojeski et al., 2021; Chiang et al., 2017; Ge et al., 2021; Mohd Ali et al., 2015; Qin and Chiang, 2019) . These techniques use experimental data to develop data -driven models (or statistical models ) to predict the physical phenomena described in the dataset used for the training. Despite the augmented prediction accuracy and the interest from the chemical in dustry, the application of ML modelling is still limited to specific industrial -scale processes because of a data quality problem. The data should respect the 5Vs rule to apply ML techniques (BBVA, 2020; Demchenko et al., 2014) , namely: 1. Volume – higher amount of data allows for higher model accuracy . 2. Variety – an observed phenomenon should be describ ed using different kind s of data. 3. Velocity - the data should be available in real -time to refine the parameters of the model . 4. Veracity – the uncertainty related to the data should be available and low . 5. Value – the information contained in the data should be help ful and not redundant. The data obtained from chemical processes lacks Value and often Veracity . The p rocess industry can benefit from the high amount of data coming from production processes. However, they usually have little use in building prediction models since 1) they cover only a narrow parameter space and 2) the 2 isolated effect of each proce ss parameter on the final output is often not captured by the data (Reis et al., 2016) . In addition, the information obtained from the plant is highly intercorrelated (McBride et al., 2020) and often suffers missed information because of communication and storage problems (Fisher et al., 2020) . The hybrid modelling approach is a solution to apply statistical modelling to cases with low availability and quality of data . It employs first-principle model ling and statistical model ling to make its predictions . The statistical model can be used either to refine the prediction of the first -principle model or to predict the parameter s utilised by it to perform the predictions. It increase s the overall accuracy and extrapolation capabilities of the model (McBride et al., 2020; Sansana et al., 2021; von Stosch et al., 2014) . For this reason , the hybrid modelling overcomes the abovementioned issu es related to the narrow investigation space that characteri ses data coming from a chemical process . It is well known that the statistical models can perform their predictions only within the boundaries with in which the data are located. On the contrary, the hybrid models can extend their predictions range beyond the data boundaries since they can rely on the presence of a first-principle model . Moreover , the integration of the first-principle model allows th e incorporation of the first -principle correlation between the variables before the model training (Glassey and von Stosch, 2020; von Stosch et al., 2014) . This reduces the amount of data required to train the model. For these reasons, we hypothesi se that using a hybrid model is a solution to Value and Veracity insufficiencies in chemical p rocess data. This work utili ses the constant volume solvent switch (CSS) as a case study . It is one of the primary downstream operation s required to produce active pharmaceutical ingredient s (API). I t aims to switch the reaction solvent to the crystalli sation solvent allow ing the purification of the API (Elgue et al., 2006; Papadakis et al., 2016) . Developing a model with high prediction capabilities of the CSS operation by using the first-principle approach is challenging to execute . The system has a significant dynamic behavio ur, making it too complex to accurately model using thermodynamic equations . In addition, the interaction of the API with the various solvents utili sed in this pr ocess is hard to predict with high accuracy. As a matter of fact, the first-principle model used in this work return s poor prediction accuracy ( Table 4). This paper proposes a new methodology to develop hybrid models in the dynamic system running under uncertaint ies, focusing on the separation systems. We intend to overcome the Value and Vera city issues mentioned above . The hybrid model utilised in this paper is composed of one first - principle model and two statistical models. Eight combinations of statistical modelling techniques have been investigated in this work . These are the combinations of two statistical functions, two optimi sers, and two training strategies. The statistical functions are namely multivariate rational function (MRF) and polynomial . The optimi sers are particle swarm optimi sation (PSO) and differential evolution (DE). The training of the models w as performed using parallel and serial strateg ies. The usage of MRF has been briefly explored in chemical engineering for the approximation of functions not described by first -principle models (Willis and von Stosch, 2017) , while it has show n excellen t capabilities in other fields (Austin et al., 2021; Hong et al., 2018) . The DE and PSO algorithms have already been applied to chemical engineering problems to perf orm both parameter identification and process optimi sation (Aguitoni et al., 2018; Angira and Santosh, 2007; Biazi et al., 2020; Dragoi and Curteanu, 2016; Mariano et al., 2011; Schwaab et al., 2008) . They are valid alternatives to the traditional gradient -based optimi sation approaches due to their ability to escape from the local minimum and converge to the global optimum in the search area (Dragoi and Curteanu, 2016) . We hypothesi se that combining these powerf ul techniques leads to a hybrid model with high prediction accuracy and robustness. This way, we resolve the Value and Veracity issues with high extrapolation capabilities. The new methodology relies on the capability of MRF and polynomial to predict and h andle non -linear behavio urs and the augmented quality of the parameter tuning of the 3 differential evolution approach compared to the traditional gradient -based approach es. This new hybrid model methodology focuse s on batch and fed -batch operations. However , it applies to all the separations and reaction processes. The model developed with this approach is validated t hrough data from an industrial plant provided by a p harmaceutical partner company. 2. Material and methods 2.1 The process The data used in this study were obtained from a process in development. The process aim s to swap the solvent from dichloromethane (DCM) to acetone (Ac). A low concentration (below 10% mol) of methanol (MeOH) is present within the reboiler at the beginning of the operation. Figure 1 shows the schematic representation of the equipment used for this operation. The CSS operation is divided in to two phases according to the optimi sation strategy suggested in the literature (Elgue et al., 2006) . These are , namely, the first evaporation and constant volume phases. During the first evaporation step , the reaction solvent (DCM in this study) is evaporated to a minimum volume optimi sed during the design phase. In the constant volume phase, the crystalli sation solvent (Ac) is fed while the evaporation of the solvent blend (DCM +Ac) takes place at the designed minimum volume. Figure 1: Schematic representation of the solvent switch process utili sed as a case study in this work. At the initial state , the reboiler contains only dichloromethane (DCM). The acetone (Ac) is contained in a tank and fed to the reboiler. During the solvent switch operation , Ac and DCM evaporate from the reboiler while the Ac is fed to keep the liquid volume within the reboiler constant . This incre ases the concentration of Ac in the solution within the reboiler and decrease s the DCM concentration . 2.2 Data description and input -output of the final model Three batches are utili sed in this study . All the batches were executed in the same reactor, with the same equipment and following the same control strategy. In addition, a supplementary batch was utilised to investigate the extrapolation capabilities of the model. This batch data was used following the training phase to test the extrapol ation abilities of the model. This extrapolation batch was executed with half the volume compared to the training batches. A list of all the process variable s contained in the datasets and utilised to develop the model is available in Supplementary Informa tion (SI) section 1.2 . The model utili ses as input the flowrate of Ac, the temperature of the heat transfer fluid within the jacke t, the reflux direction and the stirring speed. It predicts all the other variables contained in the dataset; however, the mai n ones are the temperature inside the reboiler, the volume inside the reboiler and the evaporated amount of solvent leaving the reboiler . 2.3 The first-principle model 4 The first-principle model contains material balance, phase equilibrium equations and heat balance. It was built with the following assumptions: 1) Thermodynamic equilibrium between liquid and vapo ur, 2) The pressure of the system is constant during the operation, 3) The viscous dissipation heating due to the st irring is negli gible, 4) The fluid within the reboiler is perfectly mixed, 5) The condenser is ideal (all the vapo ur entering into it is condensate without any vapo ur-liquid equilibria) , 6) No chemical reaction s take place during the process, 7) Inert gas is not dissolved within the liquid, 8) The gas phase behaves as ideal. A detailed description of the first-principle model is available in SI section 1.3. The energy balance was written by considering the evapor ation of the solvent blend driven by the service fluid. The overall heat -balance equation is reported in (1): 𝑑(𝑀𝐿∗ℎ𝐿+𝑀𝑉∗ℎ𝑉) 𝑑𝑡= =𝐹𝑖𝑛∗ℎ𝐹+𝑅∗𝐹𝑟𝑒𝑓∗ℎ𝑟𝑒𝑓−𝐹𝑜𝑢𝑡∗ℎ𝑉+𝐹𝑔𝑎𝑠𝑖𝑛𝑒𝑟𝑡∗ℎ𝑔𝑎𝑠+𝑈∗𝐴∗(𝑇𝑠𝑒𝑟𝑣𝑖𝑐𝑒 −𝑇) (1) where 𝑀𝐿 [𝑚𝑜𝑙 ] is the molar amount of liquid within the reboiler, ℎ𝐿[𝐽∗𝑚𝑜𝑙−1] is the liquid enthalpy, 𝑀𝑉 [𝑚𝑜𝑙 ] is the molar amount of vapo ur within the reboiler, ℎ𝑉[𝐽∗𝑚𝑜𝑙−1] is the liquid enthalpy, 𝐹𝑖𝑛[𝑚𝑜𝑙 ∗𝑠−1] is the molar flowrate of crystall isation solvent flowing to the reboiler, ℎ𝐹 [𝐽∗ 𝑚𝑜𝑙−1] is the enthalpy of the crystalli sation solvent flowing to the reboiler, 𝑅[−] is the state of the reflux (0 or 1), 𝐹𝑟𝑒𝑓[𝑚𝑜𝑙 ∗𝑠−1] is the flowrate of the reflux stream, ℎ𝑟𝑒𝑓[𝐽∗𝑚𝑜𝑙−1] is the enthalpy of the reflux liquid stream, 𝐹𝑜𝑢𝑡 [𝑚𝑜𝑙 ∗𝑠−1] is the flowrate of vapo ur leaving the reboiler, 𝐹𝑔𝑎𝑠𝑖𝑛𝑒𝑟𝑡[𝑚𝑜𝑙 ∗ 𝑠−1] is the flowrate of inter gas flowing into the system, ℎ𝑔𝑎𝑠[𝐽∗𝑚𝑜𝑙−1] is the enthalpy of the inert gas, 𝑈 [𝐽∗𝑠−1∗𝐾−1∗𝑚−2] is the overall heat transfer coefficient between the liquid within the reboiler and the heat transfer fluid, 𝐴 [𝑚2] is the contact area between the liquid within the reboiler and the heat transfer fluid, 𝑇𝑠𝑒𝑟𝑣𝑖𝑐𝑒 [𝐾] is the temperature of the heat transfer fluid, 𝑇 [𝐾] is the temperature of the fluids (liquid and vapo ur) within the reboiler. The heat -transfer coefficient from the jacket to the reboiler (𝑈) was modelled utili sing the he at resistance chain model (Bird et al., 2006) . For the study case , the wall is composed of two solid layers (steel/glass enamel) and the convection layer inside the reboiler; the resistance related to the convection layer jacket -side was considered negligible. The convection heat transfer coefficient inside the reboiler was calculated by using the Nusselt equation correlation (2). 𝑁𝑢=ℎ∗𝐷 𝑘=𝑓(𝑅𝑒,𝑃𝑟,𝑉𝑖)=𝑘𝑈𝐴∗𝑅𝑒𝛼∗𝑃𝑟𝛽∗𝑉𝑖𝛾 (2) where 𝑁𝑢 [−] is the Nusselt number, ℎ [𝐽∗𝑠−1∗𝐾−1∗𝑚−2] is the thermal convection coefficient inside the reboiler, 𝐷 [𝑚] is the characteristic dimension of the system, 𝑘 [𝑊∗𝑚−1∗𝐾−1] is the thermal conductivity of the liquid inside the reboiler, 𝑘𝑈𝐴[−], 𝛼[−], 𝛽[−] and 𝛾[−] are case - dependent parameters, 𝑅𝑒[−] is the Reynolds number, 𝑃𝑟[−] is the Prandtl number , and 𝑉𝑖[−] is the ratio between the viscosity of the liquid on the wall of the reboiler and in the liquid bulk. To compute (2) four case -dependent parameters should be identified. In the literature , the following parameter values are suggested 𝑘∈[0.3,1.5],𝛼=0.66,𝛽=0.33,𝛾=0.14 (Mohan et al., 1992) . 2.4 The structures of the statistical models Applying a first -principle model with the parameter optimi sed for the analysis case returned poor prediction accuracy and low generali sation capabilities among the batches. Two main deviating parameters were identified from the sensitivity analysis on the parameters of the first -principle model . Namely, the heat transfer coefficient (3) and the liquid enthalpy (4). Two statistical functions were employed t o correct the identified parameters and model their deviat ions: 𝑘𝑈𝐴𝐻𝑀=𝑘𝑈𝐴𝐼𝐷𝐸𝐴𝐿+𝑘𝑈𝐴𝑑𝑒𝑣(3) ℎ𝐿𝐻𝑀=ℎ𝐿𝐼𝐷𝐸𝐴𝐿+ℎ𝐿𝑑𝑒𝑣(4) 5 where 𝑘𝑈𝐴𝐻𝑀 [−] is the corrected proportional value in the Nusselt equation (2), 𝑘𝑈𝐴𝐼𝐷𝐸𝐴𝐿[−] is the ideal proportional value in the Nusselt equation, 𝑘𝑈𝐴𝑑𝑒𝑣[−] is the correction function for the proportional factor in the Nusselt equation, ℎ𝐿𝐻𝑀[𝐽∗𝑚𝑜𝑙−1] is the corrected value of the liquid enthalpy, ℎ𝐿𝐼𝐷𝐸𝐴𝐿[𝐽∗ 𝑚𝑜𝑙−1] ] is th e ideal value of the liquid enthalpy, ℎ𝐿𝑑𝑒𝑣[𝐽∗𝑚𝑜𝑙−1] is the correction function for the liquid enthalpy. We hypothesi se that t he heat transfer coefficient deviates because of the effect of the liquid volume change on the mixing efficiency . For this reason, 𝑘𝑈𝐴 was computed as a function of the liquid volume inside the reboiler. The deviation of the liquid enthalpy from the first-principle model was computed as a function of the concentration of Ac and the temperature of the liquid within the reboiler . In this case , the ℎ𝑑𝑒𝑣𝐻𝑀 is not only the enthalpy of the mixture containing Ac, DCM and MeOH. However, it includes 1) the enthalpy excess of the mixture Ac -DCM, 2) the interaction between the solvents an d the API , and 3) the deviation from the ideality of the model. MRF and polynomial functions were employed to model the abovementioned deviations . The explicit expressions of the functional form s are reported in Table 1. In the equation contained in this table, 𝑉[𝑚3] is the volume of liquid within the reboiler, 𝑇 [𝐾] is the temperature within the reboiler, 𝑥𝐴𝑐[𝑚𝑜𝑙 ∗𝑚𝑜𝑙−1] is the molar -fraction of Ac within the liquid contained in the reboiler, 𝑘𝑈𝐴𝐼𝐷𝐸𝐴𝐿,𝐴1,𝐵1,𝐴2,𝐵2 are the fitting coefficients identified from the optimi ser Table 1: Explicit form of the functional structure used to perform the black -box function estimation. MRF Polynomial 𝑘𝑈𝐴𝐻𝑀 𝑘𝑈𝐴𝐼𝐷𝐸𝐴𝐿+𝐴1∗𝑉 1+𝐵1∗𝑉 𝑘𝑈𝐴𝐼𝐷𝐸𝐴𝐿+𝐴1∗𝑉+𝐵1∗𝑉2 ℎ𝐿𝑑𝑒𝑣 𝐴1∗𝑥𝐴𝑐+𝐴2∗𝑥𝐴𝑐∗𝑇 1+𝐵1∗𝑥𝐴𝑐+𝐵2∗𝑥𝐴𝑐∗𝑇 𝐴1∗𝑥𝐴𝑐+𝐴2∗𝑥𝐴𝑐∗𝑇+𝐵1∗𝑥𝐴𝑐2+𝐵2∗𝑥𝐴𝑐2∗𝑇2 2.5 The search algorithms Differential evolution (DE) (Storn and Pric e, 1997) and particle swarm optimi sation (PSO) (Kennedy and Eberhart, 1995) were employed to identify the coefficients of the black -box models. Both the algorithms were trained with 20 particles initially distributed on the coefficient space using a Latin hypercube sampling strategy. 30 iterations were used to perform the training of the models. Each search was repeated 10 times with different positions of the initial population. A more detailed description of the search algorithms is in the SI section 1.4. For the DE, the best/2/bin was used as mutation strategy with a mutation coefficient F=0.5. The cross -over probability factor was set to be 0.7. For the PSO, the inertia coefficient 𝜔=0.85, the individual coefficient 𝑐𝑝=2 and socia l coefficient 𝑐𝑔=1. 2.6 The training approaches Two different training approaches were used to train the coefficients of the two statistical models (𝑘𝑈𝐴𝐻𝑀 and ℎ𝐿𝑑𝑒𝑣) with the different optimi sers. The statistical models were trained simultaneously (Parallel approach) or consecutively (Serial approach). The serial approach was possible due to the structure of the statistical models contained in the hybrid model. The ℎ𝑑𝑒𝑣𝐻𝑀 model does not have any impact during the first evaporation (see Section 2.1) since the concentration of Ac in this phase is zero. For this reason, the 𝑘𝑈𝐴𝐻𝑀 model could be trained separately on th is first step of the process. Once its coefficient s were identified, the ℎ𝑑𝑒𝑣𝐻𝑀 model w as trained on the constant volume step of CSS (see Section 2.1) . The authors hypothesi se that the performance of the training approaches will change with the functional form and the optimi ser employed during the training. In addition, the authors hypothesi se that the serial approach will return model s with higher prediction accuracy and 6 robustness. Eight different combination s of training approach es, optimi sers and functional form s were tested and compared ( Table 2). Table 2: Abbreviation of the eight combination s of training approach, optimi ser and functional form tested and compared in this work Abbreviation Approach Optimi ser Functional form PSO/MRF/Parallel Parallel PSO MRF PSO/poly/Parallel Parallel PSO Polynomial DE/MRF/Parallel Parallel DE MRF DE/poly/Parallel Parallel DE Polynomial PSO/MRF/Serial Serial PSO MRF PSO/poly/Serial Serial PSO Polynomial DE/MRF/Serial Serial DE MRF DE/poly/Serial Serial DE Polynomial The training of the statistical models was performed utilising two batches ; the fitness functions were calculated for each of them , and then the two values were averaged and fed back to the optimi ser. The training loop and the variables used in it are shown in Figure 2. Figure 2: Hybrid model training structure. The coefficients set proposed by the optimi ser is used to run the simulation driven by experimental input. The model returns the estimated time profile of the output variables. These are compared with the experimental output from the plant , and the fitness function is calculated using these deviations. The calculated value of the fitness function is given to the optimi ser, which will propose a new set of coefficient s to test . 2.7 Fitn ess function The fitness function was composed of the most representative and stable variables of the process. Their experimental values were compared to the model predictions at each time step. In case of premature interruption of the DAE resolution, the fitness function returns 1014; on the contrary , the equation (5) is applied to compute its value . 𝐶𝑓𝑖𝑡𝑛𝑒𝑠𝑠 (𝝓)=𝐶𝑇𝑒𝑚𝑝𝑅𝑒𝑏 (𝝓)+𝐶𝑅𝑒𝑐𝑇𝑎𝑛𝑘 (𝝓)+1𝑒10∗ 𝐶𝑇𝑖𝑚𝑒 (𝝓) (5) where 𝝓 is the vector associated with the statistical model coefficients, 𝐶𝑇𝑒𝑚𝑝𝑅𝑒𝑏 (𝝓) is the fitness function part associated with predicting the temperature inside the reboiler ; it is computed as the mean squared error (MSE) of the temperature. 𝐶𝑅𝑒𝑐𝑇𝑎𝑛𝑘 (𝝓) is the fitness function part associated with predicting the volume inside the receiver tank . This variable has a high uncertainty due to process - related unquantified losses ( e.g. sensor noise, vent gas purge, sloshing ); 𝐶𝑅𝑒𝑐𝑇𝑎𝑛𝑘 (𝝓) is thus computed taking into account these uncertainties as described in SI section 1.5. 𝐶𝑇𝑖𝑚𝑒 (𝝓) is fitness 7 function part associated with the prediction of the end time of the process . A more detailed description of the cost function is available in SI section 1.5. 2.8 Definition of the succes s criteria The prediction accuracy of the hybrid models was evaluated utili sing the fitness function values returned from each training run. Since each model was trained 10 times , the median and the variance of the fitness function values were employed as descriptor s for each training group. Lower values of the variance indicate higher robustness of the training technique to the initial conditions. Lower values of the median indic ate the better prediction accuracy of the technique . A twofold approach was applied to evaluate the quality of the black -box models that drive the parameters of the first -principle model . The quality of the models that define 𝑘𝑈𝐴𝐻𝑀 is evaluated utili sing both the physical boundaries suggested in the literature (Mohan et al., 1992) and a reference trend proposed by the DE/poly/Serial . The ℎ𝑑𝑒𝑣𝐻𝑀 statistical model was evaluated only utili sing the trend proposed by DE/poly/Serial . It was used as a reference for two main reasons: 1) in the preliminary tests , it return ed both the more robust fitness score median and variance in moving from the training set to the test set and 2) the statistical models estimated by this technique are less sensitive to the initial population distribution. The index of dispersion (IoD) (6) was utilised to assess the quality of the black -box models returned by the various techniques and training approach 𝐼𝑜𝐷 =𝜎2 𝜇(6) where 𝜎 is the standard deviation of the training group and 𝜇 is the average of the training group. For each training group, the 𝐼𝑜𝐷 value was calculated at fixed independent variables of the black -box models. The average and maximum values of this descriptor are used as indicator s: the lower their values , the more robust the technique with the initial population. 3. Results and discussion 3.1 Parallel training 3.1.1 Comparison of the complete hybrid model prediction in the parallel training In the parallel training , the hybrid models based on the polynomial structure show higher accuracy and robustness than the MRF ( Figure 3a). The presence of the denominator in the MRF structure increases the non -linearity of the coefficient search and narrows down the validity domain of the functions. It led to some runs to numerical instabilities or premature interruption of the hybrid mod el on the test set. In the polynomial form, the absence of the denominator increases the range of the solution, decrease s the non -linearity of the search and reduce s the instabilities of the system. All these factors allow the polynomial -based models to ha ve higher accuracy and robustness than the MRF -based models ( Figure 3). The optimi ser in the parallel training has a marginal effect on the overall hybrid model perfor mances. In the literature , the DE algorithm is well -known to be more robust and performant in the coefficient search than the PSO (Civicioglu and sBesdok, 2013; Vesterstrom and Thomsen, 2004; Zhang and Wei, 2014). However, in this case , its results do not differ so much from the one of the PSO . We hypothesise that t he robustness of the DE is not enough to stabili se the solution proposed by the MRF statistical functions for the case study . The models obtained using DE/poly/Parallel return a lower median value of the fitness function on the training and test sets than the models trained with the PSO/poly/Parallel. This is possible because the DE utili ses a more complex search procedure than t he PSO . The DE has the possibility to both crossover and mutate the members of the population , facilitating the search and allow ing escape from local 8 minima (exploration capabilities). On the other hand, the PSO has a simpler particle selection algorithm based only on the position of the particle and its score. It decreases the possibility of escaping from local minima (Stacey et al., 2003) . Therefore, the search performed by the PSO is more inclined to examine the convexity in which the current best particle is allocated. This behavio ur makes th e search executed with the DE more incline d to overfit the training set when the MRF is employed as the statistical function within the hybrid model . On the contrary, when the polynomial functional form is employed as statistical function, this behavio ur increase s the prediction capabilities of the obtained hybrid model. Figure 3: a) Boxplot of the fitness function for diffe rent approaches, optimi sers and functions obtained on the train set. b) Boxplot of the fitness score for different approaches, optimi sers and functions obtained on the test set. In these plots , numerical instabilities are detectable when the fitness functi on values are equal to 1014. For the plot reporting the scores on the test set, the overall view is reported at the top, and the zoom of the lower fitness functions area is on the bottom . The polynomial function and the serial approach are the better performing techniques in this study. They return the lower median values and the lower variance on the test set. No significant differences are detectable in changing the optimi ser. However, the usage of D E slightly increases the model accuracy and robustness. 3.1.2 Comparison of the black -box model parameter s prediction in the parallel training From Figure 4 and Table 3, one can observe that the polynomial functions show higher robustness with the initial condition and higher identification accuracy than the MRF when predicting 𝑘𝑈𝐴𝐻𝑀. The MRF structure shows numerical instabilities. The role of the optimi ser is neglig ible when the MRF s are used as the statistical functional form. On the contrary, when the polynomial structure is used, the usage of DE increases the robustness of the prediction s (Table 3) and their accurac ies. 9 Figure 4: Profile of the 𝑘𝑈𝐴𝐻𝑀 predicted with the parallel approach with the various combination of optimi ser and black -box structure . MRF return poor prediction capabilities and low robustness since the numerical instabilities and wrong trend. The usage of the polynomial function returns better prediction capabilities than the MRF. In this case, when DE is used as the optimiser , the parameter estimation shows higher robustness with the initial population. From Figure 5, one can observe the lower prediction accuracy and robustness of the MRF structure. It returns high ly scattered solutions with different trend s without any physical meaning , such as 𝑘𝑈𝐴 values lower than 0 . The solutions delivered by the PSO/MRF/Parallel do not follow the trends of the reference for both 𝑘𝑈𝐴𝐻𝑀 and ℎ𝑑𝑒𝑣𝐻𝑀, furthermore the predicted values of 𝑘𝑈𝐴𝐻𝑀 are outside the acceptable physical boundaries , or they show negative values . Despite these deficiencies, the fitness score of the PSO/MRF/Parallel on the training set is still low both in variance and median since the technique overfits on the training set. The overfitting is due to the synergistic training of the two models (parallel). To compensate for the deviation brought from one black -box model (e.g. 𝑘𝑈𝐴𝐻𝑀), the other (e.g. ℎ𝑑𝑒𝑣𝐻𝑀) should adapt its values to reduce the value of the fitness function. The search carried out by the PSO is more affected by this synergistic effect as none of the solutions proposed by the PSO follow s the reference trend , and most of the 𝑘𝑈𝐴𝐻𝑀values are negative. The usage of the DE optimi ser slightly mitigates the effect of the non -linearity brought by the MRF , yet the majority of the solutions proposed by the DE/MRF/Parallel do not follow the reference. Only the models associated with the random state 8 return a solution that follow s the reference tr end. This unreliability of MRF trained with parallel training limits its usage and, in this case, prohibits it. The usage of the polynomial increases the accuracy and robustness of the solution. When the polynomial structure is employed, both the optimi sers return functions that follow the reference. However, the DE is more robust to the initial population than PSO because of the lower scattered profile s of the black -box estimated parameters . Using the polynomial structure drastically increases the capabili ties to identify the coefficient s of the black -box model. It allows the domain of the black - box models to be ℝ. Thus, the majority of the solutions proposed by the PSO are now able to follow the reference trends in 𝑘𝑈𝐴𝐻𝑀 values except for the solution associated with random states 3, 4 and 7 (Figure 4). This a lso affect s the quality of the models that describe ℎ𝑑𝑒𝑣𝐻𝑀 that now have the same trend of the reference. However, the PSO/poly/Parallel approach return s more scattered solutions than DE/poly/Parallel. In Table 3, the average and the maximum values of the IoD related to the 10 PSO/poly/Parallel technique are higher than the DE/poly/Parallel. This behavio ur is due to the increased robustnes s of the DE algorithm to the initial population compared to the PSO optimi ser (Vesterstrom and Thomsen, 2004) . All the black -box models obtained with the DE follow the trend of the reference , and all the values of 𝑘𝑈𝐴𝐻𝑀 are within the literature boundaries. Figure 5: Pro files of the ℎ𝑑𝑒𝑣𝐻𝑀 predicted using the parallel approach with the various combination of optimi ser and black -box structure calculated at 50°C. The polynomial functions show higher robustness and accuracy than the MRF. The models obtained with the MRF show numerical instabilities and poor accuracy regardless of the used optimi ser. The polynomial function shows higher prediction accuracy and robustness with the random state with both the optimi sers. The superior prediction accuracies and robustness offered by the DE/poly compared to the other techniques is related to a combination of several factors: 1) the DE search is very robust to the variations in the initial population; 2) the DE has better capabilities in finding the global minimum compared to the PSO (Vesterstrom and Thomsen, 2004) and 3) the coefficient search for the polynomial structure is linear and less complex to perform compared to MRF. Table 3: The IoD average and max absolute value for the various techniques and training approach es. The calculation is made at the same independent variable values for the different random states Parallel training Serial training 𝒌𝑼𝑨𝑯𝑴 𝒉𝑳𝒅𝒆𝒗 𝒌𝑼𝑨𝑯𝑴 𝒉𝑳𝒅𝒆𝒗 Technique Average Max Average Max Average Max Average Max PSO/MRF 636.77 2.87e4 7.01 7.23 1.39e -2 2.80e -2 449.02 687.13 PSO/poly 0.43 0.74 194.49 592.99 2.96e -3 4.24e -3 86.58 217.71 DE/MRF 3.15 566.66 3.39e3 7.78e3 2.01e -2 3.47e -2 466.34 1328.65 DE/poly 8.57e -2 1.98e -1 153.75 266.02 6.16e -4 1.95e -3 5.24 7.14 3.2 Serial training 3.2.1 Comparison of the hybrid model prediction in the serial training Due to mathematical instabilities, none of the hybrid models trained with the serial training approach show any premature interruption on the test set (Figure 3) . The PSO/MRF/Serial returns a higher median fitness score on the training set than the PSO/poly/Serial , thus a lower accuracy on the training 11 set. However, on the test set , the result is t he opposite favo uring the MRF from the point of accuracy . This is related to the synergic effect of the lower exploration capabilities of the PSO technique and the non-linearity in the coefficient search provided by the MRF structure. These two factors do not allow the PSO/MRF/Serial to overfit the training set and have a better score on the test set compared to the PSO/poly/Serial. The functional form marginalises when the DE is used as an optimi ser. DE/MRF/Serial and the DE/poly/Serial return comparable variance and median fitness scores . This is linked to the robustness of the DE technique that mitigates the non -linear effect s in the coefficient search of the MRF and let s it behave as a polynomial. 3.2.2 Comparison of the black -box model coefficient prediction in the serial training All the techniques , when trained with the serial approach , return predicted 𝑘𝑈𝐴𝐻𝑀 values within the boundaries proposed in the literature . The solution s obtained with the PSO technique are less accurate and robust than those proposed by the DE (Figure 6 and Table 3). All the solutions of ℎ𝑑𝑒𝑣𝐻𝑀 proposed from this training approach return values that follows the reference trend (Figure 7). The MRF show s lower robustness compared to the polynomial . In Table 3, their IoD is one order of magnitude higher than the one associated with the polynomial . This is linked to the linearity of th e coefficient search in the polynomial function; it makes the coefficients search less complex to perform . The usage of DE has higher robustness than the PSO because of the more complex algorithm that the DE uses to find the solution. Like the 𝑘𝑈𝐴𝐻𝑀, also for ℎ𝑑𝑒𝑣𝐻𝑀 the models proposed by the PSO are more scattered than the solution s proposed by the DE when using the same black -box form. In addition, using the same optimi ser, the polynomial functions show fewer scattered solutions than the M RF (Table 3). The serial training results in the lower median and variance fitness function in the test set values when the solution is searche d with the PSO/MRF . However , this technique is not optimal for predicting the coefficient black -box model where the DE/poly return s the best coefficient search. Figure 6: Predicted profile of 𝑘𝑈𝐴𝐻𝑀obtained with the serial app roach with various combinations of optimi sers and black -box structure. All the solutions are within the acceptable range and follow the correc t trend of the predicted variables. However, 12 the polynomial show higher robustness than the MRF since the lower scattering of the solutions. The DE report higher robustness to the initial condition in the optimi sation than the PSO. Figure 7: Predicted profiles of ℎ𝑑𝑒𝑣𝐻𝑀 obtained with the serial approach using various combination s of optimi sers and statistical functions . The solutions are calculated at 50°C. All the solutions follow the reference trend and show comparable prediction accuracy. The MRF show s lower robustness with the initial condition than the polynomial function. T he DE increase the robustness of the solution. The Poly/DE show the lower scatting degree of this work and the higher robustness with the initia l populations. 3.3 Comparison between the two training approaches In general, serial training shows a lower variance and median of the fitness function for all the technique s (Figure 3). It increases both the accuracy and the mathematical stability of the solut ions. In addition, it reduces the influence of the initial population on the model training for the DE and PSO. The serial training has higher accuracy and robustness than the parallel approach in predicting the behavio ur of the overall system and the tre nd and the values of the black -box model. The primary deficiency observed in the parallel training is the mathematical instabilities in the test set linked to the synergistic effect of the two models in the coefficient search. The parallel training perform s the coefficient search in a space with seven degrees of freedom. Moreover, these variables belong to different sub -models (enthalpy and heat transfer coefficient models ). This increases the complexity of the system , and the overall fitness function can h ave more local minima compared to the serial training approach. On the contrary, the serial approach performs the coefficient search at first on a coefficient space with three degrees of freedom during the first evaporation (see Section 2.1) , then on a coefficient space with four degrees of freedom during the constant volume phase (see Section 2.1) . In addition, the coefficient search in fluences only one black -box model in this case . Given that, the DE performs a better coefficient search comp ared to the PSO thanks to its augmented exploration capabilities , allowing the search to escape from local minimum easily . The black -box functional form helps in decreasing the complexity of the cost function shape. The usage of the MRF increases the non -linearity due to the presence of the coefficient s at the denominator. On the contrary, the polynomial has a linear coefficient search. Therefore, using the polynomial black -box structure decreases the complexity of the fitness function with 7 degrees of freedom and makes the coefficient search less complex . 13 A step toward reducing the complexity of the fitness functions is splitting the train ing of the two statistical models in to two different moments. In this case, the coefficient search is more accurate than the one carried out with the parallel training. Here , all the techniques converge to a solution and return models with high accuracy and robustness. This is possible thanks to the reduced complexity of the fitness function that return s a local minimum from which is easier to escape. The technique that better performs in the serial training is the DE/poly. In this case , there is a combination of the high exploration capabiliti es and robustness of the optimi sation algorithm, the reduced non -linearity provided by the polynomial structure and the reduced complexity obtained with the training splitt ing. However, the performance difference between the serial and the parallel trainin g is negligible in this cas e. In other words, parallel training performs as well as serial training when DE/poly is employed as the training technique . 3.4 Extrapolation capabilities of the trained models As discussed in section 2.2, t he extrapolation capabilities of the model were evaluated on a n extra batch . It was executed utilising half of the volume compared to the batches on which the model s have been trained. From Table 4, it can be detected how the hybrid model outperforms a case -optimi sed first-principle model. In addition, from this table , one can see how serial training has superior results compared to parallel training. However , the difference is highly remarkable when PSO and MRF are employed to train the model. On the contrary , the accuracy of the models obtain ed by DE and polynomial function has little change when the training approach has changed. Table 4 proves how the hybrid model developed with the techniques reported in this study can generali se and extrapolate very well on the variable area that was never experienced. This is a n essential step when dealing with data that lacks variety. Table 4: Determination coefficient (R2) of the profiles within the reboiler predicted for the batch obtained with half of the volume than the b atches used for the training. Reboiler Temperature Reboiler Volume Technique Parallel Serial Parallel Serial PSO/MRF 0.973 0.976 0.781 0.881 PSO/poly 0.986 0.986 0.894 0.916 DE/MRF 0.974 0.986 0.789 0.926 DE/poly 0.986 0.986 0.931 0.918 FP model 0.975 0.707 4 Conclusions This paper introduces a novel approach to develop ing hybrid models for chemical processes. The methodology aims to predict the behavio ur of chemical systems under uncertainty (see SI section 1.5 and SI section 1.2 ). The proposed technique can work with data that lacks variety and value. Moreover, the resulting model ha s been shown to have enhanced extrapolation capabilities . The study case utilised in this paper is the constant volume solvent swi tch; however, the proposed methodology applies also to reactive systems and continuous processes. Besides developing a robust hybrid model, the methodology shows promising capabilities in estimating the physical coefficient s of the system and letting their value s be within acceptable physical boundaries. This work investigated the performances of two heuristic optimi sers (DE and PSO), two black -box functional structures (MRF and polynomial) and two training approaches (Serial and Parallel). The hybrid model s based on the polynomial functions achieve d better results than the MRF . The linearity of the polynomial form with the coefficient allows a more accurate coefficient identification than the MRF . The MRF structures report mathematical instabilities and overfitting issues when trained in 14 parallel. On the contrary, the hybrid models based on the polynomial functions do not show any mathematical instabilities. Although the serial approach usage reso lved the mathematical stability issues , the prediction capabilities of the models based on MRF do not increase. The role of the optimi ser is marginal to the hybrid model accuracy. Even though the accuracy of the solutions proposed by the DE is higher than the one obtained with the PSO, the quality of the two models is comparable. On the contrary, the optimi sers play a crucial role in estimating the coefficients of the black -box functions . The coefficient searches obtained with the PSO are more scattered an d sensitive to the initial population positions than the ones proposed by the DE. For these reasons, the authors suggest using the PSO when there is no need to perform physical parameter estimation and the quality of the hybrid model is evaluated only on i ts final predictions. By contrast, the DE is suggested when the black -box coefficient estimation drives crucial information about the process. The training approach has a crucial impact on the quality of the solutions. The serial approach stabili ses the b lack-box solutions makin g them less sensitive to the initial population of the optimi sers compared to the parallel training. This effect is linked to the lower complexity of the benchmark function when training the system in a serial manner. The hybrid mod el obtained with the serial approach employing the polynomial function as a statistical function trained with the differential evolution algorithm (DE/Poly/Serial) returns a hybrid model with both the highest prediction capabilities and the most stable black -box function estimation. The high robustness of the optimi ser, the linearity of the coefficient search and the serial training reduce the complexity of the cost function . This combination of search algorithm /model/ training scheme 1) reduces the dispersion of the black -box solutions , 2) let s the estimated parameters be with in the acceptable physical boundaries and 3) returns a hybrid model able to predict the system behavio ur with high accuracy . Acknowledges This work was supported by the DA P2CHEM ICON project funded by VLAIO and Catalisti (VLAIO HBC.2020.2455). The resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government. Competing Interests KU Leuven submitted a patent application including these results on 16th December 2022 to the European Patent Office. The patent application number is EP22214177.2. Ulderico Di Caprio, Min Wu, Furkan Elmaz, Siegfried Mercelis, Peter Hellinckx and M. Enis Leblebici are included in the inventors list.