Manon Granger-Delacroix

Nadine Leconte

Cyril Grassin

Françoise Le Goff

Fabienne Garnier-Lambrouin

Marieke Van Audenhaege

Geneviève Gésan-Guiziou
email: genevieve.gesan-guiziou@inrae.fr

Skimmed milk microfiltration in diafiltration mode: Impact of solvent nature and concentration factor on spiral-wound membrane performance operated at low temperature

Keywords:

Cold skimmed milk microfiltration using polymeric spiral-wound membranes has a growing interest to separate serum proteins from casein micelles. Diafiltration using reverse osmosis water (ROW-DF) is needed to increase the transmission of serum proteins. The objective of this study was to investigate the performance of skimmed milk microfiltration in ROW-DF mode when operating with single industrially sized standard SW modules at low temperature (~12 • C). As expected during ROW-DF, the permeation flux increased with the rate of diafiltration, due to both the decrease of permeate viscosity and decrease of ionic strength. It is shown for the first time, that the transmission of protein (mainly β-Lactoglobulin, β-LG, major serum protein of bovine milk) presents different scenarios (increase, decrease and 'parabolic' shape) as function of diafiltration rate and operating conditions. These performances are quite different from the ones obtained when operating diafiltration with permeate from milk ultrafiltration, PUF-DF, and can be explained by the modifications induced by ROW on the properties of the fluid (feed) and on the properties of the concentrated layers of casein micelles accumulated at the membrane. The quantity of β-LG which can be released into the permeate was linked to the concentration factor applied before ROW-DF and the ability of deposit to swell was linked to the transmembrane pressure applied as critical conditions could lead to a cohesive deposit unable to swell. If the increase of repulsive electrostatic interactions (induced by the increase of diafiltration rate) between casein micelles and β-LG supported the release of β-LG entrapped in the deposit, it impaired the crossing of β-LG remaining in the retentate. The result was a progressive decrease of β-LG transmission with the increasing diafiltration rate. In summary, transmission of β-LG was maximized for low value of transmembrane pressure (0.7 bar), a high value of concentration factor (3.2) and an intermediate diafiltration rate (2.1).

Introduction

Fractionation of skimmed milk proteins using crossflow microfiltration (MF) is widely used in the dairy industry [START_REF] Fauquant | Microfiltration du lait sur membrane minerale[END_REF][START_REF] Saboya | Current developments of microfiltration technology in the dairy industry[END_REF][START_REF] Brans | Membrane fractionation of milk: state of the art and challenges[END_REF]. The increasing number of MF equipment is explained by a large panel of applications using retentate and permeate fractions which have advantageous composition and functionality [START_REF] Singh | Functional properties of milk proteins[END_REF]. So far, separation of casein micelles and serum proteins (SP) by MF is commonly operated at a temperature above 50 • C using ceramic membranes with a 0.1 µm mean pore diameter. However, the number of MF plants equipped with spiralwound (SW) membranes is aimed to increase in dairy industry due to the low prices of polymeric membranes (800 kDa -0.1 µm) and their high compactness. These membranes are expected to be mainly conducted at low temperature (~10 • C) to limit bacterial growth in the casein micelles retentate, as it is the case in skimmed milk ultrafiltration performed with SW membranes.

To date, MF with SW membranes was reported to have poor performances (permeation flux (Jp) and serum protein transmission) compared to the ones obtained with ceramic membranes [START_REF] Zulewska | Efficiency of serum protein removal from skim milk with ceramic and polymeric membranes at 50 degrees C[END_REF]. Lower performances can be partly explained by the higher value of transmembrane pressure (TMP) applied in MF with SW membranes compared to ceramic ones [START_REF] Hartinger | Milk protein fractionation by means of spiral-wound microfiltration membranes: effect of the pressure adjustment mode and temperature on flux and protein permeation[END_REF]. The low temperature used for MF with SW membranes could also be suspected to be responsible for a lower serum protein transmission but results obtained by Hartinger et al., (2019a) at 10 • C and 50 • C are not clear about the impact of temperature because TMP also varied during their experiments.

To increase serum protein transmission during MF operating with SW membranes, filtration can be carried out in diafiltration mode (DF). DF is a washing procedure where a solvent -typically reverse osmosis water (ROW) or permeate from milk ultrafiltration (PUF) -is added to the retentate, either in batch and intermitted mode or in continuous mode, to drag residual serum proteins through the membrane. In addition to the transmission of serum proteins, DF with ROW eliminate the minerals and the lactose from retentate whereas DF with PUF keep their concentrations constant. Lactose is known to be highly responsible for the viscosity of the soluble phase and minerals are partly responsible for the structure of the casein micelles, which accumulated at the membrane surface [START_REF] Reitmaier | Impact of different aqueous phases on casein micelles: kinetics of physicochemical changes under variation of water hardness and diafiltration conditions[END_REF]. Then, different performances of filtration are expected when operating MF in DF mode with PUF or ROW.

Microfiltration carried out in diafiltration mode with PUF has already been shown to be efficient, both with ceramic and SW membranes. Heidebrecht and Kulozik [START_REF] Heidebrecht | Fractionation of casein micelles and minor proteins by microfiltration in diafiltration mode. Study of the transmission and yield of the immunoglobulins IgG IgA and IgM[END_REF] studied the transmission of serum proteins during continuous DF of skimmed milk retentate (concentration factor, CF = 2) with PUF, using Isoflux ceramic membranes (Tami industries, 0.14 µm, 50 • C) at various TMP (from 0.6 to 3.0 bar). They showed that, the transmissions of native serum proteins were stable and constant up to a DF rate of 7 with no dependence on the TMPs (except for the lower TMP of 0.6 bar). The transmissions of native serum proteins were satisfactory with an order of transmission depending on the size of serum proteins and electrostatic interactions between the individual proteins and the casein micelles accumulated at the membrane surface (the transmission of α-lactalbumin (55 %) > β-lactoglobulin (50 %) > IgG (47 %) > other minor proteins). However, during their experiments, the transmissions of the overall serum protein were impaired by the long-term process exposure at 50 • C, which was responsible for a denaturation of serum proteins. Moreover, this study does not bring data on permeation flux obtained during the experiments performed.

Hartinger and Kulozik [START_REF] Hartinger | Milk protein fractionation by spiral-wound microfiltration membranes in diafiltration mode -Influence of feed protein concentration and composition on the filtration performance[END_REF] studied the performance of skimmed milk microfiltration in continuous PUF diafiltration mode performed with SW membrane (Synder V0.1) at low temperature (10 • C). They showed that regardless of the CF applied (1 -4) the permeation flux was stable during DF up to the diafiltration rate of 7. During DF performed at CF 1 and TMP 0.7 bar, the transmission of β-lactoglobulin was stable up to DF rate of 7. According to the results obtained, it could be expected that protein transmission is dependent of the CF value chosen for the DF step (from 1 to 4). However, as described by Heidebrecht and Kulozik, [START_REF] Heidebrecht | Fractionation of casein micelles and minor proteins by microfiltration in diafiltration mode. Study of the transmission and yield of the immunoglobulins IgG IgA and IgM[END_REF], Hartinger and Kulozik, [START_REF] Hartinger | Milk protein fractionation by spiral-wound microfiltration membranes in diafiltration mode -Influence of feed protein concentration and composition on the filtration performance[END_REF] operated the MF and the DF at constant axial pressure drop of 1.0 bar⋅m -1 meaning that the cross-flow velocity was not similar between the different concentration factors. In these conditions, it is not possible to conclude about the role of operating conditions (concentration factor and crossflow velocity) on protein transmission during DF operated with PUF.

Diafiltration mode with water is the most used mode in dairy industry. It has already been investigated with ceramic membranes [START_REF] Tremblay-Marchand | A process efficiency assessment of serum protein removal from milk using ceramic graded permeability microfiltration membrane[END_REF] and SW membranes at high temperature (~50 • C) [19] and ambiant temperature (~24 • C) [START_REF] Marella | Optimization of spiral-wound microfiltration process parameters for the production of micellar casein concentrate[END_REF]. These studies carried out batch and intermittent DF (consisting in succession of dilutions and concentrations) with ROW during MF at 50 • C and CF 3 [START_REF] Tremblay-Marchand | A process efficiency assessment of serum protein removal from milk using ceramic graded permeability microfiltration membrane[END_REF]19] or at 24 • C and CF 4 [START_REF] Marella | Optimization of spiral-wound microfiltration process parameters for the production of micellar casein concentrate[END_REF]. They confirmed the increase of permeation flux during DF with ROW as previously observed by Ng et al (2008) when operating ROW DF with ultrafiltration polymeric membranes. They also showed that the purity of retentate (casein/true proteins) and the transmission of serum proteins were higher when the rates of DF increased in the range from 0.7 to 2.0. However, it is difficult from these studies to understand the role of dilution and resulting physico-chemical modifications on structural changes of casein micelles deposit layer and on transmission of serum protein as DF was performed in intermittent mode resulting in valleys and peaks in the flux and protein removal graphs.

Reitmaier et al. [START_REF] Reitmaier | Impact of different aqueous phases on casein micelles: kinetics of physicochemical changes under variation of water hardness and diafiltration conditions[END_REF][START_REF] Reitmaier | Effect of changes in ionic composition induced by different diafiltration media on deposited layer properties and separation efficiency in milk protein fractionation by microfiltration[END_REF] investigated the impact of different aqueous phases (including PUF and ROW) on casein micelles characteristics, deposited casein layers properties and separation properties in milk protein fractionation by MF. They showed that ROW induces a significant increase of deposited casein layer hydration, thereby increasing flux levels and separation process. However, they operated at 50 • C with a single ceramic membrane and the transferability of the observed effects to common application of SW membranes at ~ 10 • C has not been demonstrated so far. This information is crucial for dairy industry, as SW membranes at low temperature have a large potential.

The objective of this study was to investigate the performance of skimmed milk microfiltration in ROW-DF mode when operating with single industrially sized standard SW modules at low temperature (~12 • C). In order to identify optimal operating conditions and to maximize serum protein transmission from the retentate, the permeation flux and transmissions of the two major individual serum proteins (α-lactalbumin and β-lactoglobulin) were measured up to a DF rate higher than 5.0. DF was performed in continuous mode, on skimmed milk concentrated at CF ranging from 2.5 to 3.5 and at two TMP values. The role of structural changes in the deposit layer was examined by comparing results with those obtained with a diafiltration mode operating with PUF in similar operating conditions.

Materials and methods

Fluids

Thermized skimmed milk was provided by Entremont Alliance (Bretagne, France) (Table 1) and stored at 4 the aqueous phase of skimmed milk at 12 • C, the PUFs showed a composition with lower soluble calcium and phosphate contents (Table 1). This difference is attributed to the lower solubility of phosphate and calcium at 45 • C (temperature of PUF preparation) compared to 12 • C resulting in lower release of calcium and phosphate from the casein micelles. Lactose was slightly higher for PUFs than for skimmed milk.

