Regression Metamodel-Based Digital Twin for an Industrial Dynamic Crossflow Filtration Process

1. Introduction

Since the term Industry 4.0 was introduced at the Hanover Fair in 2011, the importance of the development of cyber-physical production systems (CPPS) is growing faster than ever also in the bio-manufacturing industry. CPPS are physical production plants equipped with computation processing units used to automatically control the process. The realization of a cyber-physical system for manufacturing is generally understood to require, among other things, a virtual replica of the process that bidirectionally interacts with the physical process, the collection and usage of process data and the optimization of the process with the help of intelligent algorithms [1]. The virtual part of a CPPS including the digital process model and the algorithms is called digital twin (DT) [2,3].

Concepts and requirements for digital twins in manufacturing have frequently been presented and reviewed in literature [4,5,6,7,8,9,10,11,12,13,14,15], whereas practical use cases where a digital twin has been developed for real-life processes are few to be found in literature. Lopez and co-workers built a hybrid model for a fermentation process that is updated with data during the process [16], but they did not close the information loop from the digital part back to the physical process. A digital process replica with such a “one-way” information flow from the physical to the digital part is called a digital shadow [4,16]. Trunzer et al. used data-driven models trained with expert knowledge to identify the operational states of an industrial chemical reaction plant from the process data [17]. However, the models are not used for plant process optimization. Chew et al. used data-driven models to control a filtration process [18]. The optimization algorithm was tested on a laboratory filtration unit and a benefit for control of industrial filtration units is predicted.

Digital twins rely on digital process models by definition, the backbone that makes process monitoring, simulation and optimization possible [19,20,21]. Regarding the literature, mathematical models used to describe bioprocesses are divided into mechanistic and data-driven models [19,22]. Mechanistic models include models based on first principles namely mass and energy balances and previous process knowledge. Data-driven models are built up on process data without incorporated previous knowledge. The trend of the recent years is to combine mechanistic and data-driven models to hybrid model structures in order to profit from the extrapolation capability of mechanistic models and the flexibility of data-driven models at the same time [22,23]. Mechanistic and data-driven models can be arranged in different serial or parallel configurations to hybrid models, depending on which modeling goal is to be aimed for [24,25]. The most commonly used model configuration in biochemical engineering is a serial structure of a data-driven model and a mechanistic model wherein the data-driven model functions as filler for overall model parts that cannot be expressed based on first principles [24].

Regarding the model description of filtration processes, a spectrum of different approaches and use cases is presented in literature [26]. Krippl and coworkers connected a neural network with mass balances in a serial model configuration to predict the permeate flux over time in a conventional crossflow filtration unit [27]. Piron et al. applied the serial model configuration to yeast filtration [28]. Díaz et al. used a serial neural network and Darcy law configuration for flux prediction, but did not use the model for process control and optimization [29]. Chew et al. used a similar model structure for the prediction of fouling and filtration resistances of an industrial water clarification process [30]. Also here, the model was not yet further used for process optimization.

A special type of model interconnection is the concept of metamodels or surrogate models [19,22,31,32]. A surrogate model is defined as a superordinated, simplified model based on the knowledge of the original, underlying model. The metamodel reduces complexity and increases calculation time to enable model calling during real-time optimization of a process. Franzreb et al. used a model-based design-of-experiments (DoE) approach to make statements about the economic variables of a complex antibody production process [33]. Wang et al. trained a neural network with data generated from a mechanistic model of protein chromatography and he recognized advantages of the metamodel to quickly predict model parameters to adapt the mechanistic model to new feed streams [34]. Reports on a metamodel built for a filtration use case could not be found in the literature.

This paper is aimed to focus on the process of dynamic cross-flow filtration (DCF). DCF is a high performance filtration technique that relies on rotating discs reducing the deposition of fouling components on the membrane surface [35,36]. Modeling of rotating membrane discs requires consideration of physical phenomena occuring only with DCF such as cake reduction and fluid backpressure due to rotation of the membranes [37,38]. In this publication, the reader is shown how we built a mechanistic, semi-empiric, hybrid process model of membrane fouling during DCF of grape must. We present a concept of training a data-driven metamodel banking on model-based DoE and response surface modeling and showcase how to straightforwardly apply the resulting metamodel for the time-critical forecasting of optimum process productivity. The model collective is implemented into a digital twin framework that reads the torque and TMP parameter at regular time intervals from the physical DCF plant during the process run and returns the permeate flow rate setpoint calculated for optimum productivity to the process actuator. The added value of introducing the digital twin to the industrial process is examined by comparing the overall productivity between 2022 and 2023 grape must filtration campaigns.

2. Materials and methods

In this section, the plant set-up of the industrial filtration runs and an examplary set of process data, that is applied for model construction, are presented. The structure and the elements of the digital model collective are explained, followed by a description of the superstruction of the digital twin framework.

