Aldjia Ait Chekdid
email: aitchekdida@univ-lorraine.fr

Cyril J F Kahn
email: cyril.kahn@univ-lorraine.fr

Eulalie Prévot
email: eulalie.prevot@sthubert.fr

Mathieu Ferrières
email: ferrièresm@sthubert.fr

Béatrice Lemois
email: lemoisb@sthubert.fr

Catherine Choquet
email: catherine.choquet@sthubert.fr

Michel Linder
email: michel.linder@univ-lorraine.fr

Ait Chekdid Aldjia

Ferrières Mathieu
email: mathieu.ferrieres@sthubert.fr

Lemois Béatrice
email: beatrice.lemois@sthubert.fr

Choquet Catherine

Linder Michel

Mixture design applied for formulation and characterization of vegetal-based fermented products

Keywords: flour, fermentation, texture, mixture design

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Today, people are aware of the health benefits of plant-based diets, the improvement and sustainability of these foods and their impact on the environment (Jeske, Zannini, & Arendt, 2018) . There are different diets around the world, and there is an increasing trend towards diets with less or no dairy products for health reasons such as lactose intolerance or allergy problems or simply as a lifestyle choice.

According to de Gavelle (2019), the main protein intake in western diets is essentially of animal origin and it should be noted that milk has long been considered as the only food containing all the essential substances for human nutrition. [START_REF] Davoodi | Effects of Milk and Milk Products Consumption on Cancer : A Review[END_REF] had grouped a number of studies that revealed the fact that certain constituents of milk and common contaminants (such as pesticides, estrogens and insulin-like growth factor 1) may be responsible for adverse health effects. Widely disseminated nutrition information campaigns in the media and emerging studies raise awareness of the importance of nutrition and reinforce the close link between food and health. The potential disadvantage of commercially available plant substitutes is the variation in nutrient profiles, with low quality or protein and mineral intakes. To promote these products as substitutes of cow's milk, manufacturers must consider protein content and quality as well as enrichment. Attention should be drawn to the possibilities to improve them by selecting raw materials based on their nutritional quality [START_REF] Mäkinen | Foods for Special Dietary Needs : Non-dairy Plant-based Milk Substitutes and Fermented Dairy-type Products[END_REF].

Flour from different grains and seeds has been widely used in preparations because of its health benefits. In addition to fibre, these flours provide vitamins, proteins, minerals and carbohydrates [START_REF] González-Montemayor | Prosopis spp. Functional activities and its applications in bakery products[END_REF]. Flours are the result from the grinding of cereal seeds, leguminous plants, certain mealy plants, or other parts of plants dried and/or de-oiled (fruits, leaves, nuts...). The Functional properties of flours reflect the complex interaction between the composition, structure, molecular conformation, these physico-chemical characteristics and the nature of the environment, in which they are associated and measured [START_REF] Kumar | Functional properties of food commodities (wheat, kidney bean, cowpea, turnip, cauliflower) flours[END_REF].

Fermentation is a natural means of improvement, without the addition of additives or preservatives in the food matrix. Thus, the industry uses defined starting cultures to better control the quality of these products and ensure their reproducibility [START_REF] Coda | Yogurt-Like Beverages Made With Cereals[END_REF]. There are various studies [START_REF] Ferragut | Physical characteristics during storage of soy yogurt made from ultra-high pressure homogenized soymilk[END_REF][START_REF] İçier | Changes on some quality characteristics of fermented soy milk beverage with added apple juice[END_REF][START_REF] Pandey | Optimization of the prebiotic & probiotic concentration and incubation temperature for the preparation of synbiotic soy yoghurt using response surface methodology[END_REF]; R. [START_REF] Wang | Effects of endogenous small molecular compounds on the rheological properties, texture and microstructure of soymilk coagulum : Removal of phytate using ultrafiltration[END_REF][START_REF] Da Silva Fernandes | Evaluation of the isoflavone and total phenolic contents of kefirfermented soymilk storage and after the in vitro digestive system simulation[END_REF][START_REF] Patrignani | Technological potential of Bifidobacterium aesculapii strains for fermented soymilk production[END_REF] on fermented plant products, which are generally soybean-based. Despite this effervescence and a multitude of studies highlighting the benefits of soybeans, the latter remains at the heart of some controversies, notably because 82% of production is derived from genetically modified crops [START_REF] James | Global Status of Commercialized Biotech/GM Crops[END_REF]. However, there are various alternatives in addition to soya, like cereals, legumes, other plants or even vegetable proteins [START_REF] Hickisch | Influence of lupin-based milk alternative heat treatment and exopolysaccharide-producing lactic acid bacteria on the physical characteristics of lupin-based yogurt alternatives[END_REF][START_REF] Akin | Functional properties of fermented milk produced with plant proteins[END_REF][START_REF] Gupta | Development of fermented oat flour beverage as a potential probiotic vehicle[END_REF][START_REF] Wang | Evaluation of chickpea as alternative to soy in plant-based beverages, fresh and fermented[END_REF][START_REF] Zannini | Development of novel quinoabased yoghurt fermented with dextran producer Weissella cibaria MG1[END_REF]. As defined by [START_REF] Coda | Yogurt-Like Beverages Made With Cereals[END_REF], the fermentation of cerealbased beverages by LAB (Lactic Acid Bacteria) improves protein digestibility, nutritional bioavailability of minerals, bio-accessibility of other nutrients (polyphenols, fiber), decreases the glycemic index, increases shelf life and improves the organoleptic quality of products. To achieve the nutritional and textural quality appreciated by the consumer, a formulation system is needed that links the impact of each compound and their interaction. In contrast to the classical approach with the study of one variable at a time, the use of a mixture design in this type of formulation allows for statistically-based optimization to better understand the dynamics of the process and to illustrate the effects and complex interactions between the constituents [START_REF] Gupta | Development of fermented oat flour beverage as a potential probiotic vehicle[END_REF]. Therefore, the answer depends on the proportions of the different ingredients in the mixture.

The aim of this study is focused on two parts: the first one is the characterization of flours in terms of composition, physico-chemical characteristics and fermentability. In the second one, the creation and realization of a mixture design with three of these flours preselected for the formulation of fermented desserts as an alternative to dairy products. Thus, the fermentability of these mixtures, their effects on the physico-chemical parameters and on the texture of the fermented products, the microbiological quality and the viability of the lactic acid bacteria during the shelf life were developed.

Materials and methods

Raw material

Flours from oat (Avena sativa L.), rice (Oryza Sativa L.), quinoa (Chenopodium quinoa) of the Celnat brand (Saint-Germain Laprade, France), coconut flour (Cocos nucifera) and lupin (Lupinus albus) of the brand Ecoidées (Soultz Sous Forets, France) and chickpea (Cicer arietinum) of the brand Priméal (Peaugres, France), has been selected for this study. All were from organic and fair farming.

Chemicals and reagents which were of analytical quality have been obtained from suppliers such as Sigma Aldrich (Schnelldorf, Germany), Lybios (Pontcharra-Sur-Turdine, France), Danisco France (Paris, France), and Biokar (Solabia group, France).

Chemical properties

The dry matter and ash content were determined by the standard method (AOAC, 1995). The content of starch and sugar was determined by enzymatic kit Total starch HK and K-SUFRG kit (Megazyme, Lybios, Pontcharra-Sur-Turdine, France). The protein content was analyzed by the Kjeldahl method (using N×6.25 as a conversion factor) (AOAC, 1995). The fat was extracted and measured by the Soxhlet method using a mixture of diethyl ether and petroleum ether (1:1 v/v) as solvent.

Physical properties

Hydration properties

To characterize the hydration properties of the flours, the water absorption index (WAI), the water solubility index (WSI) and the swelling power (SP) were determined using the method described in [START_REF] Kraithong | Physicochemical and functional properties of Thai organic rice flour[END_REF] and [START_REF] Lai | Effects of hydrothermal treatment on the physicochemical properties of pregelatinized rice four[END_REF].

A 2.5 g sample of flour was dispersed in 30 mL of distilled water and heated at two temperatures 30 °C and 90 °C independently, for 30 min in a water bath with stirring using a vortex mixer every 10 min. The cooked dough was cooled to room temperature and then centrifuged at 5732 g for 15 min (Centrifuge Heraeus Megafuge 8R, Fisher Scientific, Göteborg -Sweden). The mass of dry solids in the supernatant was recovered by evaporation overnight at 105±1 °C and the pellet was weighed. The characteristic indexes were calculated by following the Eq. ( 1), ( 2), (3).

WAI=

Mass of the pellet g Mass of the sample (g)

(1)

WSI ( g 100 ⁄ g)=
Mass of dissolved solids in the supernatant g x100

Mass of the sample (g)

(2)

SP = WAI × 1 - (3)

Analysis of the particle size

Mastersizer 3000 Laser Diffraction Scale provides particle size distributions by using the Aero S dry dispersion unit (Malvern Instruments, United Kingdom).

The flour was put into the dry dispersion unit and the measurements were performed in triplicate, with a refractive index of 1.52, 2.0*10 5 -3.0*10 5 Pa of pressure, and opening of the hopper of 3 mm. Then, the particle size distribution by decile (D10, D50, and D90) was recorded.

Differential scanning calorimetry (DSC)

The thermal characteristics of flours have been studied using differential scanning calorimetry (DSC250, TA Instruments, USA) calibrated using an empty pan.