Microfiltration set-up and experimental protocol

Pilot scale MF was carried out using single industrially sized standard polymeric 800 kDa SW membrane (FR 3A 6338, spacer of 46 mil, PVDF, total membrane area of 15.9 m 2 , Synder, California, USA). All the filtration experiments were performed with the same cleaned membrane. Membrane was first cleaned with sodium hydroxide (alkaline solution, NaOH 30 %, Quaron, France) at pH 11 and 45 • C, with a retentate recirculating flow rate of 18 m 3 ⋅h -1 and two successive steps: the first 10 min were performed without permeation and the next 10 min with permeation. Another similar sequence was then performed with nitric acid (acid solution, HNO 3 ; 58 % purity, Quaron, France) at pH 2.5 and 45 • C with a recirculating flow rate of 18 m 3 ⋅h -1 . Water used to rinse and clean the filtration pilot before and after MF was water from network filtered sequentially on 5.0, 1.0 and 0.2 μm cartridges. The hydraulic resistance, calculated according to Darcy's law, was recovered after cleaning: Rm = (3.1 ± 0.5) × 10 11 m -1 .

Before MF, water initially present into the retentate compartment was gently flushed with skimmed milk (three times the volume of the retentate compartment).

All the experiments were performed at constant temperature (12.0 ± 0.3 • C), mean retentate pressure (2.0 ± 0.1 bar), and retentate recirculating flow rate (18.0 ± 0.7 m 3 ⋅h -1 ). They were carried out at CF ranging from 2.5 to 3.5 and at two constant TMP, 0.7 ± 0.1 and 0.9 ± 0.1 bar. TMP = 0.7 ± 0.1 bar corresponded to the lowest pressure applicable in the range of CF studied to avoid retro-filtration (TMP > ½ Retentate Pressure drop). TMP = 0.9 ± 0.1 bar was a value closer to critical conditions of filtration.

Each MF experiment was divided into three phases: a concentration phase, a stabilization phase, and a MF phase operated in DF mode. Concentration phase of skimmed milks was realized by concentrating the retentate either in the feed tank (by returning retentate into the feed tank while permeate is extracted) or in the retentate loop. During the stabilization phase (15 to 90 min), MF was performed at a constant CF and retentate was continuously extracted at a flowrate fixed by the targeted CF. During the diafiltration phase, known volume of retentate was continuously diafiltrated with ROW (ROW-DF) or PUF (PUF-DF). The chosen diafiltration solvent was hold at 12 • C for 30 min before the experiment. DF was carried out at constant flowrate, meaning that the DF flowrate was maintained equal to the permeation flowrate.

During the three phases of the experiment, Jp and SP transmission were evaluated as performance indicators.

Analyses

Skimmed milk, retentate and permeate were collected for analyses. Fat content was measured by the Gerber method [21]. Dry matter was obtained after desiccation of the sample at 105 • C for 7 h. Ashes were measured after combustion of 5 g of sample at 550 • C for 5 h. Total and soluble calcium were determined by atomic absorption spectrophotometer (AA300, Varian France) after dilution of the samples in a solution containing 10 % (v/w) of 6 g⋅L -1 lanthanum chloride and 10 % (v/w) of N/50-hydrochloric acid [START_REF] Brule | Etude de la teneur en éléments minéraux des produits obtenus lors de l'ultrafiltration du lait sur membrane[END_REF]. The phosphate content was measured by chromatography method (Dionex DX 500 HPLC, Thermo Fisher Scientific, Les Ulis, France) as described by Gaucheron et al. [START_REF] Gaucheron | Determination of anions of milk by ion chromatography[END_REF]. Total cations contents were determined from ashes solubilized in 1 Nhydrochloric acid. Soluble ions were extracted from samples by ultrafiltration on 10 kDa membranes (Vivaspin VS2002, Sartorius, Aubagne, France). Conductivity was determined at 12 • C (samples were maintained 20 min at 12 • C before measurement) with a waterproof four-ring ATC conductivity meter (HI 9033, HANNA instruments, measurement error 1 %).

Lactose content was determined by high-performance liquid chromatography (Dionex HPLC) according to a protocol adapted from Thierry et al. [START_REF] Thierry | Conversion of l-Leucine to Isovaleric Acid by Propionibacterium freudenreichii TL 34 and ITGP23[END_REF]. Briefly, lactose was extracted from samples by ultrafiltration (12 • C, 1 800 g, 30 min) on 10 kDa membranes (Vivaspin VS2002, Sartorius) before dilution with 0.01 N-sulphuric acid. HPLC were realized using a column packed with an ion-exchange resin (Aminex A-6, Bio-Rad, France).

Total nitrogen (TN), non-casein nitrogen (NCN) and non-protein nitrogen (NPN) were determined by the Kjeldahl method according to the standard [START_REF]Milk and milk products -Determination of nitrogen content -Part 1: Kjeldahl principle and crude protein calculation[END_REF], (International Dairy Federation, 2014). For protein content calculation, 6.38 was used as conversion factor. For protein content calculation in NCN and NPN filtrates, a correction factor was calculated to consider the weight of precipitate. TN was determined with an experimental error of 1 % and NCN and NPN with an experimental error of 5 % (Experimental errors were calculated as the ratio between standard deviation and the mean value).

Concentrations of α-lactalbumin (α-LA) and β-lactoglobulin (β-LG) were determined in skimmed milks, retentates and permeates by reversed phase high-performance liquid chromatography (RP-HPLC) according to a protocol adapted from Resmini et al. (1989). RP-HPLC was performed with a PLRP-S column (PL1912-3801, 300 Å, 8 μm, 150 mm × 2.1 mm, Agilent Technologies) on an Alliance 2695 (Waters, Saint-Quentin-en-Yvelines, France). Contents of α-LA and β-LG in milks and retentates were determined from the NCN filtrates. In both skimmed milks and retentates, measured values correspond to content of α-LA and β-LG present in the aqueous phase either in native or partially denatured forms. In permeates, contents of α-LA and β-LG were measured directly on samples and measured values correspond to the total content of α-LA and β-LG. Contents of α-LA and β-LG were calculated with an experimental error of 3 %.

All the analyses were realized in duplicate, some have been done in triplicate or quadruplicate.