2.1. Production set-up and operation

The experimental set-up for the filtration runs is shown in Figure 1. The filtration machine (F1) is a DCF type number 312/32 (Andritz Separation GmbH, Vierkirchen, Germany) equipped with four shafts and 256 ceramic membrane discs. Permeate passes through the porous membrane discs and it is discharged from the machine inside the hollow shafts. Retentate is collected inside the process chamber and can be released by opening a valve (V1). The membrane discs are mounted on top of each other on every shaft, so that they form membrane stacks. Adjacent membrane disc stacks overlap, causing increased shear rate at the overlap region, thus achieving an increased cleaning effect from solid deposits at the overlapping points. The total available filter area is 32.8 m², the membrane diameter is 312 mm.

The feedstock considered for the filtration experiments of this research study is pressed grapes containing pulp and other solids. Thus, the valuable product of the process is the clarified sweet must, i. e. the permeate. A membrane pore size of 200 $\mathrm{n}$$\mathrm{m}$ is used for the filtration process. The filtration runs were performed at a production facility of a wine manufacturer in Italy during the regular sweet must production campaigns from August to October 2022 and 2023.

A centrifugal pump (P1) transports the feed from a storage tank into the DCF chamber. The feedstock is periodically mixed in the storage tank, however, a fully uniform solids content in the feed during the entire filtration process cannot be ensured. The pumping rate of the feed pump is controlled by a programmable logic controller (PLC) using the signal of sensors recording transmembrane pressure (TMP) and permeate flow rate. TMP is measured by the difference between feed inlet pressure and permeate outlet pressure, i. e. TMP reduction due to membrane rotation is neglected. The filtration process starts in “flow controlled” operation mode, that is pumping rate of the feed pump is controlled for a constant permeate flow rate. The initial permeate flow setpoint is predefined by the plant operator at the beginning of the process run from experience dependent on estimated solids concentration in the feed and post-cleaning permeability of the membrane discs. As soon as the TMP setpoint of 0.8 bar is exceeded once during the production run, operation mode automatically switches to “TMP controlled” and the feed pump rate is controlled for keeping the TMP at a constant value.

The rotation speed of the membrane drive is controlled to a constant value of 340 min⁻¹ during all experiments. Readout of power consumption from the frequency converter of the motor rotating the membrane shafts is done to receive data about the membrane drive torque. The membrane drive torque is strongly correlated to the retentate viscosity inside the DCF chamber and provides information about the increase of solids concentration in the retentate during a filtration run. Opening of the retentate valve for release of retentate from the DCF chamber is controlled by the membrane torque. At the beginning of the process, the DCF chamber is filled with feed suspension. During concentration phase, the retentate valve is kept closed until the torque setpoint of 300 N m is reached. Subsequently, retentate discharge begins by gradual opening of the retentate valve until the torque does not increase any further. Closing of the retentate valve happens with a torque hysteresis of 10 %, i. e. it closes only when membrane torque falls below 270 N m.

The end of a filtration run is decided by the operator based on the amount of feed suspension available. A cleaning procedure with variable efficiency is followed by each filtration process. The permeate water flux is measured after every cleaning cycle for the assessment of the cleaning efficiency. All sensor and actuator setpoint data is logged 24/7 in 5 $\mathrm{s}$ intervals as time series in a server database and was used for model and digital twin validation after download in csv format and segmentation into single filtration runs.

2.2. Mechanistic-empirical process dynamics model

A dynamics model that describes the physical process as close as possible is the backbone of the digital twin. We decided to build a mainly mechanistic model complemented with semi-empirically chosen adjustment parameters in order to use existing process knowledge as a counterbalance to the limited variations in the available process data from filtration runs in industrial production.

The calculation of the process parameters TMP, permeate, retentate and feed flow rates, membrane drive torque and valve opening percentage is implemented in a time loop with a cycle time of $3.6$ $\mathrm{s}$. In order to replicate the solids concentration gradient along the shaft axis from the feed inlet to the retentate outlet, the DCF chamber was conceptually divided into segments and the permeate flow rate and membrane drive torque were calculated separately for each segment. The number of segments was set to $s_{\max }=$ 4 after a segment parameter study and the comparison of model simulation accuracy. TMP was not calculated per segment, since pressure drop along the shaft axis was assumed to be negligible. The sequence of calculation steps of the process dynamics model is presented in Figure 2 in the right-hand detail box.

The Darcy equation was used to obtain permeate flow rates during TMP-controlled operation mode or TMP during flow-controlled mode, respectively:

$$\mathrm{TMP}\left(t_i\right)=\frac{\eta ·Q_{\mathrm{per}}(t_i,s)·R_{\mathrm{tot}}(t_i,s)}{A_{\mathrm{segm}}}$$