For each flour 2±0.5 mg were dispersed in distilled water to achieve the ratio of (30/70 w/w) (flour: water), in a hermetically sealed aluminum capsule as it was described by [START_REF] Villanueva | Microwave absorption capacity of rice flour. Impact of the radiation on rice flour microstructure, thermal and viscometric properties[END_REF]. The samples were then held for 5 min at 30 °C and then heated to 90 °C at a rate of 20 °C/min.

The TRIOS software was used to determine the starting temperature (T0), the peak temperature (Tp), and the ending temperature of the gelatinization (Tf), as well as the ΔHGel (gelatinization enthalpy).

Flour fermentation

We followed the kinetics of acidification of the different flours individually, controlling the decrease of pH of a solution (flour + water) heat-treated before being inoculated with lactic acid bacteria (LAB) strains composed of the species Streptococcus thermophilus and Lactobacillus delbrueckii subsp. bulgaricus (YO-MIX™ 305, Danisco, Paris, France).

The inoculated tubes were placed in a water bath at a temperature of 45 °C and the pH variations were recorded by a D291 consort module (Consort bvba, Turnhout, Belgium), with a 30 min interval between two measurements over a total duration of 21 h.

Experimental plan

Mixture design

The design of the mixture plan was carried out with the NemrodW software. Thirteen mixtures were obtained in the plan with 3 constituents (Table 1) according to a Scheffé matrix with lower ratio constraints (40 % for chickpea, 10 % for oat, 20 % for coconut) and higher ratio constraints (60 % for chickpea, 20 % for oat, 40 % for coconut) from dry matter. The centered cubic reduced mathematical model allows us to predict the values of the experimental responses from the iso-response curves, for water retention, color difference, and texture (hardness and cohesion). Eq. ( 4) shows the cubic model that was fitted to each response, with bi representing the coefficients, and xi representing the values of the factors.

y = b * x + b * x + b * x + b * x * x + b * x * x + b * x * x + b * x * x * x (4)

Preparation of the mixtures

The flours were rehydrated with tap water at 50 °C in the proportion indicated in the mixture design table to obtain a 15 g/100 mL solution of dry matter in the final formulation. The heat treatment applied to the different samples was pasteurization at 90 °C for 10 min, using a

Stephan device (Stephan Machinery GmbH, Hameln, Germany). The products were then rapidly cooled to 45 °C, which is the inoculation temperature. A quantity of lactic acid bacteria (0.3 DCU/kg) was regenerated in a sterile bottle with 10 mL of plant-based juice 30 min before inoculation. The flasks were mixed regularly and gently to rehydrate the LAB strains and thus obtain a homogeneous solution without residue. The pH changes were recorded by a D291 consort program (Consort bvba, Turnhout, Belgium), with a 30 min interval between two measurements over a total duration of 6 h.

The pH and acidity measurements were performed once the fermented beverages were processed and stored at 4 °C. The pH measurement was taken at room temperature with a pHmeter (Model: Mettler Toledo F20-Kit FiveEasy™ Benchtop). Acidity is the measure of the quantity of acid (mL NaOH 0.1 N/10 g, expressed in degree of acidity).

Characteristics of fermented desserts during storage

The stability of the fermented products was evaluated at 4 °C for 30 d. The samples were analyzed at 24 h, 15 d, and 30 d after manufacture.

Texture analysis

The analysis of the texture of the fermented desserts was carried out using a Lloyd traction/compression analyzer (AMETEK S.A.S, Elancourt, France) equipped with a 10 N cell force, by carrying out the TPA texture profile analysis procedure, according to the method of [START_REF] Wang | Effects of endogenous small molecular compounds on the rheological properties, texture and microstructure of soymilk coagulum : Removal of phytate using ultrafiltration[END_REF] with minor modifications.

Each sample was compressed twice to 40 % of the original height with a cylindrical probe (1 cm in diameter, 10 cm long) at a crosshead speed of 2 mm/s.

Water Holding Capacity (WHC)

Water holding capacity (WHC) is defined as the ability of the sample to contain all or a portion of its water. The WHC of the samples was determined by a centrifugation method described by Wang et al., (2017) with slight modifications. The fermented products (5 g) were centrifuged at 6368 g at 4 °C for 20 min. WHC (g/100 g) is defined according to the following Eq. ( 5).

100 ⁄ = !! "#$%&#'( !! )#* * 100 (5)

Color analysis

Color is measured by the standard method established by the "Commission Internationale de l'Eclairage" (CIE), based on the determination of tristimulus values (XYZ), from a threedimensional space, called CIELAB. This space is a non-linear transformation of the tristimulus values and each color is defined by the brightness value L* and the color channels a* (red-green) and b* (blue-yellow). The mixtures color was analyzed using a colorimeter (Konica Minolta CR-400 series, Roissy CDG Cedex) using a white ceramic plate as a reference for calibration.

The Total Color Difference (∆E76) is a measure of the visual difference between the color of each sample relative to the reference color in a CIELAB space. The values of delta (∆E76) are obtained according to the following Eq. ( 6).

Microbiological analysis

A single product was defined for the realization of these experiments. It was determined according to the optimal paths of the mixture design, by selecting the optimal parameters of the results obtained previously. The selected product corresponds to the E1 mixture of the previously prepared mixture design, which is composed in proportion of dry matter of: 60 % chickpeas, 20 % coconut and 20 % oat. A sample (1 mL) of the fermented mixture was mixed with 9 mL of sterilized tryptone salt broth (NaCl 0.85 % (w/v)). A series of appropriate dilutions of the sample was individually inoculated in triplicate on selective mediums. Man, Rogosa and Sharp agar (MRS) were used for enumeration of Lactobacillus delbrueckii subsp.

bulgaricus and M-17 agar for enumeration of Streptococcus thermophilus and then incubated at 37 °C for 48 and 72 h respectively according to the method of [START_REF] Li | Production of tofu by lactic acid bacteria isolated from naturally fermented soy whey and evaluation of its quality[END_REF].

After incubation, petri dishes with 30 to 300 colonies were counted and expressed as log10 colony forming units per milliliter (CFU/mL) of sample. The Total Aerobic Flora (Yeasts, Molds and Enterobacteriaceae) and Thermoresistant, (Spores, Thermophilic and Mesophilic Aerobic Flora) were determined on the media Petrifilm ™ 3M ™ and the Mesophilic Flora nonlactic in an agar gelysate medium (SFA medium for the enumeration of contaminating germs in dairy products) of the brand Grosseron (Grosseron, 44220 Couëron, France).

Statistical analysis

The statistical analyses were performed using NemrodW software (NemrodW ® 27 Traverse de la Malvina Bât. B, 13012 Marseille, France), which is a software package for designing and analyzing experimental and mixture design. All experiments were run at least three times and allowed for the determination of mean values and standard deviations. ANOVA tests were performed to compare the means as well as the correlation coefficients between the other parameters. The variances were considered to be significant at the threshold of 95 % (p < 0.05).

Results and discussion

Composition of the flours

The composition of the six different flours have been grouped in Table 2, the results obtained show that there was a difference in their composition. The dry matter content varies from 89.87±0.13 g/100 g in the oat flour to 96.35±0.25 g/100 g in the coconut flour. The variability of the composition is linked to the variation in the origin of the seeds that make up the flours.

The highest protein contents were obtained in the Fabaceae (lupin flour) 39.17±1,75 g/100 g.

The fat composition was variable from 2.50±0.40 g/100 g in rice to 17.0±0.40 g/100 g in coconut, and between 6 and 11 g/100 g in the other flours. In terms of carbohydrate content, the highest concentrations of starch were recorded in cereals (oats and rice) and pseudocereals (quinoa), a value that reaches 72.49 g/100 g in rice flour, in contrast to the contents of sugar, which are the lowest in rice and oat flour, the highest value 14.32 g/100 g was recorded in coconut. The overall composition corresponds to a certain extent to literature results [START_REF] Trinidad | Dietary fiber from coconut flour : A functional food[END_REF][START_REF] Nowak | Assessment of the nutritional composition of quinoa (Chenopodium quinoa Willd.)[END_REF]S. Wang & Zhu, 2016;[START_REF] Sultan | Moisture mediated effects of γ-irradiation on physicochemical, functional, and antioxidant properties of pigmented brown rice (Oryza sativa L.) flour[END_REF][START_REF] Zhang | Solvent retention capacity of oat flour : Relationship with oat β-glucan content and molecular weight[END_REF][START_REF] Olukomaiya | Effect of solid-state fermentation on proximate composition, anti-nutritional factor, microbiological and functional properties of lupin flour[END_REF] and producer claims for each flour type.

Physicochemical characteristics

The grinding process used in the production of flour can be different from one seed to another which can be the source of size differences in studied flours [START_REF] Benković | Flow properties and chemical composition of carob (Ceratonia siliqua L.) flours as related to particle size and seed presence[END_REF]. Analyses have shown a significant effect of the type of flour on the particle size (Table 2). The particle diameter varied from 11.27±0.06 to 971.67±58.24 μm, with a median diameter (D50) between 36.63±1.40 µm in chickpea and Liu (2016) for rice flour. Therefore, the ΔHGel values recorded are lower in all the flours studied compared to these publications which are between 3 and 6 J/g. This difference may be due to the flour/water ratio used. According to [START_REF] Zhu | Structures, properties, modifications, and uses of oat starch[END_REF], amylopectin content affects the reduction of gelatinization parameters as well as certain grinding parameters such as time.