Viscosity of permeates was measured at 12 • C by shear stress with a viscometer (LS400, LAMY RHEOLOGY Instruments, Champagne au Mont d'Or, France). Samples were maintained at 12 • C during 5 min before measurement.

Calculations

The rate of diafiltration (Tx), was calculated as:

Tx = V d V R0 = V p V R0 (1)
with V d , V R0 and V p the volumes (L) of solvent of diafiltration, of retentate to be diafiltered and of permeate extracted during the microfiltration in DF mode.

The transmission (Tr) of β-LG or α-LA in the aqueous phase was calculated according to:

Tr = ( C p C r ) × 100 (2)
with C p and C r the concentrations of β-LG or α-LA measured in permeate and retentate respectively at different sampling times (g⋅L -1 ). During diafiltration and in case of total transmission of component (e.g. case of lactose, soluble ions), the retentate concentration of the component can be calculated from the mass balance equation defined as [START_REF] Foley | Minimisation of process time in ultrafiltration and continuous diafiltration: the effects of incomplete macrosolute rejection[END_REF][START_REF] Oueslati | A simple methodology for predicting the performances of hyaluronic acid purification by diafiltration[END_REF]:

C r (Tx) = C r0 e -Tx (3)
with C r0 (g⋅kg -1 ) the concentration of the component in the retentate to be diafiltered and C r its concentration in the retentate at Tx.

In skimmed milk, calcium and phosphate ions are distributed between the solvent and the colloidal phase (casein micelle). The concentration (mg⋅L -1 ) of a specific ion i in the colloidal phase, C i;colloidal was calculated as suggested by Ferrer et al. [START_REF] Ferrer | Does ultrafiltration have a lasting effect on the physico-chemical properties of the casein micelles?[END_REF]:

C i;colloidal = C i;total -C i;soluble × VF (4)
with C i;total , the total concentration of the ion i; C i;soluble its concentration in the soluble phase; and VF, the volume factor that considers the effects of the dilution, the non solvent water and the volume occupied by the casein micelles. It was calculated as suggested by Ferrer et al. [START_REF] Ferrer | Does ultrafiltration have a lasting effect on the physico-chemical properties of the casein micelles?[END_REF]:

VF = (1 - ( 1.01 × C fat + 0.03 × C CN + 0.2 × C TP )) × ( ρ i ρ s ) (5)
where C fat , C CN and C TP represent respectively the content of fat, casein (CN = TN-NCN) and true proteins (TP = TN-NPN) in g⋅kg -1 . ρ i and ρ s are the densities of the sample and the soluble phase, respectively. They are estimated from the content of casein and lactose of the sample, based on a correlation established with data from Sodiaal (confidential). The Concentration Factor, CF is defined by.

CF = ( C CN-ret C CN-Milk ) (6)
where C CN-ret and C CN-Milk represent respectively the content of casein in retentate and initial milk.

Results

Concentration and stabilization phases

As expected, during the concentration phase, Jp decreased with the increase of CF. Between CF 1.0 and 3.2 and for a TMP of 0.7 bar, Jp dropped by about 50 %, from 14.2 to 7.4 L⋅h -1 ⋅m -2 . These results are consistent with the work of Hartinger and Kulozik, [START_REF] Hartinger | Milk protein fractionation by spiral-wound microfiltration membranes in diafiltration mode -Influence of feed protein concentration and composition on the filtration performance[END_REF] who observed a drop of about 60 % of Jp when CF raised from 1.0 to 3.0 during MF carried out with SW membrane in similar conditions (T = 10 • C, TMP = 0.7 bar). At CF 2.8, Jp was equal to 7.5 and 10.7 L⋅h -1 ⋅m -2 with a TMP of 0.7 and 0.9 bar respectively, confirming that the limiting flux was not reached at TMP = 0.7 bar and CF = 2.8.

Transmission of β-LG (Tr β-LG ) and α-LA (Tr α-LA ) remained constant with the increase of CF, regardless of the mode of concentration (concentration in feed tank or in the retentate loop). Despite a high variability of Tr α-LA from one experiment to another (Tr α-LA varied from 57 ± 1 % to 78 ± 2 % for a TMP of 0.7 bar) Tr α-LA was always higher than Tr β-LG , in accordance to the size of the proteins (α-LA: 14.2 kg⋅mol -1 and β-LG (dimer): 36 kg⋅mol -1 ). The lower the TMP, the higher the transmission of serum proteins: Tr β-LG was 36 ± 2 % (n s = 14) and 30 ± 1 % (n s = 4) at TMP 0.7 and 0.9 bar respectively, and Tr α-LA was 69 ± 9 % and 58 ± 3 % at TMP 0.7 and 0.9 bar respectively (n s is the number of samples considered for the calculation of the transmission. The samples used for the calculation of serum protein transmission are samples from the "concentration and stabilization phases" of experiments performed at a given applied TMP. The transmission values obtained at different CFs and time points were averaged because the transmission values of the serum proteins were constant at a given applied TMP over the CF range studied). These results are consistent with the works of Gésan-Guiziou et al. [START_REF] Gésan-Guiziou | Critical stability conditions in crossflow microfiltration of skimmed milk: transition to irreversible deposition[END_REF] and Hartinger et al. [START_REF] Hartinger | Investigation on the spatial filtration performance in spiral-wound membranes -Influence and length-dependent adjustment of the transmembrane pressure[END_REF].

When CF increased and regardless of the TMP applied, the contents of total calcium and phosphate increased in the retentate compartment. At the same time, the contents of soluble calcium and phosphate remained constant (Results not shown) suggesting that increase in the total content was induced by the colloidal calcium and phosphate. During the concentration phase, the proportions of colloidal calcium and phosphate on casein content expressed as C Ca;colloidal /CN and C phosphate; colloidal /CN ratios (Eq. ( 4)) were steady: 0.87 ± 0.05 and 0.36 ± 0.02 g⋅mol -1 respectively. These results suggested that calcium phosphate present in casein micelles was not solubilized in the aqueous phase during concentration and that only the colloidal fraction of minerals associated to casein micelles was concentrated. These results are consistent with Ferrer et al. [START_REF] Ferrer | Does ultrafiltration have a lasting effect on the physico-chemical properties of the casein micelles?[END_REF] who did not observe any solubilization of colloidal calcium phosphate during ultrafiltration of casein micelles at 40 • C and CF lower than 5.

During the stabilization phase, the permeation flux and the transmission of α-LA and β-LG were stable during 15 to 90 min, regardless of TMP and CF values used.

Diafiltration phase

4.2.1. Diafiltration with PUF (PUF-DF) at TMP 0.7 bar Regardless of the CF, Jp remained stable during the PUF-DF carried out at TMP = 0.7 bar: Jp = 9.6 ± 0.2 and 7.9 ± 0.3 L⋅h -1 ⋅m -2 at CF 2.5 and 2.8 respectively (Fig. 1a).

During the PUF-DF, and regardless of CF, Tr β-LG and Tr α-LA , calculated from experimental retentate and permeate concentrations, decreased (Fig. 2a &b). Between Tx = 0.0 and 2.0, Tr β-LG decreased from 35 to 23 % for CF = 2.5 and from 34 to 31 % for CF = 2.8 while Tr α-LA dropped from 75 to 51 % and from 57 to 39 % respectively.

During the PUF-DF, the permeate lactose concentration slightly increased (Results not shown) due to the difference of lactose content between the PUFs and the skimmed milk. Conductivity of retentate fractions remained stable at 3110 ± 78 µS⋅cm -1 (Fig. 1 b) and mineral contents of the retentate fractions remained stable. These results suggest that the slight differences of lactose and mineral contents observed between the soluble fractions of skimmed milk and PUFs (Table 1) have no consequence on the ionic environment of skimmed milk during PUF-DF nor on Jp values.

Diafiltration with ROW (ROW-DF) at TMP 0.7 and 0.9 bar

During ROW-DF, Jp increased progressively with the increase of Tx (Fig. 1c). The increase of Jp (ΔJp) was of the same order of magnitude for the experiments performed with the low CF values (from 2.6 to 3.2) (Table 2): during diafiltration operated between Tx = 0.0 and 2.0, Jp increased by 44 and 36 % at TMP 0.7 bar and CF 2.6 and 3.2 respectively and by 42 % during DF operated at TMP 0.9 bar and CF 3.0. However, ΔJp was much lower during DF carried out with the high values of CF (3.5) and TMP (0.9 bar): in these conditions, Jp increased by only 15 % between Tx = 0.0 and 2.0.

For higher rate of diafiltration (Tx > 2.0), ΔJp levelled off. For instance, Jp increased by only 28 % (from 42 to 54 L.h -1 ⋅m -2 ) between Tx = 2.0 and 4.0 with TMP = 0.9 bar and CF 3.0.

During ROW-DF at TMP of 0.7 bar, Tr β-LG followed a "parabolic shape", which was more pronounced with higher CF (3.2compared to 2.6) (Fig. 2c): for ROW-DF operated at CF 3.2, Tr β-LG increased from to 50 % when Tx ranged from 0.0 to 2.1 and then decreased to 26 % at Tx = 5.7. For ROW-DF operated at CF 2.6, Tr β-LG increased from 33 to % when Tx ranged from 0.0 to 0.5 and then decreased to 31 % at Tx = 2.0. Simultaneously, Tr α-LA regularly decreased during the ROW-DF at CF 3.2 (from 77 to 66 and 38 % at Tx = 0.0, 2.1 and 5.7 respectively) but remained stable during ROW-DF at CF 2.6 (64 ± 2 %) (Fig. 2d).

During ROW-DF at TMP of 0.9 bar, Tr β-LG slightly increased regardless of CF (from 30 ± 2 % to 37 ± 3 % between Tx = 0.0 and 4.0) (Fig. 2e). Simultaneously, Tr α-LA decreased from 61 ± 5 % to 37 ± 3 % regardless of CF (Fig. 2f).

During ROW-DF, viscosity of the permeate dropped with the increase of Tx (Fig. 3b). The viscosity of permeate decreased from about 1.45-1.50 mPa⋅s at Tx = 0.0 to 1.25 mPa⋅s after Tx = 4.0. The permeate concentration of lactose decreased progressively when Tx increased (Fig. 3a): it dropped from 48 ± 2 g.L -1 at Tx = 0.0 to 6 and 1 g⋅L -1 at Tx = 2.0 and 4.0, respectively. The drop of permeate concentration of lactose during ROW-DF was directly induced by the dilution of the retentate lactose concentration with ROW as the lactose (342 g⋅mol -1 ) is totally transmitted through the MF membrane (prediction according to Eq. ( 3)).

Conductivity of the retentate decreased when Tx of ROW-DF increased (Fig. 1d). Regardless of TMP and CF applied the decline of the conductivity was similar, from 3003 ± 140 µS.cm -1 to 449 µS.cm -1 for Tx = 0.0 and 4.0, respectively.

Simultaneously, total and soluble contents of phosphate in the retentate dropped when Tx increased (Table 3). However, the drop of total phosphate was higher than the drop of soluble phosphate. Between Tx = 0.0 and Tx = 2.0, total phosphate content dropped by 806 mg⋅L -1 whereas soluble phosphate content only dropped by 553 mg⋅L -1 at TMP = 0.7 bar and CF 2.6. Similar observations are done at TMP = 0.7 bar and CF = 3.2.

These results indicate a release of colloidal phosphate in the soluble phase induced by a partial solubilization of phosphate during ROW-DF: the ratio of C phosphate-colloidal /CN calculated with Eq. ( 4) confirmed this hypothesis as it dropped from 0.37 to 0.28 mol⋅g -1 between Tx = 0.0 and 2.0 when DF is carried out at TMP = 0.7 bar and CF = 2.6.

Discussion

Realising dairy protein fractionation with SW membranes requires improving the performance of skimmed milk microfiltration equipped with these membranes. A solution is to carry out MF in diafiltration mode with reverse osmosis water (ROW-DF). This study provides elements to define the optimal operating conditions to conduct cold MF in ROW-DF mode with SW membrane. It gives some guidelines for the industry to maximize the removal of serum proteins and minimize the volume of ROW to be used. This work shows that the performances (permeation flux and individual serum protein transmission) of skimmed milk MF operated in ROW-DF mode are quite different from the ones obtained when operating PUF-DF. The modifications observed can be explained by the modifications induced by ROW on the properties of the fluid (feed) including the characteristics of casein micelles, and on the properties of the concentrated layers of casein micelles accumulated at the membrane.

Permeation flux with ROW-DF

The increase of Jp with the increase of Tx during ROW-DF is in general agreement with previous studies. To our knowledge, this is the first time that an increase of Jp was observed when MF is operated in ROW-DF mode at 12 • C, but this result is very similar to the ones previously obtained by others when operating ROW-DF at ~ 50 • C [START_REF] Tremblay-Marchand | A process efficiency assessment of serum protein removal from milk using ceramic graded permeability microfiltration membrane[END_REF]19,[START_REF] Marella | Optimization of spiral-wound microfiltration process parameters for the production of micellar casein concentrate[END_REF].

Role of permeate viscosity

The increase of Jp observed during ROW-DF can partially be explained by the decrease of permeate viscosity. When the viscosity of the permeate is unchanged (case of PUF-DF), Jp is constant. In contrast, when the viscosity of the permeate decreases (case of ROW-DF), Jp increased. It can be assumed that the decrease in viscosity of the permeate is mainly attributed to the decrease in lactose content as it is largely accepted that the concentration of serum proteins has little impact on viscosity as their concentration is about 8 times lower than that of lactose. Taking this into account, the viscosity of permeate can be estimated as a function of lactose content (and then as a function of Tx) from our experimental data. An estimated permeation flux corrected by the viscosity of the permeate for each Tx (Jp corr ) can then be calculated (Fig. 4a &b): [START_REF] Ferrer | Does ultrafiltration have a lasting effect on the physico-chemical properties of the casein micelles?[END_REF] with Jp the permeation flux (L⋅h -1 ⋅m -2 ) measured during the DF and µ p (mPa⋅s) the viscosity of the permeate at a given Tx. Jp corr during PUF-DF was calculated by considering the slight evolution of lactose content of permeate resulting from the slight differences in lactose content between skimmed milk and PUFs used for the PUF-DF. Finally, the percentage of increase of Jp induced by the decrease of viscosity (ε) can be calculated as follows: 4b). Exception has to be made for the ROW-DF on high CF (3.5) and high TMP (0.9 bar) (open triangle -Fig. 4b) for which 100 % of the Jp increase is explained by the reduction of the permeate viscosity. This is an interesting result which indicates that the permeability properties of casein micelles fouling formed in high TMP (0.9 bar) and CF (3.5) are not modified during ROW-DF whereas properties of accumulated layers formed at lower TMP or CF changed.

Jp corr (T x ) = Jp(T x ) × μ p (T x ) μ p (Tx=0)

Role of ionic strength