(1)

where $\eta$ is the dynamic viscosity of water in $\mathrm{Pa}\mathrm{s}$, $Q_{\mathrm{per}}$ is the permeate flow rate in L h⁻¹, $R_{\mathrm{tot}}$ is the total filtration resistance in m⁻¹ and $A_{\mathrm{segm}}$ is the filtration surface area of one segment. The dependencies of the variables on the time $t_i$ and the segment number s are marked. In case of flow-controlled operation mode, the calculation is not trivial, because although the total permeate flow rate is set by the operator or the digital twin, its distribution between the different segments of the DCF model is unknown at this stage of the time step. Therefore, in an intermediate calculation for each segment a reference flow rate was calculated for a reference TMP. Subsequently, the sum of these reference permeate flow rates was compared to the desired total permeate flow rate and the TMP was adjusted in order to match the flow rates. Knowing the required TMP, the real segment-wise permeate flow rates were calculated from the Darcy equation.

The resistance-in-series approach was chosen for the modeling of the filtration resistances:

$$R_{\mathrm{tot}}(t_i,s)=R_{\mathrm{m}}+R_{\mathrm{pore}}(t_i,s)+R_{\mathrm{cake}}(t_i,s)$$

(2)

where $R_m$ is the intrinsic membrane resistance in $\mathrm{m}$, $R_{\mathrm{pore}}$ is the pore resistance in m⁻¹ and $R_{\mathrm{cake}}$ is the cake resistance in m⁻¹. The pore resistance was defined to follow an exponential relation:

$$R_{\mathrm{pore}}(t_i,s)=R_{\mathrm{pore},\mathrm{ref}}·\left(exp\left(k_{\mathrm{pore}}·k_{p,i}(t_i,s)\right)-1\right)$$

(3)

where $R_{\mathrm{pore},\mathrm{ref}}$ is the reference pore resistance in m⁻¹, $k_{\mathrm{pore}}$ a model adaption parameter in L g m⁻¹. The kinetic parameter $k_{p,i}(t_i,s)$ is calculated as follows:

$$\frac{\mathrm{d}{k}_{p,i}(t_i,s)}{\mathrm{d}{t}_i}=c_{\mathrm{fine}}·\frac{Q_{\mathrm{per}}(t_i,s)}{A_{\mathrm{segm}}}$$

(4)

where $c_{\mathrm{fine}}$ is the concentration of fine particles in the DCF chamber in g L⁻¹. Fine particles are assumed to pass through the membrane and fine particle concentration is set to 20 % of the feed concentration. A rational equation of the Langmuir type was used to model the cake resistance striving asymptotically to a maximum during filtration runs:

$$R_{\mathrm{cake}}(t_i,s)=R_{\mathrm{cake},\mathrm{SS}}(t_i,s)·\frac{k_{c,i}(t_i,s)}{k_{c,i}(t_i,s)+k_{\mathrm{cake}}}$$

(5)

where $R_{\mathrm{cake},\mathrm{SS}}(t_i,s)$ is the steady-state cake resistance in m⁻¹ and $k_{\mathrm{cake}}$ a model adaption parameter in g h L⁻¹. The increase of the parameter kinetic parameter $k_{c,i}(t_i,s)$ in time is defined as being dependent on the chamber concentration:

$$\frac{\mathrm{d}{k}_{c,i}(t_i,s)}{\mathrm{d}{t}_i}=c(t_i,s)$$

(6)

and the steady-state cake resistance follows the equation:

$$R_{\mathrm{cake},\mathrm{SS}}(t_i,s)=R_{\mathrm{cake},\mathrm{SS},\mathrm{ref}}·{\left(\frac{Q_{\mathrm{t},\mathrm{per}}\left(t_{i-1}\right)}{Q_{\mathrm{t},\mathrm{per},\mathrm{ref}}}\right)}^{n_Q}·exp\left(\frac{\mathrm{TMP}\left(t_{i-1}\right)}{{\mathrm{TMP}}_{\mathrm{compress}}}\right)·{\left(\frac{c(t_{i-1},s)}{c_{\mathrm{ref}}}\right)}^{n_c}$$

(7)

where $R_{\mathrm{cake},\mathrm{SS},\mathrm{ref}}$ in ${\mathrm{m}}^{-1}$, $Q_{\mathrm{t},\mathrm{per},\mathrm{ref}}$ in L h⁻¹, ${\mathrm{TMP}}_{\mathrm{compress}}$ in bar, $c_{\mathrm{ref}}$ in g L⁻¹ and dimensionless $n_Q$ and $n_c$ are model parameters. Variable values taken from the previous time step of the calculation are marked with index $i-1$. The composition of the steady-state resistance equation results from the analysis of the 2022 campaign’s process data and is based on physical knowledge about cake formation, cake shearing and cake compressibility in filtration.