The water absorption index (WAI), the water solubility (WS) and the swelling power (SP) characterize the interaction of the product with water and are often important indicators that affect the subsequent treatment food. The results of WAI, WS, and SP at 30 and 90 °C are shown in Fig. 1. An increase in water absorption was observed in all flours with an increase in hydration temperature, followed by an increase in swelling power and a decrease in solubility.

The results were partly in agreement with the literature [START_REF] Kraithong | Physicochemical and functional properties of Thai organic rice flour[END_REF][START_REF] Lai | Effects of hydrothermal treatment on the physicochemical properties of pregelatinized rice four[END_REF][START_REF] Maqbool | Influence of soaking and germination on physico-chemical composition and functional properties of chickpea (var. SKUAST-233) flour[END_REF], which reported WAI values of 2.38-2.75 g/g, WSI values of 9.84-20.10 g/100 g and SP values of 2.15-2.42 in rice and chickpea flours.

According to [START_REF] Kraithong | Physicochemical and functional properties of Thai organic rice flour[END_REF] a high WAI can be explained by the presence of a large number of hydrophilic groups in starch molecules, high sugar and protein content. In the case of chickpea flour, a high absorption index, greater than 7 (g/g), was observed as the presence of significant levels of starch and protein. The latter may have promoted a strong involvement of hydrogen bonds with polar or charged side chains. Also, a large particle size contributes to low absorption, just as a high lipid content would block hydration with hydrophobic parts as it can be observed for lupin flour characterized by the lowest WAI at 90 °C, with a high lipid content and a very low starch content.

Fermentability of flours

The evolution of acidification in the six flours studied was compared in Table 2. A decrease in pH was observed after the addition of the LAB strains, resulting an acidification of the medium with a final pH ranging from 4.81±0.06 to 4.40±0.07. The decrease in pH could be attributed to the increase in organic acids produced by microbial activities during the degradation of carbohydrates [START_REF] Olukomaiya | Effect of solid-state fermentation on proximate composition, anti-nutritional factor, microbiological and functional properties of lupin flour[END_REF]. Although, the use of quinoa flour, chickpea and coconut allows a reduction of the pH, in less than 16 h of incubation, oat, rice, and lupin flours have a prolonged latency phase, with pH above 4.5, even after 21 h of incubation. This delay could be attributed to the slow metabolism of the strains present in the LAB strains, the nature of the substrates consumed, or the incubation temperature [START_REF] Hickisch | Influence of lupin-based milk alternative heat treatment and exopolysaccharide-producing lactic acid bacteria on the physical characteristics of lupin-based yogurt alternatives[END_REF]. A variation which is also found on the average speeds which is at least of the order of -0.05±0.004 pH unit/hour in lupin flour, and which can expect at most -0.29±0.01 pH unit/hour in quinoa flour.

Mixture design

Choice of ingredients

The previous step specified above allowed the selection of three compounds for the formulation of a fermented vegetable product, so the choice is to focus on oat, chickpea, and coconut flour. The combination of these flours will give the products a firm texture without the addition of gelling agents because of their starch composition and their water absorption properties during gelling in the presence of water at high temperatures. The combination of these flours will allow fermentation and bring nutritional value to the final products with protein, fat, and fiber content. Oats and chickpeas are pregnant foods and among the most consumed cereals and legumes in the world for their rich nutritional intake. They provide carbohydrates, mainly in the form of starch, soluble dietary fiber, lipids, balanced proteins, and several vitamins. Coconut flour is a real economic interest as a residue generated during the extraction of coconut milk and coconut oil. This flour is not only abundant in dietary fiber but also trans fatty acids and low in carbohydrates. The consumption of oat, chickpea or coconut products has been associated with various health impacts such as the reduction of cholesterol, the risk of cardiovascular disease, gastrointestinal disorders, as well as cancer prevention [START_REF] Trinidad | Dietary fiber from coconut flour : A functional food[END_REF][START_REF] Sreerama | Nutrients and antinutrients in cowpea and horse gram flours in comparison to chickpea flour : Evaluation of their flour functionality[END_REF][START_REF] Martínez-Villaluenga | Health benefits of oat : Current evidence and molecular mechanisms[END_REF].

Acidification kinetics

Fermentation kinetics were carried out to determine the pH evolution during a 6 h incubation at 45 °C. The results are shown in Fig. 2.

The pH varies from 6.48±0.12 to 6.08±0.03 after 30 min from the start of fermentation, and from 4.76±0.04 to 4.58±0.06 after 6 h of fermentation depending on the formulation. In the case of fermentation of flour mixtures (oats, chickpeas, coconut), acidification is more accelerated, unlike fermentation of flours alone. Mixing would have reduced the time required to reach a target pH close to 4.5. This may be related to the composition of the raw materials (presence of different sugars), which are consumed by the LAB strains, resulting in a reduction of the pH of the medium [START_REF] Da Silva Fernandes | Evaluation of the isoflavone and total phenolic contents of kefirfermented soymilk storage and after the in vitro digestive system simulation[END_REF]. The presence of legumes that are a source of complex carbohydrates (fiber, starch and oligosaccharides), which also contain growth factors and prebiotic components for the bacteria, also promotes a gentle fermentation process [START_REF] Akin | Functional properties of fermented milk produced with plant proteins[END_REF]. The pH decreased after 24 h of storage, ranging from 4.42±0.01 to 4.56±0.00 and the acidity from 3.73±0.06 to 4.80±0.20, to stabilize during the month of storage at 4 °C. According to [START_REF] Gupta | Development of fermented oat flour beverage as a potential probiotic vehicle[END_REF], a very small variation in pH of fermented oat flour beverage products stored under refrigerated conditions can be recorded, attributed to the low metabolic activity of the LAB strains, and therefore to the low acid production during the storage period. By reference to the type of compound and their interactions, the pH increases and the acidity decreases with the addition of oats. An effect that is also important during the interaction of the 3 ingredients (oat, chickpea and coconut).

However, the interactions among the flours in pairs had negative effects on the pH response and positive effects on the acidity response.

Influence of flour composition on color difference

According to the results obtained, there is a significant effect (p<0.001) of the flours used on the color of the mixtures. The results (Table 3) indicate that product E1 had the highest ΔE (17.37±0.01), the furthest difference from the reference white, while E11 was the product with the lowest ΔE (14.92±0.01). The color differences obtained in the mixtures (oats, chickpea, coconut) are lower than those reported by [START_REF] Luana | Manufacture and characterization of a yogurt-like beverage made with oat flakes fermented by selected lactic acid bacteria[END_REF], (ΔE>18) in the fermented oat flake products. This shows the interaction between the three flours used and their impact in the reduction of ΔE. Consequently, the final color of the product can be predicted according to the composition using the response surfaces presented in Fig. 3.

The ΔE did not change during the month of storage, the largest variation is less than 1 unit, which is not detectable with the naked eye. Analysis after 24 h showed that the increase in oat content plays a role in the final color of the product. This influence is half as important for coconuts and up to four times less for chickpeas. The oat-coconut interaction had a negative effect which is most significant on the color difference.

Influence of flour composition on water holding capacity

After fermentation, the water retention rate varies between 58.03±0.21 g/100 g and 75.31±0.72 g/100 g in the products obtained. Statistical analysis shows that these values increase significantly with the increase in oat and chickpea flours. After 30 d of storage at 4 °C, the WHC was between 56.64±0.10 g/100 g and 66.12±0.09 g/100 g, reflecting an increase in syneresis. The result is a small water pill on the surface of the jars in products where water retention is less than 65 g/100 g, a phenomenon also common in dairy yogurt. According to Shi [START_REF] Wang | Evaluation of chickpea as alternative to soy in plant-based beverages, fresh and fermented[END_REF], a fermented plant product containing starch may undergo retrogradation during cooling, causing water to be expelled from the gel, resulting in syneresis. The WHC in mixtures can decrease by up to 13 g/100 g during 30 d storage at 4 °C.

These results are close to those obtained by [START_REF] Coda | Yogurt-Like Beverages Made With Cereals[END_REF], which indicate a decrease in WHC of 10 g/100 g under the same conditions and for the same storage time.

Influence of flour composition on texture

There was a significant effect of composition on the hardness and cohesion of processed products (Tables 2 and3). The study by [START_REF] Hickisch | Influence of lupin-based milk alternative heat treatment and exopolysaccharide-producing lactic acid bacteria on the physical characteristics of lupin-based yogurt alternatives[END_REF], on alternatives to yogurt showed that heat treatment and the choice of fermentation micro-organism had an impact on the physical characteristics of the products. It was observed that at 24 h of storage, hardness varies from 0.10±0.00 N to 0.23±0.01 N and cohesion from 0.31±0.10 to 0.69±0.03. After 30 d of storage, the hardness increases slightly from 0.02 to 0.12 N, and the cohesion decreases from 0.01 to 0.22.