As already mentioned in the literature (at 50 • C) (Ng et al. 2008, [25], the increase of Jp during ROW-DF may be attributed to the decrease in ionic strength of retentate (conductivity as indicator) that may affect the casein micelles properties and the features of accumulated layers at the membrane surface.

The decrease of ionic strength leads to modification of the internal properties of casein micelles [START_REF] Lazzaro | Comprendre les comportements des micelles de caséines dans des environnements variés, de leur équilibre minéral à leurs propriétés colloïdales et fonctionnelles : émulsion et coagulation présure[END_REF]. In this work, the « surplus » of phosphate recovered in permeate during the ROW-DF was estimated as the difference between the content of phosphate measured in the permeate at Tx and the expected phosphate content in the permeate at Tx calculated from the soluble phosphate content in the retentate (at Tx) and the mass balance equation (Eq. ( 3)). Such a calculation shows a solubilization of phosphate from the casein micelle to the soluble phase (Fig. 5) that is similar for all investigated operating conditions of ROW-

Table 2

Increase of permeation flux (ΔJp) in % during diafiltration with reverse osmosis water (ROW-DF) at a transmembrane pressure value (TMP) of 0.7 or 0.9 bar and a concentration factor (CF) from 2.6 to 3.5. Abbreviations. Tx, rate of diafiltration; Jp, permeation flux.

DF. This result is in agreement with Li and Corredig [START_REF] Li | Calcium release from milk concentrated by ultrafiltration and diafiltration[END_REF] and Alexander et al. [START_REF] Alexander | Changes in the calcium cluster distribution of ultrafiltered and diafiltered fresh skim milk as observed by Small Angle Neutron Scattering[END_REF] who showed partial solubilization of micellar calcium during skimmed milk UF (CF = 2.5) operated in ROW-DF at Tx 2.5 to 5.0. This result was supported by Ferrer et al. [START_REF] Ferrer | Changes in the physico-chemical properties of casein micelles during ultrafiltration combined with diafiltration[END_REF] who additionally showed that the structure (e.g. higher turbidity due to lower refractive index of the casein micelles) and the physical properties of casein micelle (e.g. lower susceptibility to aggregation, lower storage modulus of rennet induced gels) were modified by the solubilization of calcium and individual caseins at high rate of DF with ROW (>50). Lazzaro [START_REF] Lazzaro | Comprendre les comportements des micelles de caséines dans des environnements variés, de leur équilibre minéral à leurs propriétés colloïdales et fonctionnelles : émulsion et coagulation présure[END_REF] explained also that solubilization of colloidal minerals leads to a swelling of the casein micelles that might result in an increase of Jp. However according to our results, it is unlikely that the solubilization of colloidal minerals affects Jp. Indeed, the solubilization of colloidal minerals was similar regardless of operating conditions used (Fig. 5) but at TMP 0.9 bar and CF 3.5, Jp corr (Jp corrected by the permeate viscosity) remained constant during ROW-DF despite the increase of Tx and solubilization of colloidal minerals (Fig. 6).

The decrease in ionic strength of retentate may be more likely responsible for change in the properties of the protein accumulated layers, which themselves dependent on the operating conditions.

The decrease of ionic strength may result in an increase of electrostatic repulsions between charged proteins (casein micelles and entrapped serum proteins) and swelling of accumulated protein layers [START_REF] Jimenez-Lopez | Ionic strength dependence of skimmed milk microfiltration: Relations between filtration performance, deposit layer characteristics and colloidal properties of casein micelles[END_REF]. As a result, the resistance of the accumulated layers is reduced, and Abbreviation. Tx, rate of diafiltration. the permeation flux increases. The same conclusion was found by Ng et al. [START_REF] Ng | Influence of diafiltration on flux decline during skim milk ultrafiltration[END_REF] during ultrafiltration of skimmed milk in ROW-DF mode. However, our work shows that the swelling of the accumulated layers depends of the operating conditions. Indeed, in high CF and TMP conditions (respectively 3.5 and 0.9 bar in our work), the permeability of protein accumulated layers was not modified during ROW-DF: the stable Jp corr in these conditions could be explained by the cohesive interactions of the deposit. Qu et al., [START_REF] Qu | On the cohesive properties of casein micelles in dense systems[END_REF] demonstrated that the casein micelles deposit can shift from a solution state to a gel state under specific conditions of CF and critical TMP. We assumed that for high values of CF and TMP cohesive interactions of the casein micelles are high and modification of the ionic environment does not modify the structure of the deposit during ROW-DF.

Not surprisingly when the repulsive electrostatic charges are unmodified (case of PUF-DF), the structure of the deposit is unchanged and Jp corr is unchanged. These results are convenient with the work of Hartinger and Kulozik [START_REF] Hartinger | Milk protein fractionation by spiral-wound microfiltration membranes in diafiltration mode -Influence of feed protein concentration and composition on the filtration performance[END_REF].

Transmission of serum proteins

Transmissions of serum proteins have never been investigated in skimmed milk MF operated in ROW-DF mode at low temperature. Interestingly our study shows that the transmission of serum protein (mainly transmission of the major serum protein, β-LG) highlights various scenarios (decrease, increase, 'parabolic shape') versus Tx, depending on operating conditions (CF, TMP). These results indicate that the assumption of constant protein transmission classically considered in diafiltration [START_REF] Foley | Minimisation of process time in ultrafiltration and continuous diafiltration: the effects of incomplete macrosolute rejection[END_REF][START_REF] Oueslati | A simple methodology for predicting the performances of hyaluronic acid purification by diafiltration[END_REF] is not relevant in our conditions to establish a mass balance and predict the transmission and concentration of proteins in the permeate. The experimental values of β-LG protein transmission highlight the fact that the properties of casein micelles accumulated layers that depend on operating conditions, may have a crucial role in entrapment and/or release of β-LG during ROW-DF. We could reasonably assume that β-LG accumulates into the casein micelles deposit and that the swelling of the deposit and change of physicochemical environment during ROW-DF leads to the progressive removal of this accumulated β-LG into the permeate.

In order to better understand the different situations encountered in the study (experiments performed at low (0.7 bar) or high (0.9 bar) TMP as well as low (2.5) and high (3.5) CF), we calculated the "surplus" of protein in the permeate during a variation of ΔTx (ΔTx = Tx n -Tx n-1; with n: sample number). The "surplus" of protein in the permeate during ΔTx was calculated by the difference between the quantity of protein recovered in the permeate during ΔTx and the quantity of β-LG released from the retentate (central portion of the pipe) during the same ΔTx. It was then expressed as: Surplus of protein in the permeate during ΔTx

= C pm (ΔTx) × Vp(ΔTx) -(C r (Tx n-1 ) -C r (Tx n )) × Vr(ΔTx) (9)
with C pm , the average concentration of β-LG recovered in the permeate during ΔTx. C pm was calculated from the instantaneous C p measured values using the trapeze method; Vp (ΔTx) and Vr (ΔTx), the volumes recovered in the permeate and diafiltered in the retentate during ΔTx. As Vp (ΔTx) = Vr (ΔTx), the surplus is expressed in g per Liter of diafiltered retentate.

Two extreme scenarios can occur:

-When "surplus" > 0: More proteins than have disappeared from the retentate (pipe) are recovered in the permeate. This indicates that proteins accumulated in the surface layers are released into the permeate. -When "surplus" < 0: Proteins that are removed from the retentate do not reappear in the permeate, indicating retention of proteins in the surface layers of the membrane.

Depending of the operating conditions, three cases can be distinguished during ROW-DF at Tx ranging from 0 to 2 (Fig. 7): i. When the MF is carried out with a high CF and a low TMP (case of MF in DF mode with CF 3.2 and TMP 0.7 bar) (Fig. 7b), the "surplus" of β-LG in the permeate is high. In that case, it can be assumed that the accumulation of β-LG into the deposit before DF is high (high CF) and that the accumulated layers are "noncohesive" (low TMP, low compression). The swelling of the deposit leads to the release of high quantity of β-LG into the permeate. ii. When the MF is carried out with a low CF and a low TMP (case of MF in DF mode at CF 2.6 and TMP 0.7 bar) (Fig. 7b), only few β-LG is released into the permeate during the very first rate of diafiltration. In that case, the accumulation of β-LG into the deposit before DF is low (low CF) and the swelling of the accumulated layers leads to a release of a small quantity of β-LG. iii. When the MF is carried out with a high CF and a high TMP (case of MF in DF mode at CF 3.5 and TMP 0.9 bar) (Fig. 7c), a release of β-LG is observed up to an intermediate value of the diafiltration rate, Tx = 2.0. In that case, it can be assumed that the accumulation of β-LG into the accumulated layers at the membrane surface before DF is high (high CF). In these conditions, the casein micelles layers did not swell, the deposit stayed cohesive and Jp corr was unchanged despite the solubilization of colloidal minerals. Due to the decrease of ionic strength, repulsive electrostatic interactions increase between casein micelles and β-LG (both negatively charged at neutral pH) leading to the expulsion of β-LG out of the deposit. At high Tx, more extensive rate of ROW-DF is assumed to increase the repulsive electrostatic interactions between casein micelles and β-LG, resulting in reduction of β-LG accumulation into the deposit. When all the β-LG accumulated into the deposit before the ROW-DF is released, transmission β-LG starts to decrease.

Contrary to ROW-DF, when the structure of the deposit is unchanged and the repulsive electrostatic charges between proteins (casein micelles and serum proteins) are not modified (case of PUF-DF), β-LG accumulated into the deposit ("surplus" of β-LG in the permeate < 0; Fig. 7a) and transmission of β-LG slightly decreased with the increase in Tx (Fig. 2a). This result is different from those of Heidebrecht and Kulozik [START_REF] Heidebrecht | Fractionation of casein micelles and minor proteins by microfiltration in diafiltration mode. Study of the transmission and yield of the immunoglobulins IgG IgA and IgM[END_REF] and Hartinger and Kulozik [START_REF] Hartinger | Milk protein fractionation by spiral-wound microfiltration membranes in diafiltration mode -Influence of feed protein concentration and composition on the filtration performance[END_REF] who observed constant transmission of serum proteins during PUF-DF. However, Heidebrecht and Kulozik [START_REF] Heidebrecht | Fractionation of casein micelles and minor proteins by microfiltration in diafiltration mode. Study of the transmission and yield of the immunoglobulins IgG IgA and IgM[END_REF] carried out MF at 50 • C instead of 12 • C in this study, and it was recently shown that such high temperature leads to low cohesive deposit compared to 12 • C [START_REF] Doudiès | Major role of voluminosity in the compressibility and sol-gel transition of casein micelle dispersions concentrated at 7 • C and 20[END_REF], that could explain a higher transmission of β-LG. Moreover, Hartinger and Kulozik [START_REF] Hartinger | Milk protein fractionation by spiral-wound microfiltration membranes in diafiltration mode -Influence of feed protein concentration and composition on the filtration performance[END_REF] only focus on the transmission of serum proteins (β-LG) for the PUF-DF operated at CF = 1. For higher CF (>1), Hartinger and Kulozik [START_REF] Hartinger | Milk protein fractionation by spiral-wound microfiltration membranes in diafiltration mode -Influence of feed protein concentration and composition on the filtration performance[END_REF] assumed that transmission of β-LG was constant but we showed in this study that transmissions of serum proteins during diafiltration depend on CF.

Remarkbly, the behavior of α-LA is quite different from the one of β-LG. No 'parabolic shape' of transmission versus Tx can be observed at low TMP (0.7 bar) and high CF (3.2) during ROW-DF. This observation may be linked to the physicochemical properties of α-LA. As it is a metalloprotein and as it conformational structure depends on the presence of a calcium moiety, it can be assumed that changes of ionic strength and solubilization of calcium phosphate during ROW-DF do not only act on swelling of accumulated protein layers but also on the net charge and/or conformation of α-LA. As a consequence, electrostatic interactions between casein micelles and α-LA during ROW-DF could be different for the ones between β-LG and casein micelles.

Our work highlighted that cold MF performances in ROW-DF mode are impacted by the properties of the fluid and the properties of the accumulated layers of casein micelles. Critical conditions leading to cohesive properties of the deposit ('gel state', in our case TMP 0.9 and CF 3.5) impair the performances of the diafiltration. On the contrary, if the protein accumulated layers are not densely packed, the decrease of ionic strength increases the repulsive electrostatic interactions of casein micelles within the deposit that swells. As a result, both permeation flux and transmission of β-LG increased. The decrease of ionic strength also increases the repulsive electrostatic interactions between casein micelles and β-LG resulting in the release of entrapped β-LG in the deposit. However, when high rates of ROW-DF are applied, the repulsive electrostatic interactions between casein micelles and β-LG impairs the transmission of the β-LG into the permeate.

Conclusion

Performances (Jp and Tr β-LG ) of cold MF using spiral-wound membrane can be improved by DF mode using reverse osmosis water, ROW. In this work the effects of operating conditions (concentration factor, CF and transmembrane pressure, TMP) were studied. Regardless of the TMP, Jp increased when the rate of ROW-DF increased. On the contrary, Jp was constant with the rate of DF when DF was realized with permeate of milk ultrafiltration, PUF. When MF in DF mode with ROW was carried out at low TMP (0.7 bar), Tr β-LG had a parabolic shape, all the more pronounced that the CF was high. When the ROW-DF was realized at high TMP value (0.9 bar), Tr β-LG slightly increased. On the contrary, when DF was run with PUF, Tr β-LG decreased with the rate of DF. This work shows that the modifications of MF performances during DF mode were attributed to the modification of the properties of the fluid, of the properties of casein micelles and of the properties of the accumulated layers of casein micelles.

About 40 % of Jp increases are attributed to the drop of permeate viscosity induced by the removal of lactose. The drop of conductivity that leads to the increase of repulsive electrostatic interactions between casein micelles into the deposit, favors the swelling of deposit (if the MF operating conditions do not induce cohesive deposit), and consequently an increase of Jp. During ROW-DF, the β-LG accumulated into the deposit formed during MF, is released during the first rate of DF. If DF is extended, the increase of repulsive interactions between casein micelles and β-LG penalized its crossing in permeate. In case of cohesive deposit, the decrease of ionic strength has no impact on Jp, transmission of β-LG slightly increases before stabilizing due to the increase in repulsive interactions between casein micelles and entrapped β-LG. In that case, no more β-LG accumulated into the deposit and none are released into the soluble phase. When the deposit is unchanged because the ionic environment of milk is not modified (case of PUF-DF), the β-LG accumulated progressively into the deposit.

As a result, to increase Tr β-LG during ROW-DF, it is preferable to accumulate β-LG into the deposit before the ROW-DF (by applying high CF value) and promote conditions that limit cohesive deposit formation (low TMP). Tr β-LG will be optimal during diafiltration for a DF rate of 2. Extended diafiltration would impair the Tr β-LG and lead to a high volume of diluted permeate that will need to be further concentrated.

Fig. 1 .

 1 Fig. 1. Permeation flux (Jp) and conductivity of retentate during the diafiltration (DF) of skimmed milk retentate with either milk ultrafiltrate (PUF) (a & b) or reverse osmosis water (ROW) (c & d). CF 2.5 and TMP 0.7 bar (solid circles, •); CF 2.8 and TMP 0.7 bar (open circles, ○); CF 2.6 and TMP 0.7 bar (solid squares, ■); CF 3.2 and TMP 0.7 bar (open squares, □); CF 3.0 and TMP 0.9 bar (solid triangles, ▴); CF 3.5 and TMP 0.9 (open triangles, Δ).

  ROW-DF -Jp corr ROW-DF ) (Jp ROW-DF -Jp corr PUF-DF ) ε in two groups of operating conditions (at TMP 0.7 bar and CF = 2.5-2.6 in one hand and at TMP = 0.7 bar at CF = 2.8-3.0 in the other hand), shows that regardless of the CF, the decrease of lactose content in the permeate during ROW-DF explains about 40 % of the increase of Jp (comparison of Fig. 1 c and Fig.

Fig. 2 .

 2 Fig. 2. Transmission of β-LG (a, c & e) and α-LA (b, d & f) during PUF-DF (a & b) or ROW-DF (c-f). Operation conditions. 800 kDa polymeric spiral-wound membrane; retentate recirculating flow, Qr = 18 m 3 ⋅h -1 ; T = 12 • C. CF 2.5 and TMP 0.7 bar (solid circles, •); CF 2.8 and TMP 0.7 bar (open circles, ○); CF 2.6 and TMP 0.7 bar (solid squares, ■); CF 3.2 and TMP 0.7 bar (open squares, □); CF 3.0 and TMP 0.9 bar (solid triangles, ▴); CF 3.5 and TMP 0.9 (open triangles, Δ).

Fig. 3 .

 3 Fig. 3. Content of lactose and viscosity of permeate during diafiltration, DF with reverse osmosis water (ROW). Diafiltration is performed after concentration of skimmed milk at CF = 3.5 (solid triangles, ▴) or 3.0 (open tringles, Δ). Operating conditions. TMP = 0.9 bar.

Fig. 4 .

 4 Fig. 4. Permeation flux corrected by the viscosity of the permeate during the diafiltration (DF) of microfiltration retentate with PUF (a) or ROW (b). Operating conditions. CF 2.5, TMP 0.7 bar (solid circle, •); CF 2.8, TMP 0.7 bar (open circle, ○); CF 2.6, TMP 0.7 bar (solid square, ■); CF 3.2, TMP 0.7 bar (open square, □); CF 3.0, TMP 0.9 bar (solid triangle, ▴); CF 3.5, TMP 0.9 (open triangle, Δ).

Fig. 5 .

 5 Fig. 5. "Surplus" of phosphate in the permeate (mg⋅L -1 ) as a function of diafiltration rate, Tx with reverse osmosis water (ROW) at TMP 0.7 bar and CF 2.5 (solid square, ■) or CF 3.2 (open square, □) or at TMP 0.9 bar and CF 3.0 (open triangle, Δ) or CF 3.5 (solid triangle, ▴).

Fig. 6 .

 6 Fig. 6. Permeation flux corrected by the viscosity of the permeate as function of conductivity of retentate during the diafiltration (DF) with reverse osmosis water (ROW). Operation conditions. CF 2.6, TMP 0.7 bar (solid square, ■); CF 3.2, TMP 0.7 bar (open square, □); CF 3.0, TMP 0.9 bar (solid triangle, ▴); CF 3.5, TMP 0.9 (open triangle, Δ).

Fig. 7 .

 7 Fig. 7. Surplus of β-lactoglobulin (β-LG) in the permeate (g⋅L -1 of diafiltered retentate) as a function of rate of diafiltration (Tx) during PUF-DF (a) or ROW-DF either at TMP = 0.7 bar (b) or 0.9 bar (c).

SW and reverse osmosis membranes. Regardless of their origin, the two PUFs had similar compositions (Table 1). Compared to Table 1

  • C one day prior MF. Before MF, skimmed milk was maintained at 12 • C for 30 min to recover a stable mineral balance. Average composition and standard deviation of skimmed milks and PUF (permeate of ultrafiltration).

	Diafiltration step was carried out with either ROW or PUF. ROW was
	water from network treated by reverse osmosis and electrodialysis
	(Pretreatment pack Progard TL1 CL2, Merck-Millipore) and by UV lamp
	at 25 • C to reach a resistivity of 15 MΩ⋅cm. PUF was obtained by ul-
	trafiltration of skimmed milk at 45 • C. It was either directly prepared
	with a 5 kDa SW membrane (6338 HFK 328 VYT, Koch membrane
	systems, Massachusetts, USA) at TMP of 2.1 -2.6 bar and permeation flux, Jp of 16 -32 L⋅h -1 ⋅m -2 (STLO Dairy Platform) or obtained after 1:4
	dilution of concentrated (27 % of dry matter) PUF. The concentrated
	PUF, provided by Eurosérum (Grand Est, France) was initially obtained
	with 5 kDa Skimmed milk at	PUF produced at
		12 • C	45 • C
	pH Fat, g⋅kg -1 Dry matter, g⋅kg -1 Ashes, g⋅kg -1 Soluble calcium, mg⋅L -1 Soluble phosphate, mg⋅L -1	6.72 ± 0.04 a 1.5 ± 1.0 a 91 ± 1 b 7.6 ± 0.2 b 337 ± 16 c 970 ± 41 a	6.64 ± 0.04 d 0 d 58 ± 3 d 4.9 ± 0.3 d 282 ± 1 d 813 ± 4 d
	Lactose, g⋅L -1 Total Nitrogen, C TN , g⋅kg -1 β-lactoglobulin, β-LG, g⋅L -1	48 ± 2 33.3 ± 1.4 a 3.14 ± 0.13 b 0.94 ± 0.08 b	52 ± 3 1.6 ± 0.1 d 0 d 0 d
	α-lactalbumin, α-LA,
	g⋅L -1
	a n = 6; b n = 5; c n = 4;

d n = 2.

Table 3

 3 Evolution of total phosphate (ΔTotal phosphate) and of soluble phosphate (Δsoluble phosphate) in retentate during diafiltration with reverse osmosis water (ROW-DF) at a concentration factor (CF) of 2.6 and 3.2. Operating condition. TMP 0.7 bar.

	TMP	CF	Tx	Total	Soluble	ΔTotal	ΔSoluble
	(bar)	(-)	(-)	phosphate (mg⋅L -1 )	phosphate (mg⋅L -1 )	phosphate (mg⋅L -1 )	phosphate (mg⋅L -1 )
	0.7	2.6	0.0	2940	867	-	-
			2.0	2134	314	806	553
	0.7	3.2	0.0	3749	888	-	-
			5.7	2305	151	1444	737

Acknowledgements

The authors would like to thanks Gilles Garric and Jean-Luc Thomas (STLO, UMR 1253, INRAE, Institut Agro) for experimental support on the STLO dairy Platform, and Maksym Loginov for fruitful discussions.

Funding

This work was supported by a grant from the Brittany Region (grant no. 16006734, INRA convention 30001292) and from FEDER (grant no. EU000171, INRA convention 30001293).

Data availability

The authors do not have permission to share data.

CRediT authorship contribution statement

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Keywords: Delay, vibration control, frequency-domain, time-varying delay, delay estimation

As in the case of the Smith predictor, when applied to dynamical systems with delay in the input, the stability of the internal model control (IMC) system can be very sensitive to the mismatch that can occur between the plant delay and the identified delay. In this paper, we propose an online identification technique for the plant with input delay, which is used to adapt an IMC controller designed to reject periodic disturbances consisting of two harmonics. The proposed method is based on wavelet analysis and cross-correlation of the measured output. Moreover, the proposed controller structure introduces a novel analytical approach which makes it suitable for the adaptive application. The practical motivation of the study is to address the delay mismatch that can be caused by approximating time-varying high-order systems by first-order models with input delay. The approach is numerically validated.

1. INTRODUCTION Some processes, mainly industrial processes, require tasks to be performed in a repetitive manner. In such scenarios, these processes may be exposed to periodic disturbances inherently occurring due to wearing over time or manufacture fault of the machine, which in return reflects on the output each cycle during the repetitive task. Moreover, control of such industrial processes may involve time-delays stemming from mainly two different sources. First, the physical setup of the process can impose a time-delay due to sensor placement and communication, for instance, in the case of hot rolling mills (Yuksel et al. (2021a)). Second, delays can be introduced artificially by approximating high-order process dynamics by first-order models with input delay. Such an approximation finds wide applications in the industry since, in general, simple time-delay models can sufficiently capture the main dynamics of complex processes required for the control design. Combined together, these two types of delays can make the overall model of the process to have a time-delay larger than the model time-constant. Therefore, the time-delay becomes essential in the design of control.

A controller that can address both the time-delay and the periodic disturbance can be obtained by modifying several different controller structures. For instance, repetitive controllers, which are renowned for their capabilities in periodic output regulations, can be modified to address input time-delays by introducing another delay component to the control-loop as in [START_REF] Omura | Attenuation of roll eccentric disturbance by modified repetitive controllers for steel strip process with transport time delay[END_REF]. Alternatively, the well-known Smith's predictor to compensate time-delays can be modified to reject periodic disturbances as in [START_REF] Zhou | Modified smith predictor design for periodic disturbance rejection[END_REF].

Another promising controller design approach is the Internal Model Control (IMC) framework. For instance, two different IMC controllers to remove periodic disturbances from systems with dead-time, one utilizing a lumped delay (Yuksel et al. (2021a)) and the other a distributed delay [START_REF] Yuksel | Internal model control with distributed-delay-compensator to attenuate multi-harmonic periodic disturbance of time-delay system[END_REF]), were recently proposed in the open literature. In contrast to the repetitive control, IMC employs the model of the controlled plant instead of the signal model. Under ideal settings, i.e. when the plant and the model perfectly match, the controller can be designed as if the control system is of feedforward.