Fluid mass balances were implemented to obtain the retentate flow rates per segment, in which feed flow rates of a segment is understood as the retentate flow rate of the previous segment:

$$Q_{\mathrm{ret}}(t_i,s)=Q_{\mathrm{ret}}(t_i,s-1)-Q_{\mathrm{per}}(t_i,s)$$

(8)

The overall feed flow rate of the process defines the inlet of the first segment:

$$Q_{\mathrm{ret}}(t_i,0)=Q_{\mathrm{feed}}$$

(9)

Chamber concentration values are calculated from segment-wise, instationary species mass balances:

$$V_{\mathrm{segm}}·\frac{\mathrm{d}c(t_i,s)}{\mathrm{d}{t}_i}=Q_{\mathrm{ret}}(t_i,s-1)·c(t_i,s-1)-Q_{\mathrm{ret}}(t_i,s)·c(t_i,s)$$

(10)

with the following boundary condition:

$$c(t_i,0)=c_{\mathrm{feed}}$$

(11)

where c is the coarse solids concentration in the DCF chamber in g L⁻¹, saying the concentration of particles and colloids which are retained by the membrane. The membrane drive torque is calculated from the solids concentration in two steps, firstly via an empirical equation relating the viscosity exponentially to the concentration:

$$\eta (t_i,s)={\eta }_{\mathrm{ref}}·exp\left(k_{\eta }·f_c\left(t_i\right)·c\left(t_i\right)\right)$$

(12)

where $\eta$ is the dynamic viscosity in $\mathrm{Pa}\mathrm{s}$, ${\eta }_{\mathrm{ref}}$ and $k_{\eta }$ are model parameters and $f_c\left(t_i\right)$ is an empirical correction function that takes into account that the suspended solids concentration in the DCF chamber does not exactly follow Equation 10, possibly because a fraction of the solids is fixed in the filter cake or attached to the chamber walls. In order to consider the effect, which results in a slowed down increase of the torque, $f_c\left(t_i\right)$ follows an exponential decay until the retentate valve is determined to open, indicating the approach of a quasi-steady state in the DCF. Secondly, the well-known first principle mechanistic equation linking viscosity and torque is applied:

$$M(t_i,s)=M_{\mathrm{seal}}+k_g·\omega ·\eta (t_i,s)$$

(13)

where $M_{\mathrm{seal}}$ is the torque already caused by the seals of the rotating shafts and $k_g$ is model parameter determined by geometry of the DCF, mainly the total filter area and the distance between the discs.

The last calculation step is to obtain the total retentate flow rate from the valve opening:

$$Q_{\mathrm{t},\mathrm{ret}}\left(t_i\right)=k_{\mathrm{valve}}·\left(\mathrm{TMP}\left(t_i\right)·X\left(t_i\right)-{\tau }_y\right)$$

(14)

where $Q_{\mathrm{t},\mathrm{ret}}$ is the total retentate flow rate in L h⁻¹, $k_{\mathrm{valve}}$ is the valve cross-section coefficient in L bar⁻¹ h⁻¹, X the valve opening in % and ${\tau }_y$ the flow limit in bar. When the critical torque is reached, the degree of valve opening is gradually increased starting from a minimum opening percentage until the torque starts to descend.

For the simulation of the model, the initial permeate flow rate setpoint, the solid feed concentration, the process duration and all model parameters have to be defined as input. After the digital process simulation, the mean productivity over a filtration run is obtained via a time averaging and space summation of the permeate flow rate:

$$P=\frac{1}{t_{\mathrm{end}}}·{\int }_{t=t_0}^{t_{\mathrm{end}}}\sum _{s=1}^{s_{\max }}{Q}_{\mathrm{per}}(t_i,s)\mathrm{d}{t}_i$$

(15)

where P is the average productivity of the production batch in L h⁻¹, which we simply call productivity in the following, $t_{\mathrm{end}}$ the process duration in $\mathrm{h}$ and $s_{\max }$ the dimensionless number of segments.

2.3. Regression metamodel

A series of in silico experiments was performed with the process dynamics model to obtain the data for building the data-driven metamodel. A fullfactorial experimental design was chosen and the initial permeate flow rate setpoint, the feed concentration and the reference steady-state cake resistance were varied in the boundaries of [400, 1200]L h⁻¹, [20, 100]g L⁻¹ and [6 · 10¹², 2.2 · 10¹³]m⁻¹, respectively. Five levels were calculated per factor. All in silico experiments were carried out with a process duration of 10 $\mathrm{h}$. Subsequently, a second-degree polynomial was fitted to the obtained values of the average productivities. Because of the strongly differing magnitude of their numerical values, the factors had to be scaled to the interval [0, 1] before performing the fit. The calculations were done using Python including pyDOE2 and the preprocessing.PolynomialFeature and linear_model.LinearRegression packages from the scikit-learn library.