The addition of oats resulted in a firmer gel due to the ability of this flour to retain water. On the other hand, the cohesion of the gel increased considerably with the oat-coconut interaction. The firmness of the yogurt is directly related to its total solids content, particularly its protein content and type (Mishra [START_REF] Pandey | Optimization of the prebiotic & probiotic concentration and incubation temperature for the preparation of synbiotic soy yoghurt using response surface methodology[END_REF]. A higher content of different proteins, but also the effect of pH would lead to a higher degree of cross-linking of the gel network, which would result in an increase in structural strength, potentially leading to an increase in water holding capacity, firmness and viscosity of the gels (Mishra [START_REF] Pandey | Optimization of the prebiotic & probiotic concentration and incubation temperature for the preparation of synbiotic soy yoghurt using response surface methodology[END_REF][START_REF] Akin | Functional properties of fermented milk produced with plant proteins[END_REF][START_REF] Wang | Evaluation of chickpea as alternative to soy in plant-based beverages, fresh and fermented[END_REF].

Microbiological quality and viability of lactic strains

The heat treatment applied on the mixture before fermentation, confirmed the microbiological quality during storage. Indeed, an acidification kinetics of the E1 mixture during fermentation (shown in fig. 4) compared to a control mixture without LAB strains, which didn't significantly acidify. Through the first week of storing at 4 °C, there was no microbial contamination. The load is as low, after 30 days, except for mesophilic non-lactic flora, moulds and yeast, which was 41, 11 and 10 germs/mL respectively, a very low microbial charge. [START_REF] Olukomaiya | Effect of solid-state fermentation on proximate composition, anti-nutritional factor, microbiological and functional properties of lupin flour[END_REF] had indicated that, the number and increase of bacteria, and more particularly yeasts and moulds, can be attributed to variable reactions and changes in the pH and other biochemical parameters in the fermented product.

According to [START_REF] Jeske | Past, present and future : The strength of plantbased dairy substitutes based on gluten-free raw materials[END_REF] the enumeration of Streptococcus thermophilus and Lactobacillus bulgaricus used in the production of yoghurt from cow's milk reaches a population of 20×10 8 cells/mL, which corresponds to about 9.30 log10 CFU/mL. The enumeration of lactic acid bacteria strains was carried out during four weeks of storage (Table 5). The quantity of both strains of bacteria at 24 h after the end of fermentation was 8.04±0.04 log10 CFU/mL for L. delbrueckii subsp. Bulgaricus and 8.10±0.06 log10 CFU/mL for S.

thermophilus. These values were stable during one month of storage at 4 °C, with a slight decrease in L. delbrueckii subsp. Bulgaricus count to 7.92±0.05 log10 CFU/mL. Considering that a probiotic must be a fermented product containing a minimum number of 6 log10 CFU/mL microorganisms throughout its shelf life [START_REF] Da Silva Fernandes | Evaluation of the isoflavone and total phenolic contents of kefirfermented soymilk storage and after the in vitro digestive system simulation[END_REF].

Conclusion

The aim of the present study was to highlight the influence of a given concentration of three ingredients individually as well as their interaction on the physicochemical and textural properties of fermented vegetable desserts. It was shown that there is a significant variation between the different flours studied, especially for composition such as starch content, which has a direct impact on properties such as gelatinization temperatures and water absorption. The results showed that the choice of the applied mixture plays a decisive role on the acidification of the mixture, the color, the water retention and the texture of the final products.

Thus, acidification is faster in a mixture than individual flour fermentation. The response surface showed a significant effect of oat-coco interaction on color, and more specifically the effect of oat concentration on water retention which impacts firmness and gel stability during storage.

The next step of this study is the selection and optimization of an ideal product, combining the information obtained, in order to provide the best solutions for the formulation of new vegetable alternatives to dairy products.

  the (ref) index indicates the sample values used as a reference and the (s) index for the sample values.

  313.33±4.73 µm in oat. The distribution was monomodal in coconut and rice flours, bimodal for oat and lupin, and multimodal for quinoa and chickpea flours.[START_REF] González-Montemayor | Prosopis spp. Functional activities and its applications in bakery products[END_REF] have reported that particle size affects the physico-chemical properties of flours, increasing the surface area per unit volume, and the distribution of macronutrients and functional components is expected to be different.The DSC analysis provided the starch gelatinization parameters. According to[START_REF] Ndangui | Production et caractérisation de farine de patate douce (Ipomoeabatatas.Lam) : Optimisation de la technologie de panification[END_REF], the presence of a single endothermic peak during DSC analysis is typical of starch gelatinization in the presence of excess water. This continuous endothermic transition indicates a swelling of the granule and a progressive disruption of the starch crystal lattice in the gelatinization temperature range. This phenomenon was not observed on flours without starch (i.e. lupin and coconut), unlike oat, quinoa, chickpea and rice flours. Enthalpy, which corresponds to the internal energy of the system resulting from the heat treatment, is proportional to the starch content. The lowest value is recorded with chickpea (0.80±0.09 J/g) which has the lowest starch content, followed by oats (1.36±0.15 J/g), quinoa (1.99±0.10 J/g) and the highest value with rice (2.04±0.02 J/g) which has the highest starch content. The results of characteristic temperatures (T0, Tp, and Tf) obtained correspond to what is reported in the literature by[START_REF] Kaur | Studies on functional, thermal and pasting properties of flours from different chickpea (Cicer arietinum L.) cultivars[END_REF] for chickpea flour and Ye,Wang, Wang, Zhou, &

  The individual fermentation tests of the flours with the LAB strains were positive with a decrease in pH down to 4.40. The comparison of the composition, hydration parameters, starch gelatinization and acidification kinetics allowed the choice of three ingredients: oat, chickpea and coconut flours. An optimization in term of physico-chemical term was then done by mixture design with constraints.

Fig

  Fig. 1. Hydration properties, water absorption index (WAI), water solubility index (WSI) and swelling power (SP) of flours.

Fig. 3 .Fig. 4 .

 34 Fig. 3. (a) 3D and (b) 2D contour plot Variation of responses, color difference (ΔE), water holding capacity (WHC) and texture (hardness, cohesion) depending on the amount of flour at 15 d of storage. (Color). A: Oat, B: Coconut, C: Chickpea

Table 1

 1 Mixture design of 13 formulations based on 3 compounds with constraints.

	No.exp.	Name	Chickpea	Coconut	Oat
	1	E1	0.6000	0.2000	0.2000
	2	E2	0.5000	0.4000	0.1000
	3	E3	0.4000	0.4000	0.2000
	4	E4	0.6000	0.3000	0.1000
	5	E5	0.5000	0.3000	0.2000
	6	E6	0.6000	0.2500	0.1500
	7	E7	0.4500	0.4000	0.1500
	8	E8	0.5500	0.3500	0.1000
	9	E9	0.5250	0.3250	0.1500
	10	E10	0.5625	0.2625	0.1750
	11	E11	0.5125	0.3625	0.1250
	12	E12	0.4625	0.3625	0.1750
	13	E13	0.5625	0.3125	0.1250

The mixture results were then analyzed by NemrodW software with a cubic model.

Table 3

 3 Color difference (ΔE), water holding capacity (WHC) and texture (hardness, cohesion) of mixture design products at different storage times.

		Storage time ΔE	WHC (g/100 g) Hardness (N) Cohesion
	E1	h	17.37±0.01	70.62±0.43	0.23±0.01	0.51±0.06
		d	16.76±0.00	67.05±0.29	0.26±0.03	0.41±0.06
		d	17.37±0.00	65.92±0.32	0.30±0.01	0.34±0.11
	E2	h	15.72±0.00	58.03±0.21	0.10±0.00	0.69±0.03
		d	15.26±0.01	54.57±1.25	0.12±0.00	0.37±0.08
		d	15.67±0.05	57.20±0.03	0.14±0.00	0.56±0.04
	E3	h	16.24±0.01	75.31±0.72	0.14±0.01	0.41±0.03
		d	15.63±0.00	63.93±0.34	0.11±0.01	0.56±0.07
		d	16.11±0.01	63.02±0.87	0.14±0.01	0.35±0.10
	E4	h	15.62±0.00	68.07±0.39	0.11±0.00	0.31±0.10
		d	15.35±0.02	58.46±1.19	0.19±0.01	0.42±0.08
		d	15.76±0.01	56.64±0.10	0.15±0.01	0.48±0.04
	E5	h	16.29±0.01	74.28±0.22	0.14±0.01	0.53±0.06
		d	15.84±0.00	64.87±0.57	0.11±0.01	0.52±0.07
		d	16.52±0.01	60.97±0.27	0.12±0.01	0.58±0.02
	E6	h	16.72±0.01	73.02±0.40	0.22±0.01	0.48±0.06
		d	16.53±0.07	65.33±0.31	0.26±0.01	0.44±0.03
		d	16.78±0.01	66.12±0.09	0.31±0.02	0.34±0.08
	E7	h	15.57±0.01	66.66±0.63	0.10±0.00	0.59±0.18
		d	15.30±0.01	59.09±0.91	0.09±0.00	0.40±0.08
		d	15.47±0.01	58.02±0.72	0.16±0.03	0.38±0.04
	E8	h	15.40±0.00	64.73±0.26	0.10±0.00	0.49±0.17
		d	15.29±0.01	58.11±0.50	0.15±0.01	0.42±0.12
		d	15.46±0.00	61.18±1.83	0.16±0.01	0.47±0.06
	E9	h	15.95±0.00	69.40±0.30	0.14±0.01	0.48±0.04
		d	16.19±0.01	63.96±0.55	0.17±0.01	0.53±0.02
		d	16.49±0.01	62.72±0.22	0.22±0.01	0.46±0.03
	E10 h	16.03±0.01	72.75±0.48	0.16±0.01	0.43±0.06
		d	16.31±0.01	65.24±0.45	0.20±0.02	0.48±0.02
		d	16.72±0.01	64.64±0.65	0.29±0.02	0.42±0.06
	E11 h	14.92±0.01	66.15±0.36	0.11±0.01	0.50±0.05
		d	15.19±0.01	61.97±0.38	0.16±0.00	0.40±0.09
		d	15.54±0.00	61.93±0.95	0.17±0.01	0.49±0.04
	E12 h	15.64±0.00	67.90±0.53	0.13±0.01	0.44±0.06
		d	15.89±0.00	62.29±0.59	0.17±0.01	0.45±0.07
		d	15.89±0.00	65.04±0.63	0.17±0.00	0.44±0.01
	E13 h	15.52±0.00	65.17±0.28	0.13±0.01	0.54±0.03
		d	15.85±0.00	61.32±0.42	0.12±0.00	0.48±0.06
		d	15.92±0.01	60.47±0.42	0.23±0.12	0.40±0.11
	All analyses were performed in triplicate.

Table 4 .

 4 Estimation of the coefficients of the mathematical model for the answers: Color difference (ΔE), water retention capacity (WHC) and texture (hardness, cohesion) at different storage times.

		Storage time ΔE	WHC (g/100 g) Hardness (N) Cohesion
	Variance analysis	24 h 15 d 30 d	<0.01*** <0.01*** <0.01***	<0.01*** <0.01*** 0.04***	<0.01*** <0.01*** <0.01***	3.32* 0.49** <0.01***
	b1	24 h 15 d	54.34** 0.37	174.71 -206.03**	0.10 -0.01	-5.26 -2.50
		30 d	16.00	-349.05*	-4.04	8.68**
	b2	24 h 15 d	102.07* -20.28	29.94 -641.09***	-1.03 -1.10	-17.39 -11.75
		30 d	14.71	-866.42*	-9.52	20.74**
	b3	24 h 15 d	230.58 -153.23	369.79 -1855.05***	-12.62 -13.87	-39.30 -11.29
		30 d	5.02	-2857.33**	-49.85*	91.74***
	b12	24 h 15 d	-240.66* 85.21	-282.78 1756.12***	1.13 1.71	44.62 28.35
		30 d	-4.62	2500.65**	24.50	-52.39**
	b13	24 h 15 d	-416.47 315.73	-927.02 3832.90***	20.48 22.88	82.41 27.68
		30 d	39.55	5908.72**	90.75*	-171.49***
	b23	24 h 15 d	-632.67 414.73	-144.77 6001.88***	29.61 30.71	147.03 75.95
		30 d	23.96	8693.57**	123.11	-260.44***
	b123	24 h 15 d	935.01 -558.82	1569.31 -9196.49***	-36.40 -38.56	-253.91 -133.60
		30 d	-51.90	-14004.66**	-172.13	410.18***
	R 2	24 h 15 d	0.87 0.88	0.91 0.95	0.86 0.74	0.42 0.60
		30 d	0.91	0.70	0.76	0.87
	b1: Chickpea; b2: Coconut; b3: Oat; *: significant effect (p<0.05); **: significant effect (p<0.01), ***: significant
	effect (p<0.001); R 2 : the coefficient of determination.

Table 5 :

 5 Control of the viability of the lactic acid bacteria strains in the optimum product during one month of storage

	Storage time	Viable counts of strains [log10 CFU/mL]
	L. delbrueckii subsp. Bulgaricus	S. thermophiles
	24 h	8.04±0.04	8.10±0.06
	7 d	8.09±0.06	7.99±0.02
	15 d	8.05±0.07	8.21±0.05
	30 d	7.92±0.05	8.10±0.05
	All analyses were performed in triplicate

Acknowledgements:

We acknowledge Saint Hubert and Université de Lorraine for their support in the realization of this study and Mrs AIT CHEKDID's thesis.

Author contributions:

A.C.A., L.M., designed the study. A.C.A., carried out the experiment, wrote the manuscript with support from K.C.J.F., and L.M. K.C.J.F., provided critical feedback and helped shape the research, analysis and manuscript. L.M., designed the model and the computational framework, analysed the data and supervised the work. P.E., F.M., L.B., C.C., helped supervise the project. All authors final approval of the version to be submitted.


Sasinee Pruekprasert
email: sasinee@nii.ac.jp

Toru Takisaka
email: takisaka@nii.ac.jp

Clovis Eberhart
email: eberhart@nii.ac.jp

Ahmet Cetinkaya
email: cetinkaya@nii.ac.jp

Jérémy Dubut
email: dubut@nii.ac.jp

Moment Propagation of Discrete-Time Stochastic Polynomial Systems using Truncated Carleman Linearization

Keywords: Stochastic systems, nonlinear systems, probabilistic safety analysis, moment computation, Carleman linearization

We propose a method to compute an approximation of the moments of a discretetime stochastic polynomial system. We use the Carleman linearization technique to transform this finite-dimensional polynomial system into an infinite-dimensional linear one. After taking expectation and truncating the induced deterministic dynamics, we obtain a finite-dimensional linear deterministic system, which we then use to iteratively compute approximations of the moments of the original polynomial system at different time steps. We provide upper bounds on the approximation error for each moment and show that, for large enough truncation limits, the proposed method precisely computes moments for sufficiently small degrees and numbers of time steps. We use our proposed method for safety analysis to compute bounds on the probability of the system state being outside a given safety region. Finally, we illustrate our results on two concrete examples, a stochastic logistic map and a vehicle dynamics under stochastic disturbance.

INTRODUCTION

Ensuring safety in cyber-physical systems under uncertainties is an important and challenging task. For instance, automated driving systems are required to take into account the uncertainties in the environment and those in control operations to avoid collision with obstacles [START_REF] Bry | Rapidly-exploring random belief trees for motion planning under uncertainty[END_REF]. Motion planning algorithms in automated driving typically compute probabilities of collision in different potential paths of a vehicle by identifying possible future states of the vehicle under the effects of measurement and actuation uncertainties.

To quantify how uncertainties affect the state of the system over time, one often uses the idea of mean and covariance propagation. Particularly, if the dynamical system is linear and subject to additive Gaussian noise, then this state follows a Gaussian distribution, and the evolution of its mean and covariance can be described by a linear equation. Given the initial mean and covariance, one can then use this linear equation to compute mean and covariance iteratively at each future time instant. This idea has been used in conjunction with a Kalman-filtering approach in motion planning by [START_REF] Bry | Rapidly-exploring random belief trees for motion planning under uncertainty[END_REF] and [START_REF] Banzhaf | From footprints to beliefprints: Motion planning under uncertainty for maneuvering automated vehicles in dense scenarios[END_REF]. Here, we are interested in whether similar iterative techniques can be used to compute mean, covariance, and even higher moments when the dynamical system is

The authors are supported by ERATO HASUO Metamathematics for Systems Design Project No. JPMJER1603, JST; J. Dubut is also supported by Grant-in-aid No. 19K20215, JSPS. nonlinear and the noise is not necessarily additive and Gaussian. We provide a positive answer to this question.

In this paper, we propose a method to compute approximations of the moments through truncated Carleman linearization [START_REF] Steeb | Non-linear autonomous systems of differential equations and Carleman linearization procedure[END_REF]. Specifically, we consider a discrete-time nonlinear stochastic polynomial system with a model that allows us to describe additive and multiplicative stochastic noise in a unified manner through randomly varying coefficient matrices. We use the Carleman linearization technique on this system to obtain an infinite-dimensional linear one that characterizes the evolution of all Kronecker powers of the state variables. We observe that the expected dynamics of this linear system relates the moments of the random coefficient matrices to the moments of the state. By using a truncation approach, we then obtain a finite-dimensional linear system, which we use to compute an approximation of the mean and covariance of the system state and its higher moments up to a given truncation limit.

We then propose a framework to verify the stochastic safety of the system via tail probability estimation techniques [START_REF] Gray | A general method for approximating tail probabilities[END_REF]. This framework relies on the efficient computation of error bounds for our moment approximations, which we provide. The exact error is the sum of many terms whose number grows quickly with respect to the dimension of the system and time, and we give an efficient way to compute a bound by splitting the sum into two parts: the part (e. size in the dimension of the system) of the sum that strongly contribute to the error and on which we make a fine computation, and the rest on which we do a cruder but faster computation.

We demonstrate the utility of our approach with two case studies. First, we compute moments for the scalar stochastic logistic map with uniformly distributed random growth/decay rates. We then look at a higherdimensional polynomial discrete-time vehicle model obtained via second-order Taylor expansion of a kinematic bicycle model. We use our approach to investigate how the uncertainty in the acceleration affects the vehicle's future positions. In particular, our approach uses the moments of the initial state of the vehicle, together with the moments of its noisy acceleration. The moments of the initial state of the vehicle is obtained by using a noisy measurement of the state and the information on the statistical properties of the measurement error. Our framework allows moment computation in a receding horizon fashion. Specifically, at each time step, new measurements can be used for specifying the moments of the new "initial" state, which can then be used to compute the moments associated with time instants that are further in the future.

Carleman linearization is a well-known approach in nonlinear dynamical systems literature. It has been used in several pieces of work for approximation [START_REF] Bellman | On some questions arising in the approximate solution of nonlinear differential equations[END_REF][START_REF] Steeb | Non-linear autonomous systems of differential equations and Carleman linearization procedure[END_REF][START_REF] Al-Tuwaim | Discretization of nonlinear models using a modified Carleman linearization technique[END_REF], and more recently for optimal control [START_REF] Amini | Carleman state feedback control design of a class of nonlinear control systems[END_REF] and model predictive control of continuoustime deterministic nonlinear systems [START_REF] Hashemian | Feedback control design using model predictive control formulation and Carleman approximation method[END_REF]. Results on Carleman linearization for stochastic systems are relatively scarce. For Ito-type stochastic differential equations with bilinear noise, [START_REF] Wong | Carleman linearization and moment equations of nonlinear stochastic equations[END_REF] used Carleman linearization to obtain differential equations describing the evolution of the moments. Moreover, [START_REF] Rauh | Carleman linearization for control and for state and disturbance estimation of nonlinear dynamical processes[END_REF] used a series-expansion approach together with the Carleman linearization technique to approximate a continuous-time nonlinear stochastic system with a discrete-time linear system under additive noise. They then use the linear system for state estimation via a Kalman filter. Recently, [START_REF] Cacace | A Carleman discretization approach to filter nonlinear stochastic systems with sampled measurements[END_REF][START_REF] Cacace | Optimal linear filter for a class of nonlinear stochastic differential systems with discrete measurements[END_REF] developed Carleman linearization-based filters for nonlinear continuous-time stochastic systems with Wiener noise and periodic measurements. Although the discretetime representation of the sampled-data filter equations in those pieces of work involves multiplicative noise, the models and the problem formulations there differ from those in our paper: on the one hand, they can describe non-polynomial dynamics, but on the other hand, they can only handle Wiener noises. We note that our truncation error analysis is partly motivated by the analysis in [START_REF] Forets | Explicit error bounds for Carleman linearization[END_REF], which provides tight approximation error bounds for the truncated Carleman linearization approach in deterministic continuous-time systems. We also note that in addition to Carleman linearization, Koopman operators are also used for deriving infinite-dimensional linear equations to describe finite-dimensional nonlinear systems [START_REF] Goswami | Global bilinearization and controllability of control-affine nonlinear systems: A Koopman spectral approach[END_REF]; [START_REF] Mesbahi | On modal properties of the Koopman operator for nonlinear systems with symmetry[END_REF], although the techniques for analysis of these two methods are quite different.

The rest of the paper is organized as follows (see also Figure 1 for the specific subsections where the main steps of the method are presented). In Section 2, we start from a finite-dimensional stochastic polynomial system, on which we apply the Carleman linearization technique to get an infinite-dimensional linear stochastic system, which we then turn into a deterministic one by taking expectation. We present our truncation-based moment approximation method and study the error it introduces in Section 3, in which we also give an application of this bound. We provide two numerical examples to demonstrate our approach in Section 4. Finally, in Section 5, we conclude the paper.

Notations. We use N for the set of non-negative integers.

We denote by M 1 ⊗ M 2 the Kronecker product of two matrices M 1 and M 2 . Moreover, we use M [k] to denote the kth Kronecker power of the matrix M , given by M

[0] = 1 and M [k] = M [k-1] ⊗ M for k > 0.
The dimensions of the matrices we use grow rapidly, so we will use n ≤k as a shorthand for k i=0 n i when manipulating dimensions of vectors or matrices.

CARLEMAN LINEARIZATION FOR STOCHASTIC POLYNOMIAL SYSTEMS

In this section, we use the Carleman linearization technique to turn a finite-dimensional discrete-time stochastic polynomial system into an infinite-dimensional one, then take expectation to get a deterministic system that describes the evolution of all moments of the system state.

Discrete-Time Stochastic Polynomial Systems

We focus on computing the moments of the state of the following finite-dimensional discrete-time stochastic polynomial system

x(t + 1) = dS i=0 F i (t)x [i] (t), t ∈ N, x(0) = x 0 , (1)
where x(t) ∈ R n is the state vector and

F i (t) ∈ R n×n i , i ∈ {0, 1, . . . , d S }
, are randomly-varying coefficient matrices.

To be more precise, we assume given a probabilistic space Ω, then x 0 (resp. x(t), F i (t)) is actually a measurable function from Ω to R n (resp. R n , R n×n i ). We will omit to mention Ω in this paper.

We remark that (1) is useful for modeling dynamics with both additive and multiplicative noise terms: the vector F 0 (t) ∈ R n represents additive noise, as x [0] = 1, while the terms F 1 (t), . . . , F dS (t) characterize the effects of multiplicative noise.

Carleman Linearization

We use Carleman linearization to obtain an infinitedimensional linear system that describes the evolution of the Kronecker powers of the state x(t). By defining

F (t) [F 0 (t) F 1 (t) F 2 (t) • • • F dS (t)] ∈ R n×n ≤d S and y(k, t) 1 x(t) x [2] (t) • • • x [k] (t) ∈ R n ≤k , (2)
the dynamical system (1) can be rewritten as

x(t + 1) = F (t)y(d S , t).
Therefore, for all j ∈ N, we have

x [j] (t + 1) = (F (t)y(d S , t)) [j] .

By using the mixed-product property of the Kronecker product (i.e., (A⊗B)(C ⊗D) = AC ⊗BD; see Section 13.2 of Laub ( 2005)), we obtain

x [j] (t + 1) = F [j] (t)y [j] (d S , t) = jdS k=0   (i l ) l≤j ∈H j,k F i1 (t) ⊗ • • • ⊗ F ij (t)   x [k] (t),
for j ∈ N, where

H j,k (i l ) l≤j j l=1 i l = k and i l ≤ d S .
We thus get the infinite-dimensional linear system

x [j] (t + 1) = jdS k=0   (i l ) l≤j ∈H j,k j m=1 F im (t)   x [k] (t), x [j] (0) = x [j] 0 , (3)
which describes the evolution of all Kronecker powers x [0] (t), x [1] (t), . . . of the state x(t).

In order to give a simpler description of the system, we introduce the matrix A j,k (t) ∈ R n j ×n k given by

A j,k (t) (i l ) l≤j ∈H j,k F i1 (t) ⊗ • • • ⊗ F ij (t).
Note in particular that A 0,0 (t) = 1 since there is only the empty sequence in H 0,0 , and A j,k (t) = 0 if j = 0 and k > 0 or k > jd S since H j,k is then empty. We also introduce, for all N, M ∈ N, the matrix A(t; N, M ) defined by blocks as

A(t; N, M ) = [A j,k (t)] j≤N k≤M , that is, A(t; N, M )    A 0,0 (t) • • • A 0,M (t) . . . . . . . . . A N,0 (t) • • • A N,M (t)    ∈ R n ≤N ×n ≤M .
Then, for any k ∈ N, we have

y(k, t + 1) = 1 x(t + 1) x [2] (t + 1) • • • x [k] (t + 1) = A(t; k, kd S ) 1 x(t) x [2] (t) • • • x [kdS] (t) = A(t; k, kd S )y(kd S , t).
(4)

Moment Equations

We now derive the deterministic system that describes the evolution of the moments of x(t) by taking expectation in (4). This gives

E[y(k, t + 1)] = E[A(t; k, kd S )y(kd S , t)].
We make the following assumptions concerning the coefficient matrices and the random initial state x 0 . Assumption 1. The matrices {F (t) | t ∈ N} are independent. Moreover, x 0 is independent of F (t), t ∈ N. Assumption 2. The matrices {F (t) | t ∈ N} are identically distributed.

Note that Assumptions 1 and 2 are not overly restrictive and they hold in fairly general situations as we discuss in Section 4. Notice also that under Assumption 1, matrices F i (t) and F j (t) are still allowed to statistically depend on each other. Moreover, Assumption 2 allows us to obtain a "time-invariant" method to compute moments, further yielding computational advantage.

By iteration of (4), we get that y(kd S , t) is given by

y(kd S , t) = A(t -1; kd S , kd 2 S ) • • • A(0; kd t S , kd t+1 S ) • y(kd t+1 S , 0), t ∈ N. (5)
It follows from Assumption 1 that A(t; k, kd S ) and y(kd S , t) in ( 4) are mutually independent. To see this, observe that A(t; k, kd S ) is composed of the matrices F i (t), which are independent of x 0 and F (t -1), . . . , F (0), which determine y(kd S , t) as given by ( 5).

It then follows that

E[y(k, t + 1)] = E[A(t; k, kd S )]E[y(kd S , t)].
Here again, to give a simpler description of the system, we introduce new matrices. Notice that the coefficients of the matrix E[A(t; N, M )] are products of the moments of coefficients of F (t) and thus independent of t by Assumption 2.

We denote this matrix by E(N, M ) ∈ R n ≤N ×n ≤M , emphasizing the fact that it is independent of the time. Similarly, we denote by E i,j ∈ R n i ×n j the matrix E[A i,j (t)]. The equation above can thus be rewritten to give the following system

E[x [j] (t + 1)] = jdS k=0 E j,k E[x [k] (t)], t ∈ N, E[x [j] (0)] = E[x [j] 0 ]. (6)
Since the matrices E j,k only depend on the moments of the matrices F (•), they can be computed offline.

MOMENT APPROXIMATION THROUGH TRUNCATION

We will now introduce an approach that allows us to compute approximations of the moments of x(t). This truncation approach is critical, as an exact computation of the moments is impossible from a practical point of view. Indeed, by (2), computing the first k moments at time t amounts to computing E[y(k, t)]. So, by iteration of ( 6),

E[y(k, t)] = E(k, kd S ) • • • E(kd t-1 S , kd t S )E[y(kd t S , 0)].
As this indicates, exact computation of the first k moments requires the knowledge of matrices E(k, kd S ), E(kd S , kd 2 S ), . . ., E(kd t-1 S , kd t S ) of exponentially increasing dimensions, making any practical computation unrealistic. We thus resort to a truncation approach where we fix a truncation limit and consider a matrix of fixed size to approximately compute the moments of x(t).

Approximate Moments and the Truncated System

In this section, we define the system that we use to compute approximations of the moments of x(t).

We fix N T ∈ N and define

x i (t) ∈ R n i , i ∈ 1, . . . , N T , by 1 x 1 (t) • • • x NT (t) = E(N T , N T ) t 1 E[x 0 ] • • • E[x [NT] 0 ] . (7)
Here, the vector x i (t) represents an approximation of the moment E[x [i] (t)] that is computed using only our knowledge of the first N T moments of x 0 .

By letting y(t

) [ 1 x 1 (t) x 2 (t) • • • x NT (t) ] ,
we obtain the so-called "truncated system", which is a discrete-time linear time-invariant system given by

y(t + 1) = E(N T , N T ) y(t), t ∈ N, y(0) = [ 1 E[x 0 (t)] • • • E[x [NT] 0 (t)] ] . (8)
The truncated system allows us to iteratively compute approximations of the moments of x(t) at consecutive time instants. Moreover, the approach only requires an offline computation of the matrix E(N T , N T ).

Approximation Error Bounds

We now investigate the approximation error introduced by truncation. Let e i (t) ∈ R n i denote the approximation error of the ith moment, that is,

e i (t) E[x [i] (t)] -x i (t).
In what follows, we provide upper bounds to e i (t) (the error on the ith moment introduced by the approximation) for any norm • induced by an inner product. These bounds allow us to use various techniques to study the distribution of the state at future time steps. We illustrate this in Section 3.3 by using tail probability approximations [START_REF] Gray | A general method for approximating tail probabilities[END_REF] to compute an upper bound of the probability to be outside of a given safe region.

Let j 0 ∈ {0, . . . , N T }. Our goal is to obtain an upper bound to e j0 (t) . First, by (6), we have

E[x [j0] (t)] = j0dS j1=0 E j0,j1 E[x [j1] (t -1)] = j0dS j1=0 E j0,j1 j1dS j2=0 E j1,j2 E[x [j2] (t -2)] = j0dS j1=0 E j0,j1 j1dS j2=0 E j1,j2 • • • jt-1dS jt=0 E jt-1,jt E[x [jt] 0 ],
and similarly, by (7),

x j0 (t) = NT j1=0 E j0,j1 NT j2=0 E j1,j2 • • • NT jt=0 E jt-1,jt E[x [jt] 0 ].
By observing that E j,k = 0 if k > jd S , we obtain from the two equations above,

e j0 (t) = E[x [j0] (t)] -x j0 (t) = j0dS j1=NT+1 E j0,j1 j1dS j2=0 E j1,j2 • • • jt-1dS jt=0 E jt-1,jt E[x [jt] 0 ] + NT j1=0 E j0,j1 j1dS j2=NT+1 E j1,j2 • • • jt-1dS jt=0 E jt-1,jt E[x [jt] 0 ] + • • • + NT j1=0 E j0,j1 • • • NT jt-1=0 E jt-2,jt-1 jt-1dS jt=NT+1 E jt-1,jt E[x [jt] 0 ]. ( 9
)
As an immediate consequence, we get the following: Proposition 3. Consider the truncated approximation of the moments of system (1) with truncation limit N T ∈ N.

If j 0 d t S ≤ N T , then x j0 (t) = E[x [j0] (t)].
Proof. We show that, if j 0 d t S ≤ N T , then e j0 (t) = 0. To this end, it is enough to show that, for all lines of (9) and sequences (j 1 , . . . , j t ) of relevant indices, at least one E j k-1 ,j k is equal to 0. Let us pick the ith line and any sequence as above. In particular, j i-1 d S > N T , hence j i-1 > j 0 d t-1 S ≥ j 0 d i-1 S , so there exists k ∈ {1, . . . , i -1} such that j k > j k-1 d S , which proves the desired result. 2 This shows that, for large values of truncation limit N T , the proposed method computes exact moments E[x [j0] (t)] for j 0 and t small enough. Note that this is due to the discrete-time nature of the finite-dimensional polynomial system (1). In the continuous-time case, approximation errors cannot be avoided in general [START_REF] Forets | Explicit error bounds for Carleman linearization[END_REF].

However, the exact value of e j0 (t) is generally hard to compute. Indeed, since E[x [i] (t)] = x i (t) + e i (t), if e i (t) could efficiently be computed, then so would E[x [i] (t)], and there would be no need in using the truncated system. In the rest of this section, we thus come up with several upper bounds for e j0 (t) that can be efficiently computed.

First, observe that ξ = max

0≤j≤j0d t S E[x [j]
0 ] can be efficiently computed in some cases. An obvious situation is when the position x 0 is determined, in which case we have ξ = max{1, x 0 j0d t S }. Another case is when x 0 obeys a wellknown distribution whose moments are easy to compute, such as uniform or normal distributions. Nevertheless, another case is when the system satisfies x(t) ∈ R and 0 ≤ x 0 ≤ 1, in which case E[x

e j0 (t) = j0d t S j=0 E j E[x [j] 0 ], (10)
where E j ∈ R n j 0 ×n j is constructed by additions and multiplications of E i,j , and thus can be computed offline (note that E j is dependent on t, but we keep this implicit for readability). From this, we can derive a global bound

e j0 (t) ≤ ξ j0d t S j=0 E j ,
where j0d t S j=0 E j can be computed offline. We can further refine (10) to consider a single line i ≤ n j0 of this equation, for which we get (e j0 (t)

) i = j0d t S j=0 v j,i E[x [j]
0 ], where v j,i ∈ R 1×n j is the ith row of E j . By repeated application of triangle and Cauchy-Schwarz inequalities, we obtain

|(e j0 (t)) i | ≤ ξ j0d t S j=0 v j,i
, where j0d t S j=0 v j,i can also be computed offline. This bound can, however, be crude in practice as the norm gets distributed over all sums and products. Here we show how to compute tighter bounds while maintaining a reasonable computational cost. For any subset J ⊆ {0, . . . , j 0 d t S }, we have the following:

|(e j0 (t)) i | ≤ j∈J v j,i E[x [j] 0 ] + ξ J j ∈J v j,i , (11)
where ξ J = max

j ∈J E[x [j] 0 ] .
The idea is that J is a set of indices where one should avoid distributing the norm over the sum. One should pick J to consist of those indices where the distribution is too crude and makes the error bound loose. In order for this method to be computationally efficient, one should J that is of relatively small size, e.g., |J| = O(t). One possible way to choose J is to fix a size k and return the set of k indices such that v j,i is largest for those indices; another way is to return the set of k indices such that E[x

[j] 0 ] is largest. Example 4. Suppose x 0 is drawn from a truncated normal distribution over the interval [0, 1]. Then ξ = 1 and ξ J = E[x [jmin] 0
] , where j min is the smallest number that is not in J. In particular, if J is chosen as the set of indices such that E[x

[jmin] 0

] is largest among j = 0, . . . , j 0 d t S , then we have J = {0, . . . , |J| -1} and

ξ J = E[x [|J|] 0 ] .

Tail Probability Approximation

In this section, we use the error bound introduced by the truncated system, which is derived in the previous section, to give a lower bound on the probability of staying inside a given safe region.

Suppose we know some approximations ( x 2 (t)) ni+j (henceforth denoted ( x 2 (t)) i,j ) of expectations E[x i (t)x j (t)] for all i, j ≤ n (where n is the dimension of the state space, as in (1)), with respective error bounds ε i,j , as well as approximations ( x 1 (t)) i of E[x i (t)] with error bounds ε i . Finally, assume that we know a global error bound x 1 (t) -E[x(t)] ≤ ε. These different bounds can be found by applying the techniques in the previous section.

We can then bound the probability to be outside of a "safe" region centered on x 1 (t) by the following proposition. Proposition 5. Let • be the Euclidean norm. For α > ε,

P( x(t) -x 1 (t) ≥ α) ≤ n i=1 ( x 2 (t)) i,i + ε i,i -max 0, (|( x 1 (t)) i | -ε i ) 2 (α -ε) 2 .
Proof. Since α > ε, by Chebyshev's inequality:

P( x(t) -x 1 (t) ≥ α) ≤ P( x(t) -E[x(t)] ≥ α -ε) ≤ E[ x(t) -E[x(t)] 2 ] (α -ε) 2 = n i=1 E[x 2 i (t)] -E[x i (t)] 2 (α -ε) 2 ,
and we conclude using

E[x 2 i (t)] ≤ ( x 2 (t)) i,i + ε i,i and E[x i (t)] 2 = |E[x i (t)]| 2 ≥ max 0, (|( x 1 (t)) i | -ε i ) 2 .
The safety region discussed in Proposition 5 corresponds to a ball in R n with radius α. We note that safety regions can be generalized to ellipsoids by changing the norm

• to a weighted norm • P induced by a positivedefinite matrix P ∈ R n . Moreover, semi-norms induced by positive semi-definite matrices can be used for problems where the safety of only certain entries of the state vector is explored. A similar argument may be used to prove a general formula for any semi-norm induced by a positive semi-definite matrix, but we omit it as it will not be used in the examples presented in this paper.

NUMERICAL EXAMPLES

In this section, we provide two numerical examples and experimental results to illustrate our techniques.

Stochastic Logistic Map

Consider the stochastic logistic map as studied by [START_REF] Athreya | Random logistic maps[END_REF], which is given by

x(t + 1) = r(t)x(t)(1 -x(t)), t ∈ N, x(0) = x 0 ,
where x 0 , r(0), r(1), . . . are mutually independent random variables; x 0 takes values in [0, 1] and r(t) all take values in [0, 4]. The scalar x(t) ∈ [0, 1] represents the population of a species subject to growth rate r(t). This system can equivalently be represented by (1) with d S = 2, F 0 (t) = 0, F 1 (t) = r(t), and F 2 (t) = -r(t).

For experiments, we chose all r(t) to be uniformly distributed over the interval [0.3, 0.7], and x 0 to follow a normal distribution of mean 0.5 and standard deviation 0.1 truncated to [0, 1].

Moment Approximation via Truncated System. We first compare our moment approximations for different truncation limits to the true value of the moments (computed using our method with N T = 256). In Figure 2, we plot the first and second moments of the truncated system with different truncation limits N T . In most cases, we observe that the system with the higher truncation limit gives a better approximation. There is one exception in the second moment approximation at t = 6; the system with the lower limit gives a better approximation. This behavior seems to be similar to the one in [START_REF] Forets | Explicit error bounds for Carleman linearization[END_REF], where a truncated system with a higher N T diverges to infinity more quickly, while it gives a better approximation over a longer time interval. (In this case, this is due to the fact that the actual system is converging to 0, and a system with truncation limit N T diverges to infinity a polynomial of degree N T .) We also observe that larger N T is required to obtain good approximations of higher moments. This is a natural consequence, as truncation discards more information on the dynamics of higher moments.

Error Bound on Moment Approximations.

Next, we evaluate our approximation method of error bounds for moments. Figure 3 shows error bounds given by the inequality (11) with different sizes of J, with parameters N T = 16, t = 4, and j 0 = 2. The index set J contains the indices j where E[x

[j]
0 ] is the largest, i.e., J = {0, . . . , |J| -1}. We observe that the error bound quickly decreases as |J| increases. This supports our expectation that we can use our error bounds for tail probability analysis with larger parameters (t, N T , and d S ).

Tail Probability Analysis. Lastly, we provide a result on tail probability analysis via the method in Section 3.3. We computed the error bound for 0 ≤ t ≤ 5 using ( 11 We observe that the size of safety intervals given by our tail probability analysis is reasonably small. It becomes cruder in later time steps, especially at t = 5. This is expected, as the approximation error of moments, which is a bottleneck in refining the error bounds, becomes larger as time progresses (cf. approximate 2nd moment in Fig. 2).

There are two major advantages of our method compared to Monte Carlo simulation. One is that our technique computes moment approximations much faster than Monte Carlo (even for large N T and a small number of samples) because we do not rely on generating random numbers. This advantage is highlighted in Table 1, which contains the online computation times for Monte Carlo simulations and our approach, averaged over 100 runs. Another advantage is that our safety interval gives a theoretical guarantee on probabilistic safety that cannot be achieved by Monte Carlo simulations.

Vehicle Dynamics

Our second example is an application to autonomous driving. For safety guarantees, autonomous vehicles need to predict their future positions. One way to achieve this is set-based reachability, as advocated by [START_REF] Althoff | Online verification of automated road vehicles using reachability analysis[END_REF]. To use their method, they must consider systems with bounded disturbances and use linearization around an equilibrium, that is, they approximate a polynomial system by a linear one, using Lagrange remainders. Carleman linearization allows taking the effect of higher dimensions of the system into account more precisely than Lagrange remainders. Moreover, our approach is probabilistic, while theirs is set-based, so the two approaches give different types of guarantees.

We consider a scenario in which, at each time step, the autonomous vehicle measures its current position with some known sensor error distributions, computes moments of its current state based on this measurement and the known error distribution, and predicts its future positions (and moments) up to t steps ahead in time by applying the truncated system to these moments. This process is repeated forever, leading to the framework in Figure 5.

More precisely, we consider the kinematic bicycle model of a vehicle [START_REF] Kong | Kinematic and dynamic vehicle models for autonomous driving control design[END_REF], which we rewrite as the following equivalent polynomial system ṗx (t) = v(t)c(t), ṗy (t) = v(t)s(t), The second-order Taylor expansion of the model above gives the following discrete-time approximation:

ψ(t) = v(t) sin β, v(t) = a(t), ċ(t) = - s(t)v(t) sin β , ṡ(t) = c(t)v(t) sin β ,
p x (t + ∆) = p x (t) + ∆c(t)v(t) + ∆ 2 2 a(t)c(t) - s(t)v 2 (t) sin β , p y (t + ∆) = p y (t) + ∆s(t)v(t) + ∆ 2 2 a(t)s(t) + c(t)v 2 (t) sin β , ψ(t + ∆) = ψ(t) + ∆ v(t) sin β + ∆ 2 2 a(t) sin β, v(t + ∆) = v(t) + ∆a(t), c(t + ∆) = c(t) -∆ s(t)v(t) sin β - ∆ 2 2 c(t)v 2 (t) sin 2 β 2 + a(t)s(t) sin β , s(t + ∆) = s(t) + ∆ c(t)v(t) sin β + ∆ 2 2 - s(t)v 2 (t) sin 2 β 2 + a(t)c(t) sin β ,
where ∆ > 0. To describe the evolution of the states of the vehicle at times 0, ∆, 2∆, . . ., we write this system in the form of (1). In particular, consider the discrete-time instant t ∈ N corresponding to the continuous time t∆. By we obtain (1) where d S = 3 and the coefficients F 0 (t), . . . , F 3 (t) depend on ∆, β, , and a(t). We consider the setting where the acceleration values a(0), a(1), . . . are independent uniform random variables over [0.9, 1], ∆ = 0.1, β = π/8, = 2.5, and, for the initial state, p x (0), p y (0), v(0), ψ(0) are independent Gaussian random variables with mean 0 and standard deviation 0.1, and c(0) = cos(ψ(0) + β) and s(0) = sin(ψ(0) + β).

Experimental Results. Figure 6 shows the expected trajectory of the vehicle as approximated by our method for different truncation limits, as well as the empirical distribution computed by 10000 runs of Monte Carlo simulation and the mean of that distribution. Figure 7 shows the distance d( x 1 (t), MC(t)) between x 1 (t) and the mean of the empirical distribution for the same truncation limits.

Figures 6 and7 show that larger truncation limits give systems that follow the empirical distribution better. It also shows that, for a fixed truncation limit, the distance to the empirical system grows larger with time. The only exception is around t = 6 for N T = 2, which can be explained by the fact that the trajectory of the truncated system crosses that of the empirical one. Notice that, N T = 8 gives an exact computation at t = 2, which means that the error at this point (corresponding to the horizontal line in Figure 7) is entirely induced by Monte Carlo. Therefore, it is pointless to consider errors that are of the same order as the one introduced by Monte Carlo at t = 2.

Also note that the largest truncation limit we use in this experiment is N T = 8, which makes the online computation very fast.

CONCLUSION

We have explored the problem of computing moments of discrete-time polynomial stochastic systems. Through a truncated Carleman linearization approach, we have obtained an iterative procedure for approximating the moments of the state of the system across different time instants. Furthermore, we have provided tail-probability bounds to assess the safety properties of the system.

In our analysis, we have used spherical safety regions characterized using the Euclidean norm. By changing the norm to weighted norms induced by positive-definite matrices, ellipsoidal safety regions can also be considered. We note that safety analysis for a non-ellipsoidal region can be conducted by using ellipsoids that are inside that region. Finding the largest such ellipsoid allows us to do safety analysis in the least conservative way. This is possible by solving a matrix optimization problem, which we leave for future work.

One of our future research goals is to use the approximate moments to reconstruct an approximation of the distribution of the state of the system. There exist different techniques to achieve this, as given for example by [START_REF] Schmüdgen | The Moment Problem[END_REF] or [START_REF] John | Techniques for the reconstruction of a distribution from a finite number of its moments[END_REF].

Another direction for future work is to make the algorithm more efficient by using the symmetry of Kronecker powers to reduce the sizes of the matrices involved in the computations, which would allow us to compute even higher moments and more time steps in the future.

  decreasing and we have ξ = E[x [0] 0 ] = 1. Using ξ, we can derive bounds for e j0 (t) by first rewriting (9) as

  ) with N T = 16 and |J| = 6t, where J is taken in the same way as the previous section. The results of the analysis are depicted in Figure 4. Orange intervals indicate the 95%probability neighborhoods of x 1 (t) computed by using Proposition 5. Blue intervals with a solid line indicate the region where 95% of 10000 Monte Carlo simulations closest to its mean belong. Dotted intervals indicate the range of 10000 Monte Carlo simulations.

Fig. 3 .

 3 Fig. 3. Error bound on moment approximations.

  for t ≥ 0, where p x (t) ∈ R and p y (t) ∈ R represent the X-Y coordinates of the mass-center of the vehicle, v(t) ∈ R denotes its speed, ψ(t) its inertial heading, and a(t) ∈ R its acceleration. The constants β ∈ R and > 0 respectively denote the angle of velocity and the distance from the vehicle's rear axle to its mass-center. The scalars c(t) and s(t) are auxiliary variables that are introduced to obtain the polynomial model above from the original model of[START_REF] Kong | Kinematic and dynamic vehicle models for autonomous driving control design[END_REF] (which involves trigonometric terms), using the same techniques as[START_REF] Carothers | Some properties of solutions to polynomial systems of differential equations[END_REF].

Fig. 6 .

 6 Fig. 5. Moment propagation framework for autonomous vehicles.

Fig. 7 .

 7 Fig. 7. Distance to the mean of the empirical distribution.letting

Table 1 .

 1 Comparison of online computation times.

	Method	Monte Carlo	Moment propagation
	Parameters	num. samples 10 10 4	N T 4 16 64	256
	Time (µs) 2.9e10 3 3.4e10 6 11 14 30	93