Under non-ideal conditions, however, as similarly observed for Smith predictors [START_REF] Morȃrescu | On the geometry of stability regions of smith predictors subject to delay uncertainty[END_REF]), the characteristic equation of the IMC applied to a plant with input delay corresponds to a quasipolynomial with two distinct exponential terms i.e.

∆(s) = N (s) + M 1 (s)e -sτ1 + M 2 (s)e -sτ2 = 0, (1)
where the polynomials N (s) and M i (s) satisfy deg (N (s)) ≥ deg (M i (s)) and τ 1 , τ 2 > 0, i = 1, 2. Such type of characteristic equations corresponds to those of time-delay systems and the stability of the system is assured if all the roots of (1) are located on the left hand plane of the complex plane [START_REF] Michiels | Stability and stabilization of time-delay systems: an eigenvalue-based approach[END_REF]). It is a known fact that the location of the roots can be very sensitive to the relative difference of the delays [START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF]), leading to a questionable robustness for application. For this reason, delay identification methods and delay adaptive controller designs received attention as an alternative to robust design methods.

As stated in [START_REF] Bjorklund | A review of time-delay estimation techniques[END_REF], delay identification methods can be grouped under four main classes: (i) methods that derive delay from an approximated model, under which cross-correlation [START_REF] Azaria | Time delay estimation by generalized cross correlation methods[END_REF]) and function estimation [START_REF] Gawthrop | Identification of time delays using a polynomial identification method[END_REF]) can be found, (ii) methods treating delay as an explicit parameter [START_REF] Orlov | Adaptive identification of linear time-delay systems[END_REF]), (iii) methods that exploit properties of certain responses [START_REF] Wang | Direct identification of continuous time delay systems from step responses[END_REF]), and finally, (iv) methods with higher-order statistics [START_REF] Nikias | Time delay estimation in unknown gaussian spatially correlated noise[END_REF]).