Including the first- and second-order parameter terms, the interaction terms and the constant bias term, the resulting response surface equation approximating the achievable productivity in dependence of the process parameters feed concentration, permeate flow rate setpoint and steady-state cake resistance has the form:

$$\begin{array}{cc}\hfill P=w_0& +w_1{c}_{\mathrm{feed}}+w_2{Q}_{\mathrm{per},\mathrm{set}}+w_3{R}_{\mathrm{cake},\mathrm{SS}}+w_4{c_{\mathrm{feed}}}^2+w_5{c}_{\mathrm{feed}}{Q}_{\mathrm{per},\mathrm{set}}\hfill \\ \hfill & +w_6{c}_{\mathrm{feed}}{R}_{\mathrm{cake},\mathrm{SS}}+w_7{Q_{\mathrm{per},\mathrm{set}}}^2+w_8{Q}_{\mathrm{per},\mathrm{set}}{R}_{\mathrm{cake},\mathrm{SS}}+{w_9{R}_{\mathrm{cake},\mathrm{SS}}}^2\hfill \end{array}$$

(16)

where $w_i$, $i\in$ [0, 9], are the parameters of the response surface, that are identified via polynomial regression.

The aim of the digital twin is to optimize the productivity of the DCF. Therefore the extrema of the predicted productivity of the meta model are determined via the first-order derivative with respect to the permeate flow rate setpoint equals zero:

$$\frac{\mathrm{d}P}{\mathrm{d}{Q}_{\mathrm{per},\mathrm{set}}}=w_2+2w_7{Q}_{\mathrm{per},\mathrm{set}}+w_5{c}_{\mathrm{feed}}+w_8{R}_{\mathrm{cake},\mathrm{SS}}\stackrel{!}{=}0$$

(17)

The second-order derivative was checked to verify that the extremum found is a productivity maximum:

$$\frac{{\mathrm{d}}^{2}P}{\mathrm{d}{Q_{\mathrm{per},\mathrm{set}}}^2}=2w_7\stackrel{!}{<}0$$

(18)

Equation 17 allows for real-time calculation of the permeate flow rate setpoint reaching for the maximum productivity of the process. Since the feed concentration is included, the optimum permeate flow rate setpoint can be easily recalculated when a change in the feed composition occurs during the filtration process run. If an unmodeled physical effect in the process causes a mismatch between the real and the modeled process dynamics, the reference steady-state cake resistance parameter is used to adapt the mechanistic model during industrial DCF runs in 5 $\min$ intervals, if necessary. Since also this parameter describing physical properties of the solids in the feed is included as factor in the response surface equation, changes of the natural feed solution sweet must during or between industrial runs do not require the repetition of the DoE pre-studies.

2.4. Digital twin framework

The digital twin framework is built to use the model components and to communicate with the physical asset. The algorithm evaluates the plant status every 60 $\mathrm{s}$ by watching the PLC internal status variable. As soon as a process run is detected, the digital twin framework’s calculations are started. In Figure 2 the flow sheet of the digital twin framework is shown in the left-hand detail box. Once at the process startup the intrinsic membrane resistance is calculated from the permeate water flux according to Darcy Equation 1. Due to detectable cake formation during DCF chamber filling already, the cake resistance parameter $k_{c,i}$ is initialized once at the moment regular filtration starts using measured TMP and permeate flow rate.

During the filtration process, the membrane drive torque is read every 5 $\min$ and a permeate flow rate setpoint recommendation for the optimum productivity is returned to the DCF process control system. In every optimization cycle, the current feed concentration is estimated using the mechanistic-empirical correlation between the feed concentration and the torque from Equation 12 and Equation 13 in the sense of a soft sensor. Secondly, the function based on first order derivative of the regression metamodel from Equation 17 is called and the permeate flow rate recommendation is obtained for the current feed concentration and reference steady-state cake resistance. Thirdly, the mechanistic-empirical model is used to simulate the past 5 $\min$ starting from the simulated status of the DCF obtained by the respective 5 $\min$ simulation in the previous optimization cycle of the digital twin. The simulation used the real experimental permeate flow rate setpoint and the current feed concentration. The resulting TMP from the simulation is compared to the actual TMP in the process. In case the TMP difference between simulation and experiment exceeds $0.1$ bar, the reference steady-state cake resistance is adapted by 10 % for model refinement.

All program parts are written in Python and uploaded on the data server platform of the filtration plant. The platform serves as an interface between the digital twin framework and the local DCF process control system and allows reading of sensor data and writing of actuator instructions during the process. Due to safety reasons, changes of the permeate flow rate setpoints were manually authorized by the operator.

3. Results

The results section is divided into the presentation of the process dynamics model functions, the model validation with industrial experiment data and the results of the metamodeling. The impact of applying the digital twin’s permeate flow rate recommendations on the 2023 process campaign is evaluated in the last section.

3.1. Mechanistic-empirical process dynamics model

A mechanistic-empirical model has been created to calculate the time course of all important process variables. Exemplary results of the predicted process dynamics can be seen in Figure 3. The simulation was defined to start with a permeate flow rate setpoint of 415 L h⁻¹ and the target permeate flow rate is realized by the model system adjusting TMP over time. The TMP increases during the initial $2.5$ $\mathrm{h}$ of the modeled process indicating a progressive membrane blocking. The slope of the TMP decreases during the initial $2.5$ $\mathrm{h}$ of the modeled process because the build-up of the cake filtration resistance, that tends towards a maximum, is dominant over the increase of the pore filtration resistance. The process time $2.5$ $\mathrm{h}$ after the start of the simulated process is characterized by the opening of the retentate valve, therefore TMP increase is stopped and the TMP level is held. The membrane rotation speed is defined to remain constant. At the simulated process time of 5 $\mathrm{h}$, the permeate flow rate setpoint was defined to increase by 80 L h⁻¹ to showcase the switch from flow controlled to TMP-controlled filtration mode. Until the end of the simulated process, the TMP is held at the pre-defined setpoint of $0.8$ bar and the permeate flow rate is calculated to slowly decrease because of increasing pore filtration resistance.

Due to the constant permeate flow rate and a constant feed concentration of 60 g L⁻¹, a linear increase of the lump chamber concentration calculated by simple mass balance is observed until the retentate valve opens for the first time. Thereafter, the simulated concentration decreases and again increases dependent on the opening state of the valve. During the retentate discharge phase, the higher the valve opening percentage is, the higher retentate volume is discharged and thus the faster the chamber concentration falls. The increase and decrease of the membrane drive torque correlates with those of the chamber concentration, following the modeled relationship between torque, viscosity and concentration.

According to the results, the mechanistic-empirical model is capable of a simplified replication of the physical effects of membrane blocking and the control behavior of the filtration plant.

3.2. Process dynamics model validation

During the process runs with the digital twin, the current TMP and torque from the experiment and the model simulation were recorded after every optimization cycle. The plots resulting from two process runs are presented in Figure 4. A visual assessment shows a good agreement between simulation and experiment. Differences between simulated and real TMP are smaller than the measurement noise. The increase in TMP within the first hour of the process run 2-22 and the subsequent stabilization to a constant TMP level is recognized by the model. The initial, real TMP of run 2-22 is met by the simulation with a deviation of less than 10 %. During the concentration phase of run 2-22, the increase in torque is adequately modeled for the initial 2 $\mathrm{h}$ of the process. However, the subsequent gradient of torque increase is limitedly overestimated by the model, so that the setpoint torque for the valve opening is reached 1 $\mathrm{h}$ earlier within the simulation. In run 2-22, the valve opening percentage was precisely set, so that alternating opening and closing of the valve does not occur. This characteristic is properly simulated from the process run time of $5.5$ $\mathrm{h}$. Overall, the torque profile during the concentration phase and the retentate discharge phase is satisfactorily modeled.

Considering the TMP course validation of run 14-23, although the initial TMP is simulated to be 0.15 bar lower than the actual experimental value, the correction function of the digital twin adjusted the simulation within the first 30 $\min$ of the run. The shift in the TMP course at 4 $\mathrm{h}$ results from a change in the permeate flow rate setpoint from 400 L h⁻¹ to 500 L h⁻¹ and is reproduced by the model simulation. The constant TMP of 0.8 bar due to the change to TMP-controlled operation mode at the time of the permeate flow rate setpoint change is simulated correctly. In the torque course of run 14-23, the increase during the concentration phase, the torque level during the retentate discharge phase and the fluctuations due to the changes of the valve opening percentage are adequately approximated by the simulation.

All validation results of 2022 and 2023 campaign’s production runs can be viewed in Figure S1 and Figure S2 in the Supplementary Information. The validated process dynamics model sets the basis for training the metamodel and for the recommendations for optimum permeate flow rate setpoints of the digital twin.

3.3. Regression metamodel

125 in silico experiments were carried out with variation of the feed concentration, the initial permeate flow rate setpoint and the reference steady-state cake resistance. Using the productivity output from the in silico experiments, a regression metamodel was created which allows inverse calculation to the optimum initial permeate flow rate setpoint in the digital twin framework. The relationship between the initial permeate flow rate setpoint and the resulting experiment’s productivity in the metamodel is shown together with the regression residuals in Figure 5. With one exception, the absolute residuals from the polynomial regression are lower than 60 L(per) h⁻¹, indicating a satisfactory goodness of fit of the simple regression metamodel to the complex mechanistic model (Figure 5 d).