Forming an adaptive controller for systems with varying or uncertain delays was shown to be possible by replacing the delays with their rational function approximation [START_REF] De Souza | An adaptive control algorithm for linear systems having unknown time delay[END_REF]) or with transport partial differential equations [START_REF] Bresch-Pietri | Adaptive control scheme for uncertain time-delay systems[END_REF]). A delay adaptive Smith predictor was also discussed in [START_REF] Bahill | A simple adaptive smith-predictor for controlling time-delay systems: A tutorial[END_REF].

Problem Statement

The control problem treated in this paper is to design an adaptive IMC controller that allows a plant, for which its only time-varying property is its input delay, track a given step input meanwhile it is subjected to an output disturbance composed of two harmonics i.e. d(t) = a 1 sin(ω 1 t + ϕ 1 ) + a 2 sin(ω 2 t + ϕ 2 )

(2) where a 1 , a 2 are amplitudes, ω 1 , ω 2 are commensurable excitation frequencies and ϕ 1 , ϕ 2 are phase shifts. The interest in such a disturbance model stems from the fact that, in the hot rolling process, the disturbance caused by the eccentricity in the rolls can be usually characterized by the sum of two sinusoidal signals.

The paper is organized as follows. Section 2 introduces the structure of the IMC controller and the adaptive modification rule for tuning the parameters. In Section 3, the simulation results of the proposed adaptive control for a prescribed delay change are presented, and in Section 4, the general points of implementation and design are discussed. Section 5 concludes the paper.

DELAY ADAPTIVE IMC CONTROLLER

Internal Model Control is the control arrangement depicted in Figure 1, in which the model of the controlled plant is employed in the control loop to compare the actual output of the plant with the estimated one.

IMC was initially and extensively studied by [START_REF] Garcia | Internal model control. a unifying review and some new results[END_REF]. Further extensions of the IMC with various adaptive methods were introduced in [START_REF] Gao | An adaptive fuzzy smith control of time-varying processes with dominant and variable delay[END_REF] and also in [START_REF] Marino | Robust adaptive compensation of periodic disturbances with unknown frequency[END_REF] for periodic disturbance rejection and in [START_REF] Hahn | Adaptive imc control for drug infusion for biological systems[END_REF] for delay compensation. A robust filter-based IMC design for removing periodic disturbance from systems approximated by first-order models with dead-time was proposed initially in [START_REF] Vyhlídal | Control system design based on a universal first order model with time delays[END_REF]. We consider the plant and its model to be composed of two parts, namely, invertible and non-invertible parts. The invertible part, denoted by G p and G m for the plant and the model respectively, encapsulates the first-order dynamics. The noninvertible part corresponds to the delays τ p and τ m of the plant and the model, respectively. The controller is treated in a similar manner where the transfer function is denoted by Q and the user-definable time-delay by θ.

The sensitivity transfer function S(s) of the scheme in Fig. 1 is given as

S(s) = 1 -Q(s)G m (s)e -s(τm+θ) 1 + Q(s)(G p (s)e -s(τp+θ) -G m (s)e -s(τm+θ) ) . (3)
Following the proposed design by [START_REF] Vyhlídal | Control system design based on a universal first order model with time delays[END_REF], the transfer function of the controller is given by

Q(s) = 1 G m (s) F (s), (4)
where F (s) is the filter that needs to be strictly proper.

Non-adaptive Controller Design

To form a non-adaptive IMC controller that asymptotically tracks the reference r and rejects a periodic disturbance of the form (2), we make the following assumptions:

(i) plant and model transfer functions G p (s) and G m (s) are identical i.e.

G p (s) = G m (s) = K T s + 1 (5)
where K and T are the gain and the time-constant of the dynamics, respectively. (ii) plant and model time delays τ p and τ m are identical (denoted only by τ ) (iii) controller Q(s) is of the form (4) If all three conditions mentioned above are met, the sensitivity (3) can be simplified to

S(s) = 1 -F (s)e -s(τ +θ) .
(6) In order to track a step input and cancel periodic disturbance (2) we need to make sure that S(0) = S(jω i ) = 0, i = 1, 2. This can be achieved by fulfilling the following set of conditions:

lim s→0 F (s) = 1, (7)
|F (jω i )| = 1, (8) arg F (jω i ) -ω i (τ + θ) = 2kπ, (9)
Preprints, Joint SSSC, TDS, LPVS 2022 September 27-30, 2022, Montreal, Canada

where k ∈ Z, i = 1, 2. To satisfy the conditions, we propose the filter F in the form of

F (s) = p(s) -z(s) p(s) , (10)
where polynomials z(s) and p(s) are defined as

z(s) = s(s 2 + ω 2 1 )(s 2 + ω 2 2 ), (11)
p(s) = 5 i=1 (s -p i ), (12)
p i ∈ C with Re(p i ) < 0.
Consequently, the filter satisfies conditions ( 7) and (8) since

F (jω i ) = p(jω i ) -z(jω i ) p(jω i ) = p(jω i ) p(jω i ) = 1, i = 1, 2 (13)
which simplifies the condition (9) to

ω i (τ + θ) = 2l i π, i = 1, 2, ( 14
) where l 1 , l 2 ∈ N, due to arg F (jω i ) = arctan Im(F (jω i )) Re(F (jω i )) = 0, i = 1, 2. (15)
Moreover, since ω 1 and ω 2 are commensurable, i.e. ∃k 1 , k 2 ∈ Z such that ω1 ω2 = k1 k2 with gcd(k 1 , k 2 ) = 1, the base frequency defined as ω b := ωi ki yields an equation for the smallest θ satisfying ( 14):

θ = 2l b π ω b -τ ≥ 0, (16)
where l b = τ ω b 2π .

Online Delay Identification

Qe