All in silico experiments with the initial permeate flow rate setpoint of 400 L h⁻¹ result in a productivity of 400 L(per) h⁻¹ within the accuracy of the regression model. With the initial permeate flow rate setpoint of as low as 400 L h⁻¹, the maximum TMP of 0.8 bar is not reached during the 10 h run and consequently the permeate flow rate could be kept for the complete simulation. The same holds true for the experiments having an initial permeate flow rate setpoint of 600 L h⁻¹ and a reference steady-state cake resistance of 6.0 · 10¹² m⁻¹. However, even in the case of this low cake resistance and the smallest feed concentration investigated of 20 g L⁻¹, this 1 : 1 dependency between productivity and initial flow rate setpoint does not hold anymore when the initial permeate flow rate setpoint is increased to 800 L h⁻¹. At the initial permeate flow rate of 800 L h⁻¹, the maximum TMP of 0.8 bar is reached before the process time of 10 h is finished and the process has to switch to TMP-controlled operation mode. As a consequence, the permeate flow rate starts to decay after this time point and the (average) productivity of the run is less than the initially chosen permeate flow rate setpoint. In Figure 5 a, the effect is shown in an increasing deviation of the course of the curves from a linear correlation of productivity with increasing initial permeate flow rate setpoint. If the assumed reference steady-state cake resistance is set to higher values of 1.4 · 10¹³ m⁻¹ or 2.2 · 10¹³ m⁻¹ (Figure 5 b and c), the time until the TMP setpoint is reached and the system switches to TMP-controlled mode with reduced productivity gets shorter. In addition, the compressibility of the filter cake starts to show a significant negative effect on the productivity at higher cake resistances. The combination of these effects can result in a situation in which the selection of a too high initial flow rate setpoint corresponds with a productivity which is lower than the productivity which could have been achieved if the process had been operated more carefully. Consequently, Figure 5 b and Figure 5 c show the existence of an optimal initial permeate flow rate setpoint.

If the feed concentration increases, e. g. from 20 g L⁻¹ to 40 g L⁻¹, productivity also decreases, however the effect is relatively minor. The reason for the decrease in productivity can be found in the fact that higher solids concentrations accelerate the increase in filtration resistance over the process time. Consequently, TMP increases faster, resulting again in an earlier change from flow-controlled to TMP-controlled operation mode. In most cases, the influence of the feed concentration onto the optimal initial permeate flow rate setpoint is quite small, however at high reference steady-state cake resistances it can be seen that a feed concentration of 100 g L⁻¹ shifts the optimum to lower values of the flow rate setpoint (Figure 5 c).

In summary, the variation of the optimum permeate flow rate setpoint with the model parameters reference steady-state cake resistance and feed concentration confirms the relevance of a metamodel enabling the calculation of the appropriate permeate flow rate setpoint in the digital twin framework during the filtration process.

In case of varying properties of the feed entering the DCF or an adjustment of the mechanistic model caused by an observed significant deviation of the simulated TMP from the experimental data, the metamodel allows to conduct a new calculation of the optimal initial permeate flow rate setpoint in practically real-time. The frequent calculation of the optimal initial permeate flow rate setpoint would not be possible when using the mechanistic model for optimization calculations directly. To find the optimum of the productivity resulting from a 10 h filtration, several runs of the mechanistic model covering the full process time would be needed, resulting in a computation time clearly exceeding the 5 min intervals of the optimization loop of the digital twin.

3.4. Effect of the digital twin on the productivity

By start of the 2023 production campaign, the digital twin was implemented on the server platform and put into operation. At this stage of the development, the permeate flow rate recommendations from the digital twin were transferred manually to the physical process by the operator during the course of the production runs due to safety reasons. The productivities of the 2022 and 2023 campaigns’ production runs and the campaign mean productivity are presented in Figure 6.

The average productivity of the productions runs in the 2022 campaign is 466 L(per) h⁻¹ ± 88 L(per) h⁻¹. An increase in average productivity is achieved from the 2022 campaign to the 2023 campaign, as the average productivity of the 2023 campaign is 536 L(per) h⁻¹ ± 63 L(per) h⁻¹. This shows, that in addition to the increased productivity, a reduction of the fluctuations of the productivity of individual runs could be achieved by implementing the digital twin. The variations in the productivity of the individual runs during both filtration campaigns are explained by different feed characteristics of the respective batches. For instance, solids concentration, fouling particle compressibility and temperature influence the degree of membrane blocking and thus the maximum applicable permeate flow rate setpoint.

The achieved increase in productivity in the campaign applying the digital twin can be mainly attributed to fact that the digital twin suggested initial permeate flow rate setpoints which resulted in DCF operation at higher TMPs. In the manual operation during the 2022 campaign the operators strived to avoid a transition into TMP-controlled operating conditions. However, the simulation runs and metamodel showed, that the reference steady-state cake resistance of sweet must is only moderate and that a transition into TMP-controlled operation mode in the course of the filtration does not harm the productivity. This predictions could be confirmed in the 2023 campaign, as can be seen in the individual TMP time courses of all runs plot in Figure S1 in the Supplementary Information.

4. Discussion

The results are discussed with regard to the quality of the mechanistic-empirical process dynamics model and the digital twin’s faculty to identify the optimum process setpoint.

4.1. Process dynamics model quality

Generally, the mechanistic-empirical process model provides a good replication of the time dynamics of the process variables TMP and membrane drive torque. However, there are effects observed during some process runs that cannot be explained with current process knowledge. For instance, the initial TMP of process run 14-23 in Figure 4 is unexpectedly high and the subsequent TMP course stays constant within 4 $\mathrm{h}$ despite a constant permeate flow rate and a closed retentate valve. Nonetheless, the adaption functionality of the digital twin recognizes the TMP deviations and corrects the model simulation. The metamodel is trained using the fundamental mechanistic-empirical model, i. e. the metamodel does not take the individually observed, unmodeled effects in the production data into account. A further development of the digital twin with regard to the ability to learn on the basis of observed effects in current or past production runs can possibly improve the model.

The development and iterative refinement of the process dynamics model was based on the data of the complete sweet must production campaign in 2022. However, the industrial production runs in 2022 were often carried out with similar setpoint values, giving only limited variation to the process data. We tried to compensate for the limited information available from the data by considering the expert knowledge from project partners and process data obtained with other feed solution types. Nevertheless, the acquisition of data from further industrial filtration runs with variations in e. g. TMP maximum setpoint, permeate flow rate setpoint and membrane drive torque setpoint, followed by the data integration into further model development would be desirable.

4.2. Digital twin optimization capability

The capability of the digital twin to optimize the process depends on the quality of the process dynamics model, but also on the accuracy of the trained metamodel. To generate the training data for the metamodel, in silico experiments were carried out, with the duration of all experiments fixed at 10 $\mathrm{h}$. Originally, the real productions runs were planned to also be carried out with a fixed duration of 10 $\mathrm{h}$ to ensure precise comparability. However, in practice, the duration of the industrial runs is primarily determined by the amount of feed available at the time of the filtration run and it is not set by the operator in advance. Thus, the duration of the filtration runs varies considerably over a campaign. Here, the optimization potential can be further exploited by matching the duration of the industrial run and the in silico experiment duration of the metamodel. For this, the expected duration of the industrial runs needs to be known. In case of variable durations of the industrial runs, the duration can be added as an additional factor to the design space of the in silico experiments, thus extending the metamodel by one parameter.

More generally, the concept of the regression metamodel based on in silico process simulations offers the ability that further parameters of the underlying mechanistic-empirical process dynamics model can be integrated into the metamodel without experimental effort on site. Possibly, besides the reference filter cake resistance some more, currently fixed, model parameters could be changed to an adaptable form to improve the predictive power of the process dynamics model during the industrial runs. In addition, a higher flexibility of the modeled operation conditions, such as a permeate flow rate ramp can be added to open new possibilities for productivity optimization. Last not least, the optimization objective is imagined to be extended from only productivity to other objective variables, for instance product yield, energy demand of the membrane drive or membrane cleaning effort.

5. Conclusion

In this study, a digital twin for an industrial DCF process was developed and applied successfully for increasing the campaign productivity. The digital twin was constructed around a process dynamics model that we based on mechanistic filtration knowledge and fine-tuned with empirical correlations derived from historical process data. A regression metamodel was interposed between the basal process dynamics model and the digital twin framework, because we found that the metamodel approach is superior to several calculations of the mechanistic-empirical process dynamics model in terms of computational speed and, together with the step-by-step simulation of the process dynamics model, a setpoint recommendation is obtained within seconds. By comparing the TMP live data with the simulation result, the process dynamics model is continuously validated during the production run. If a non-replicated process effect occurs, the digital twin adjusts a model parameter already included in the meta model. By this, the meta model approach holds even in the case of changes in the filtration properties of the feed, and allows the robust prediction of an optimal permeate flow rate setpoint.

Given a mechanistic process model, the presented method can be used to optimize a (bio)chemical process without having to perform repeated computation-intensive calculations during the operation. The expandability of the metamodel offers the prospect of integrating further process parameters into the optimization or optimizing further target variables.