ỹ(t) = [ỹ(t-((N -1)∆t), [ỹ(t-(N -2)∆t), . . . , ỹ(t)] (18)
of the controlled system within a sliding window

T W = N • ∆t: ϑ[n φ ] = (u ⋆ ỹ)[n φ ] = N i=1 u[i] • ỹ[i + n φ ]. (19)
Where n φ ∈ [-N, N ] ⊂ Z (because we use full implementation of cross-correlation) is a discrete-time shift between the signals, N is a number of signal samples in the sliding window, and ∆t is a discretization time step that satisfies ∆t ≪ τ and T W ≫ φ(ω)/ω, where φ(ω) is a phase shift of the controlled system G p

φ(ω) = arg(G(jω)e -jωτ ) = -τ ω -arctan(T ω). ( 20
)
The estimated output ỹ(t) = y(t) -d(t) is a measured output y(t) subtracted by the identified disturbance d(t). The disturbance identification is described in the next section.

According to (20) the cross-correlation consists of the delay influence and time constant influence. For the separation of these two parameters, it is also necessary to know the frequency of the signal ω = 2πf = 2π/T ω . The method is derived for single-frequency signal. However, when the signal contains more than one frequency, the cross-correlation is influenced. To avoid the influence of the second frequency, we filter the inputs to the cross-correlation using the low-pass Butterworth filter, where the cutoff frequency ω c must satisfy ω 1 < ω c < ω 2 .

After filtering, we use cross-correlation only for the frequency ω 1 .

An example of the result of the cross-correlation of periodic signals is visualized in Fig. 3

(ω) = -τ0 • ω -arctan(T • ω),
where ω is the identified frequency, and from a period range modificator C:

A = φ0 (ω) ω -C T ω 2 1 ∆t , C ∈ (0, 1) (22)
and

B = φ0 (ω) ω + C T ω 2 1 ∆t , C ∈ (0, 1). ( 23
)
From ( 20) and ( 21) it is possible to express the identified time delay τ :

τ = ñ • ∆t + arctan(T • ω) ω . ( 24
)
The resolution of the peak position is ∆t. To improve the accuracy of the peak position, we interpolate the neighborhood of the identified τ with a cubic spline (visible in the lower graph in Fig. 3). The finer resolution is then possible to obtain by finding the maximum of such a spline on the interval (τ -∆t, τ + ∆t).

Disturbance identification

The disturbance is identified when the system runs without control. The approach is similar to the time delay tracking. It is based on tracking small inaccuracies in the initial estimate.

In our case, we assume that the disturbance has character according to (2). There are six parameters to identify: ω 1 , ω 2 , a 1 , a 2 , ϕ 1 , ϕ 2 . For the identification of disturbance frequencies ω 1 , ω 2 and phase shifts ϕ 1 , ϕ 2 , the discrete convolution of an adaptive Morlet wavelet and an output signal y(t) is used:

Γ[n, ω] = (y * Ψ)[n, ω] = N k=1 y[k] • Ψ[n -k, ω]. ( 25
)
Where N is the number of samples in the sliding window and ω is changing in the desired range of frequency identification.

The wavelet is a symmetrically damped cosine function:

Ψ[k, ω] = exp (-k∆t) 2 2σ[ω] 2 • cos (ω • k∆t). ( 26
)
The damping is proportional to the period of oscillations:

σ[ω] = 2 • T ω = 4π ω . (27)
The resulting wavelet should have a sum equal to 1 for each ω. Thus, every wavelet is divided by its sum:

Ψ[k, ω] = Ψ[k, ω] N k=1 Ψ[k, ω] . (28)
The identified frequency is considered that ω, for which the sum of squared convolution Γ[n, ω] (25) is maximized:

ω = arg max ω N k=1 Γ 2 [k, ω]. (29)
The identified time shift is calculated from the position of a peak in convolution for identified ω:

ñ = arg max n Γ[n, ω], (30)
where the shift in radians is:

k • φ = ñ∆t - T ω 4 • ω, k ∈ Z. (31)
One wavelet can identify one ω and one φ. Thus, it is necessary to use two wavelets to identify both sine components.

After identification of the frequencies and phase shifts of both sine components, the amplitudes a 1 , a 2 are identified using the discretized vector of y(t) and the discretized vector of d(t), which is calculated by substituting ω1 , ω2 , φ1 , φ2 in equation for disturbance (2). The amplitudes a 1 and a 2 are identified by sweeping around the expected values minimizing the sum of squared error between measured signal and the model of disturbance:

a 1 , a 2 = arg min a1,a2 N k=1 y[k] -d[k, a 1 , a 2 , ω1 , ω2 , φ1 , φ2 ] 2 .
(32)

Algorithm

The simplified algorithm scheme is shown in Fig. 4. To identify the initial disturbance, it is necessary to run the system without control. After identifying the disturbance, the controller is activated and it is possible to track the time shift between the input to the system from the controller u(t) and the output of the system subtracted by the estimated disturbance ỹ(t) = y(t) -d(t). If the average output values y(t) in the sliding window T W are damped below the user-defined limit ȳ(t) < y lim and the change in the identified delay τ (t) is less than the user-defined limit ∆τ (t) < τ lim , then the identified delay is held as the maximum value τ (t) = max(τ (t)), where τ is: τ (t) = [τ (t-((M -1)∆t), τ (t-(M -2)∆t), . . . , τ (t)] (33) in the sliding window of length:

M = 2 • T ω ∆t . (34)
We hold the maximum value because the identified τ (t) oscillates. The oscillations are caused by the mismatch of the length of sliding window and the integer multiple of the disturbance period T ω .

NUMERICAL VALIDATION

The proposed method was simulated in Matlab-Simulink using the ODE45 solver, with max. step size equal to ∆t and relative tolerance equal to 10 -9 . The exact parameters of the system, controller, and adaptive algorithm can be seen in Table 1.

The first five seconds of simulation are without control, as you can see in Fig. 5, then the controller is activated. The evolution of the variable delay is shown in Fig. 6. The delay is equal to the initial estimate at the beginning of the simulation and the controller is able to suppress the disturbance very quickly. The delay changes after 15 seconds of simulation by 10 % with a slope of 0.1, and after 30 seconds it increases by 20 % with the

T [s] 0.017 K [-] 0.587 τ [s] 0.211 ∆t [s] 0.01 ω = [ω 1 ω 2 ] [rad•s -1 ] [8π 16π] θ [s] 0.039 C [-] 0.2 y lim [mm] 0.05 τ lim [-] 0.005 T W [s] 5
same slope. The speed of tracking the delay change depends on the length of the sliding window T W .

The response to the delay change is promising. The oscillations in the case of a 10 % change are less than 1.5 mm, which is much smaller than the disturbance amplitude (3.5 mm), and the oscillations are damped in less than 5 seconds (Fig. 7). As can be seen in Fig. 8, the oscillations caused by the 20 % change are less than 15 mm, which is more than the original disturbance, but they are also damped quickly in less than 10 seconds.

DISCUSSION

The simulation in Matlab-Simulink provides promising results, but there are few points for further discussion. In this case, we simulated only a variable time delay, and the adaptive wavelet was used only for initial disturbance identification. The disturbance features are considered constant and without any other noise. Our approach needs the update of the disturbance parameters, because even for small inaccuracies of disturbance identification, the time delay tracking can fail in a long time Simulation with 20 % delay change Fig. 8. Response to 20 % delay change after 30 seconds of simulation. frame. Another problem can occur in case the frequencies of the disturbance are close to each other. Then it would be necessary to pay attention to Butterworth filter design. This method cannot deal with disturbance and time delay identification at the same time. Another problem is the separation of the impact of the system time constant T , which is considered known and constant. Time delay tracking would fail if the delay change is greater than the base period T ω , because the algorithm would track a different position of the peak in cross-correlation. However, the main purpose of the application of this algorithm is to handle small and slow changes in time delay because even small changes can lead to instability.

CONCLUSION

An adaptive IMC controller design to suppress two-harmonic disturbances acting on a first-order system with dead-time has been proposed and numerically validated. The proposed filter structure for the IMC controller explicitly distinguishes the controller parameters responsible for delay compensation and disturbance rejection, allowing fast adaptation of the controller parameters for a newly identified system delay.

  Two-harmonic Periodic Disturbance Compensation for Dead-Time Processes with Online Identification of Input Delay ⋆ Adam Peichl * Michal Kuchař * Can Kutlu Yuksel * , * * Silviu-Iulian Niculescu * * Tomáš Vyhlídal *

Fig. 1 .

 1 Fig. 1. Conventional internal model control scheme for periodic disturbance compensation

Fig. 2 .

 2 Fig. 2. Proposed adaptive internal model control scheme. We propose an adaptive scheme shown in Fig. 2 that tracks the changes in time delay and updates the delay parameters in the controller θ and in the model τ m . The adaptation of the time delay τ m is based on the discrete cross-correlation of discretized input: u(t) = [u(t-((N -1)∆t), u(t-(N -2)∆t), . . . , u(t)] (17) and estimated output: ỹ(t) = [ỹ(t-((N -1)∆t), [ỹ(t-(N -2)∆t), . . . , ỹ(t)] (18)

  of the simulation example which follows. The result has several peaks where the position of each peak represents the total time shift between the signals. The distance between these peaks corresponds to the base period of the signals T ω = 2π/ω 1 . For a successful identification of the time shift, it is necessary to have an a priori estimate of the time delay τ0 . The error of the time delay estimate must satisfy |τ -τ0 | ≪ T ω /2. Otherwise, the identified delay would have an error k • T ω , k ∈ N. It is possible to express the total time shift between the signals as the position of a peak: ñ = argmax nφ |ϑ[n φ ]| s.t. A ≤ n φ ≤ B (21) in a specified range of time shift n φ ∈ [A, B] ⊂ Z using the estimate of time delay, which should narrow the crosscorrelation range to identify only one position of a peak. The range is calculated from the intital total time shift estimate φ0

Fig. 3 .

 3 Fig. 3. Cross-correlation of two shifted periodic signals within sliding window T W = 1 s and time-step ∆t = 0.01 s (upper graph) and the detail of estimated range of the total time shift (lower graph).

Fig. 7 .

 7 Fig. 5. Simulation without control at the beginning. Control started after 5 seconds.

Table 1 .

 1 Parameter values for simulation

	START
	Run system without control
	Identify disturbance
	parameters using wavelet
	Run controller
	Cross-correlation of and
	for time delay tracking
	Is	No
	and	?
	Yes
	Hold max. time delay value in
	sliding window
	Yes
	End of simulation?
	No
		No
	Is	?
	Yes
	STOP the
	simulation
	Fig. 4. Simplified algorithm of time delay tracking
	Parameter:	Value:
	p = [p 1 p 2 p 3 p 4 p 5 ]	[-15 -30 -45 -60 -75]

Preprints, Joint 8th IFAC Symposium on System Structure and Control, 17th IFAC Workshop on Time Delay Systems, 5th IFAC Workshop on Linear Parameter Varying Systems. Montreal, Canada, September 27-30, 2022 © 2022 the authors. Accepted by IFAC for publication under a Creative Commons License CC-BY-NC-ND

Preprints, Joint SSSC, TDS, LPVS 2022 September 27-30, 2022, Montreal, Canada

Preprints, Joint SSSC,TDS, 2022 September 27-30, 2022, Montreal, Canada

⋆ This work was supported by the Czech Science Foundation project No. 21-07321S -Persistent problems of repetitive control. The work of the first three authors was also supported by the Grant Agency of the Czech Technical University in Prague, student grants No. SGS21/152/OHK2/3T/12 and SGS20/159/OHK2/3T/12. The work of the third and fourth authors was also supported by a public grant overseen by the French National research Agency (ANR) as part of the Investissements d'Avenir program, through the "ADI 2020" project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02.


Maxime Ryckewaert

Anis Taleb Bendiab

Fabienne Henriot

Ryckewaert Maxime

Lamour Julien

Taleb-Bendiab Anis

Bendoula Ryad

Removing spatial effects of

Keywords: Chemometrics, experimental design, maize, spectrometry, reflectance spectra

During an experiment, the presence of spatial autocorrelation can undermine plant response interpretation [START_REF] Ball | Spatial heterogeneity affects variety trial interpretation[END_REF][START_REF] Legendre | Effects of spatial structures on the results of field experiments[END_REF]. Even by using blocking design and randomization, measurement can have a spatial dependence [START_REF] Oehlert | Early drought stress detection in cereals: simplex volume maximisation for hyperspectral image analysis[END_REF]. This may reduce the utility of applying the classical analysis of variance (ANOVA). Generalized linear mixed models (GLMMs) can handle the spatial dependence as random effects [START_REF] Bolker | Generalized linear mixed models: a practical guide for ecology and evolution[END_REF]. For multiple responses, GLMMs can be used as in longitudinal analysis [START_REF] Liang | Longitudinal data analysis using generalized linear models[END_REF][START_REF] Zeger | Longitudinal data analysis for discrete and continuous outcomes[END_REF]. But GLMMs are generally inappropriate for multivariate data especially with the spectral data.

In agronomy, spectral data is increasingly being used [START_REF] Araus | Field high-throughput phenotyping: the new crop breeding frontier[END_REF]. High-spectral resolution provides more information and can be valuable for disease detection (Mahlein et al., 2012), for crops monitoring (Ecarnot et al., 2013;Mahajan et al., 2016;Vigneau et al., 2011) and early detection of water-stress [START_REF] Behmann | Detection of early plant stress responses in hyperspectral images[END_REF]Römer et al., 2012). To take advantage of an experimental design with multivariate data such as spectral data, chemometrics methods have been developed [START_REF] Brereton | Chemometrics in analytical chemistry-part I: history, experimental design and data analysis tools[END_REF]. Among these methods, ANOVA-Simulatenous component analysis (ASCA) [START_REF] Smilde | ANOVA-simultaneous component analysis (ASCA): a new tool for analyzing designed metabolomics data[END_REF] gives statistical significance of the factors and can reveal their interdependencies.

In this presentation, we describe the possibility to use such approaches in agronomy to study spatial effect of maize spectra included in an experimental design. An example is given by applying ASCA method on a dataset of 480 spectra where 10 genotypes are confronted by 2 irrigation treatment. Spectra are acquired with highspectral resolution (256 bands) in the range of 310 nm -1100 nm. The obtained results confirm that this class of methods is relevant to identify variability due to spatial effect. Describing the spectral region due to this effect extends the possibility to remove unwanted effects due to spatial dependence.


Maxime Ryckewaert

Nathalie Gorretta

Henriot Fabienne

Alexia Gobrecht

Jean-Michel Roger

Jean-Michel Roger Mul- Tivariate

Ryckewaert Maxime

Roger Jean- Michel

Multivariate analysis of variance of vegetation spectra dataset included into an experimental design by using ANOVA-SCA and ANOVA-Target Projection

Keywords: Chemometrics, ASCA, ANOVA-TP, experimental design, spectrometry, reflectance spectra, agronomy, maize

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Progress made in the design of spectrometric equipment but also in the processing of associated spectral data has democratized the use of spectrometry for experimental acquisitions. In agronomy, spectral data has thus demonstrated its value for disease detection [START_REF] Mahlein | Hyperspectral imaging for small-scale analysis of symptoms caused by different sugar beet diseases[END_REF] or even for crops monitoring [START_REF] Ecarnot | Assessing leaf nitrogen content and leaf mass per unit area of wheat in the field throughout plant cycle with a portable spectrometer[END_REF][START_REF] Mahajan | Monitoring nitrogen, phosphorus and sulphur in hybrid rice (Oryza sativa L.) using hyperspectral remote sensing[END_REF][START_REF] Vigneau | Potential of field hyperspectral imaging as a non destructive method to assess leaf nitrogen content in Wheat[END_REF]. Consequently, spectral data are expected to increasingly be associated with experimental designs. An experimental design is designed to study the variables involved, called factors. In agronomy, these factors can be genotype, treatment or time. An experimental design allows as well, to separate and evaluate the effects of these factors on variance [START_REF] Kirk | Experimental design[END_REF] with a minimum of experiments. Especially for multivariate data where all sources of variability may affect conclusions [START_REF] Massart | Handbook of chemometrics and qualimetrics, Data handling in science and technology[END_REF].

Lately, chemometrics methods have been developed to analyze multivariate data fully integrated to an experimental design [START_REF] Brereton | Chemometrics in analytical chemistry-part I: history, experimental design and data analysis tools[END_REF]. Initially used in metabolomics as NMR spectroscopy and chromatography, ANOVA-Simultaneous component analysis (ASCA) [START_REF] Smilde | ANOVA-simultaneous component analysis (ASCA): a new tool for analyzing designed metabolomics data[END_REF] and ANOVA-Target-Projection (ANOVA-TP) [START_REF] Marini | Analysis of variance of designed chromatographic data sets: The analysis of variance-target projection approach[END_REF] provide statistical significance of every factor and information about their effect on the variability of the data. We are interested here in the use of these techniques for the study of vegetation spectra associated with an appropriate experimental design.

In this presentation, we first discuss the assumptions involved in ASCA and ANOVA-TP. Then, we compare these methods to analyze a dataset of 480 vegetation spectra where 10 genotypes of maize are confronted to 2 irrigation treatments. Spectra are acquired above the crops with a high-spectral resolution (256 bands) in the range of 310 nm -1100 nm. The obtained results confirm that both methods, with their own specificities, are relevant to describe the influences of factors on the variance. ASCA and ANOVA-TP share many similarities and provide the same conclusions concerning the significance of each factor. However some differences are highlighted in terms of scores and loadings. Indeed, ASCA provides orthogonal loadings whereas ANOVA-TP gives accurate estimation of between and within-class variance.


Ryckewaert Maxime

Alexia Gobrecht

Morel Julien

Roger Jean-Michel

Henriot Fabienne

Nathalie Gorretta

The impact of the spatial resolution of highly resolved spectral data on pansharpening methods to reconstruct a hyperspectral image

Keywords: Pan-sharpening, precision agriculture, spectrometry, hyperspectral imaging, plant phenotyping

Many precision agriculture applications require acquisition and processing of high spectral and spatial data from remote sensing to capture the subtle differences caused by plant physiological responses and their variability within a crop [START_REF] Delalieux | Unmixing-Based Fusion of Hyperspatial and Hyperspectral Airborne Imagery for Early Detection of Vegetation Stress[END_REF]. This is particularly the case in plant phenotyping.

Even if satellite images are available, their spatial and spectral resolutions are not appropriate for plant phenotyping. Pan-sharpening methods have therefore been developed to increase the spatial resolution of the multispectral information by fusing a panchromatic image (i.e. high spatial /low spectral resolution) with a multispectral one [START_REF] Vivone | A Critical Comparison Among Pansharpening Algorithms[END_REF]. If hyperspectral satellite sensors are about to be available, their spatial resolution is limited [START_REF] Loncan | Hyperspectral pansharpening: a review[END_REF][START_REF] Yokoya | Hyperspectral and multispectral data fusion: A comparative review[END_REF] and the cost of the sensors is still too high for agricultural applications. In order to overcome the limit of the spatial resolution and to reduce the price of the acquisition, sensors can be used in proximal detection (i.e. at shorter distances to the scene of interest), hand-held or embedded on mobile platforms [START_REF] Araus | New Technologies for Phenotyping[END_REF][START_REF] Deery | Proximal Remote Sensing Buggies and Potential Applications for Field-Based Phenotyping[END_REF][START_REF] Sankaran | Low-altitude, high-resolution aerial imaging systems for row and field crop phenotyping: A review[END_REF]. The objective of this work is to assess the potential of using pan-sharpening methods available in the literature, which are algorithms dedicated to improve the spectral and spatial resolution of multispectral satellite images, on close-range images. And the first issue addressed here is the impact of the spatial resolution of the spectral data on the quality of the reconstructed image.

To do so, a hyperspectral image (from HySpex VNIR 1600) of a sugar beet plot approximately acquired at a height of 1 meter with high-spatial resolution 0.47 mm/pixel and high-spectral resolution (160 bands) in the range of 409 nm -987 nm has been considered as a reference image. From this reference, a high-spatial resolved panchromatic image is obtained by spectral degradation. On the other hand, the spectral information has been extracted by spatial degradation in a regular grid of 𝑁x𝑁 pixels giving a spectra number defined by 𝑁 2 , where 𝑁 belongs to the interval [2,20]. The spectrum assigned to the new pixel is the mean of the spectra contained in the reference hyperspectral image.

The fused images are then obtained from different pan-sharpening algorithms organized into three categories : Component substitution methods, Multiresolution analysis methods and subspace-based methods. To assess the quality of the fused image in comparison to the reference image, different figures of merit are computed : a spatial index with the cross correlation (CC), a spectral index with the spectral angle mapper (SAM) [START_REF] Yuhas | Discrimination among semi-arid landscape endmembers using the spectral angle mapper (SAM) algorithm[END_REF] and finally two integrated indices with the root mean-square error (RMSE) and the erreur relative globale adimensionnelle de synthèse (ERGAS) [START_REF] Wald | Quality of high resolution synthesised images: Is there a simple criterion?[END_REF].

From 𝑁= 10 (i.e. 100 spectra) where the new pixel size corresponds to 51x51 pixels of the panchromatic image, the values of the quality indices don't greatly improve. For this spatial resolution, the three best methods are the smoothing filter-based intensity modulation (SFIM) [START_REF] Liu | Smoothing filter-based intensity modulation: A spectral preserve image fusion technique for improving spatial details[END_REF], Brovey transform (BT) [START_REF] Gillespie | Color enhancement of highly correlated images. II. Channel ratio and "chromaticity" transformation techniques[END_REF] and the coupled nonnegative matrix factorization (CNMF) [START_REF] Yokoya | Coupled nonnegative matrix factorization unmixing for hyperspectral and multispectral data fusion[END_REF] with respectively 6.65°, 7.49°, 6.07° for the SAM index and 0.794, 0.835, 0.949 for the CC index. For all values of 𝑁, the CNMF systematically outperforms the other methods in the spectral and spatial domain confirming the robustness of this method when dealing with hyperspectral data. These promising results are useful to choose the best combination of available cheap sensors.