C Strub

X M Meyer
email: xuan.meyer@ensiacet.fr

C Brandam

A Lebrihi

Stoichio-kinetic scheme extraction from experimental data: Application to antibiotics production with Saccharothrix algeriensis

Almost a quarter of the deaths in the world results of infectious diseases. They are more and more caused by bacteria having developed a (or several) resistance to antibiotics [START_REF] Mukhopadhyay | The RNA polymerase "Switch Region" is a target for inhibitors[END_REF]. This phenomenon is worrying for humanity and stresses the importance of constantly renew the pool of bio-active molecules. Thus, the production of new bio-active molecules on the pathogenic currently available antibiotic resistant strains is the subject of more and more interdisciplinary research projects. In this context, a new bacterial specie Saccharothrix algeriensis NRRL B-24137 (Zitouni et al., 2004) has been isolated in 1992 from the Sahara desert (Algeria). This filamentous bacterium belongs to the actinomycetes family. It produces molecules belonging to the dithiolopyrrolones family which present antibacterial, antifungal, and interesting anticancer activities (Webster et al., 2000 ;[START_REF] Oliva | Antimicrobial properties and mode of action of the pyrrothine holomycin[END_REF][START_REF] Minamiguchi | Thiolutin, an inhibitor of huvec adhesion to vitronectin, reduces paxillin in huvecs and suppresses tumor cell-induced angiogensis[END_REF]. A previous work (Bouras, 2005;Bouras et al., 2006a) highlighted the influence of the culture medium composition on the dithiolopyrrolones production, including on the production of new molecules never described in the scientific literature [START_REF] Bouras | Precursor-directed biosynthesis of new dithiolopyrrolone analogs by Saccharothrix algeriensis NRRL B-24137[END_REF]. This capacity to produce antibiotics upon request by precursor-directed biosynthesis facts that the study of Saccharothrix algeriensis presents a fundamental undeniable interest, in addition to an obvious useful interest. As this bacterium has been recently discovered and characterised, no deep knowledge on its metabolism is available. Several studies have been undertaken in parallel to better understand the dithiolopyrrolones production pathway by Saccharothrix algeriensis. For the study presented here, a macroscopic approach has been chosen to investigate the bacterium metabolism. No data coming from "omics" are used but the "reaction engineering" approach adopted to analyse the metabolism from a qualitative and quantitative point of view enables to get more knowledge about the mechanism of dithiolopyrrolones production. Two broad approaches exist to model microbial bio-productions:

-The black box models which do not require any knowledge on the microorganism. Only the inputs (substrates) and the outputs (products) of the system are taken into account. These models do not enable to draw assumptions on the metabolic pathways. It remains purely descriptive.

-The structured models, like the "Metabolic Flow Analysis" model, which favour a more systemic approach of the microorganism. They are generally complex and require a good knowledge of the metabolism network of the bacteria. The quantification of hundreds of flux from a limited measurements number compels to assume that all intermediate compounds are in a quasi-stationary state. This drawback may be avoided if intracellular compounds can be measured. In this work, we adopted an intermediate approach based on an analysis of experimental data. It allows establishing a global reaction scheme of the bacterium metabolism. The advantage of this type of model is that, while remaining at the extra cellular scale, with experimental information relatively easy to obtain, it allows getting information on the metabolism of the cell. Stoechio-kinetic models consist of the decomposition of bacterial growth in some reactions of invariant stoechiometry with time. This "reaction engineering" approach is widely used in chemical engineering to determine the occurring reactions in a complex reaction scheme. It is based on the bilinear structure of the component mass balances in batch reactor. The first adaptation of the factorial analysis to the treatment batch reactor data is due to [START_REF] Hamer | Stoechiometric interpretation of multireaction data: application to fed-batch fermentation data[END_REF]. This method underwent several adaptations [START_REF] Bonvin | Target factor analysis for the identification of steochiometric models[END_REF]Harmon et al., 1995 ;[START_REF] Fillon | Data reconciliation based on elemental balances applied to batch experiments[END_REF][START_REF] Moreau | Contribution au contrôle de la culture de Brevibacterium linens en réacteur batch[END_REF][START_REF] Castro-Martinez | Brettanomyces bruxellensis : Etude Métabolique, Cinétique et Modélisation. Influence des facteurs environnementaux[END_REF]. This paper presents the results of this approach applied to the growth and production of thiolutine (the main dithiolopyrrolone produced) by Saccharothrix algeriensis, on a semi-synthetic medium. The consistency of the experimental data is first checked in order to further perform the qualitative and quantitative analysis on reliable data. These three parts are presented successively hereafter.

MATERIALS AND METHODS

Microoganism and cultures conditions

Sa. algeriensis NRRL B-24137 was used in this study. Microbial spores were obtained from solid culture on petri dishes filled with conservation medium. They were maintained in 25% glycerol at -20°C. A 100 ml volume of semi-synthetic medium was inoculated by 3.5 ml of this suspension and incubated on a rotary shaker (New Brunswick Scientific, Edison, NJ, USA) at 250 rotation per minute (rpm) at 30°C for 52 h. Five ml of the resulting preculture was used to inoculate each culture. Before inoculation, the pH was adjusted to 7 by addition of 1 mol.L -1 NaOH. Cultures were performed in NBS reactors containing 2 L of medium. Culture lasted for at least one week. pH was maintained at 7±0.035 by the automatic addition of 1 mol.L -1 NaOH and 1 mol.L -1 HCl solutions. Aeration rate of one vessel volume per minute (0.5 vvm) was employed. The agitation rate was controlled to keep the dissolved oxygen level above 30 % saturation with a starting rate of 150 rpm. The pH and dissolved oxygen level were monitoring using Ingold specific electrodes. Temperature was regulated at 30°C. Percentages of O 2 , CO 2 exhaust gas were determined by a gas analyzer (Servomex 4100, paramagnetic transductor for O 2 and infrared transductor for CO 2 ).

Media composition

The conservation media had the following composition (per litter of distilled water): 10 g malt extract, yeast extract 4 g, glucose 4 g and agarose 18 g. Semi-synthetic medium used as growth and production medium contained (per litter of distilled water): glucose 15 g, yeast extract 2 g, NaCl 2 g, (NH 4 ) 2 SO 4 2 g, KH 2 PO 4 0.5 g, K 2 HPO 4 1g, MgSO 4 .7H 2 O 0.2g, CaCl 2 .2H 2 O 1 g and MOPS 2 g, Uracil 20 mmol. Attempts were made to substitute yeast extract to better control the medium composition without success [START_REF] Strub | Investigations on Saccharothrix algeriensis grow synthetic media[END_REF]. The concentrations in chemical elements and free amino acids of the yeast extract used during this study (288620 Bacto TM Yeast Extract, Technical) are presented in tables 1. Its global elementary composition was evaluated to C 3.32 H 6.54 O 2.17 N 0.79 .

Analytical procedures

For estimation of dry cell weight (DCW), 3 ml samples of homogenised culture broth were centrifuged at 16000 g for 10 min in preweighed Eppendorf tubes. Pellet was washed first with 0.25 mol.L -1 NaOH solution, then with 0.35 mol.L -1 HCl solution and then with distilled water. Eppendorf tubes containing pellet were dried at 105°C for 48 h, and weight after cooled for 30 min in a dessicator (Bouras et al., 2007). The measurement relative error is 5 %. Supernatant coming from DCW determination was filtered at 0.2 µm and used to measure metabolites. The analysis of dithiolopyrrolone antibiotics was carried out by non polar chromatography (HPLC, Bio-Tek Instruments) with the method described by Bouras et al. (2006b). Glucose was quantified with a biochemical analyzer with fixed glucose oxydase enzymes on membrane (YSI2700 select). Glucose was determined with amperometric quantification after its enzymatic oxidation. The amperometer answer is linear for a concentration in glucose between 0 and 25 g.L -1 . The measurement relative error is 3 %. Ammonium ions and α-amino nitrogen were quantified using specific enzymatically methods (Diagnostics Ammonia kit from Boehringer-Mannheim, using glutamate deshydrogenase and Microdon Kit using glutamate oxydase) and an automatic multiparametric analyser (Mascott Lisabio). The signal answer is linear for concentrations in α-amino nitrogen and ammonium ion between 0 and 500 mg.L -1 . The measurement relative error is 5 %. Amino acids quantification was done with AminoQuant HP1900 (non polar column C 18 ) and associated protocol. Before analysis, proteins contained in the supernatant were precipated with methanol 75 % v/v one night at 4°C. Amino acids were automatically derivated with orthophthalaldehyde-9-fluorenylmethyloroformate (OPA-FMOC).

Data reconciliation

Experimental values were reconciliated according to the chemical elemental (C, N, O, H, …) mass balances. The objective is double: first, it enables to check the consistency of experimental data and second, it can slightly correct them so that the mass balances on chemical elements are completely satisfied. Data reconciliation consists in solving the optimization problem [START_REF] Fillon | Data reconciliation based on elemental balances applied to batch experiments[END_REF]:

         0 X M ] ) X X ( V X m X 2 1 [ Min m 1 T X ˆ) (
Where X m is the (nc, nobs) matrix of the measured values of the concentrations of the n c components on the n obs samples, is the (nc, nobs) matrix of reconciliated values and V the measurement X variance-covariance (nc,nc) matrix calculated from the measurement relative errors.

The M =0 constraint expresses the chemical element balance between two analyzed samples:

X ˆ               E 0 E 0 E 0 0 E M     0 0 E  with E the elemental decomposition matrix of the different measured compounds.
To check the data consistency, three statistical tests are performed: the Generalized Likelihood Ratio (Narasimhan and Mah, 1987), the Iterative Measurement Test [START_REF] Serth | Gross error detection and data reconciliation in steam metering system[END_REF] and the Residual Criterion [START_REF] Wang | Application of macroscopic balances to the identification of gross measurement errors Biotech[END_REF].

Stoechiometric coefficients and reaction extends determination

It is possible to extract a reaction scheme and corresponding reactions extends from the molar hold-up by exploiting the bilinear structure of the component mass balances. In the case of a batch reactor in which n r reactions are performed and n c components are measured on n obs samples, the data matrix D, consisting of n obs rows containing the cumulative molar amount of the n c measured species can be written as the product of two matrixes R and S the matrix of the reaction extend (n obs ,n r ) and the matrix of stoechiometric coefficients of the n r reactions (n r ,n c ) respectively. When a qualitative analysis of the raw data enables to assume the active reactions, a direct determination of R and S by solving the constrained optimization problem based on reconciliated data can be done [START_REF] Moreau | Contribution au contrôle de la culture de Brevibacterium linens en réacteur batch[END_REF]:

            RS D balances elemental the verifies S n 1, = j and n 1, = i for 0 s s et 0 s ) D Max( ) D D ( min calc r obs j i, j 1, + i j i, i exp j , i 2 j calc j , i exp j , i S , R

RESULTS

Data reconciliation

For the amino-acids to have a sufficient weight in the elemental balance to perform data reconciliation, the 20 amino-acids from the yeast extract have been assimilated to one "pseudo amino-acid" with the composition C 4.73 H 9.69 O 2.38 N 1.26 calculated from the weight of the different amino-acids in the yeast extract. This can be done as the main amino acids are consumed nearly at a same rate during the first phase of the experiment (fig. 1). The biomass composition has been evaluated at different time during growth and the results did not show any significant evolution. So, biomass was assimilated to a chemical compound with the formula C 3.82 H 6.75 O 1.89 N 0.8 assumed to be constant during the whole experiments. Data reconciliation enlightens a lack of carbon at the very beginning of the fermentation. This indicates that a non measured carbohydrate coming from the yeast exact is consumed by Sa. algeriensis preferably to glucose as carbon source. This compound was assimilated to glucose for the quantitative analysis with the quantity given by the data reconciliation (fig 2). No deviation was found on nitrogen balance.

Qualitative analysis of the growth of Sa. algeriensis in liquid medium

The continuous measurement of the carbon dioxide production by Sa. algeriensis shows that its growth on the semi-synthetic medium can be divided into four phases (fig. 3). The corresponding evolution of the biomass, the glucose, the -amino nitrogen, the ammonium and the thiolutine are represented on figure 4. Fig. 3: carbon dioxide production rate (g.L -1 h -1 ) versus time (h) during Sa. algeriensis growth on a semi-synthetic medium (YE 2g.L -1 , Uracil 20 mM) at 30°C, pH 7, 0.5 vvm and 150 rpm.

For all experiments, the four phases have been observed. First, the amino-acids are consumed as a carbon and nitrogen source. When this substrate is exhausted, a secondary metabolism of the bacterium is activated to produce thiolutine. No growth is observed during this second phase but some carbohydrates are consumed certainly as energy source for maintenance. Then the primary metabolism is activated again to achieve the growth on glucose and ammonium as carbon and nitrogen source respectively. During this phase, a thiolutine degradation is observed. The fourth phase begins as the ammonium is exhausted and again a thiolutine production is observed.

The qualitative analysis of the experimental results indicates a sequential consumption of the carbon source substrates, a probable diauxic growth, a preferential consumption of small molar weight aminoacids and a thioluthine production partially decoupled from the growth and initiated when a substrate deprive occurs.

Quantitative analysis of the growth of Sa. algeriensis in liquid medium

From the qualitative analysis, we suggest to represent Sa. algeriensis growth and thiolutine production through a set of four reactions: growth on amino-acids (R1), maintenance on glucose (R2), growth on glucose (R3) and thiolutine production (R4). These reactions occur during the different growth phases formerly identified as follows:

 Phase 1: growth on amino-acids and maintenance on glucose  Phase 2: maintenance on glucose and thiolutine production  Phase 3: growth on glucose  Phase 4: maintenance on glucose and thiolutine production Two more reactions could have been written to express the thiolutine degradation and the biomass decrease but as the compounds resulting from these degradations were not determined nor measured, these reactions were omitted. To identify the reaction scheme stoichiometric, thiolutine data were corrected to consider no degradation and biomass was considered to remain constant during the last phase.

For each reaction a stoechiometric relationship was established: These coefficients together with the reaction extends have been evaluated for experiments at different glucose and yeast extract initial concentrations (Table 2) with a good agreement between calculated and experimental molar content values, as shown on Figure 5 for biomass and thiolutine. Stoichiometric coefficients do not appear to be sensitive to the glucose concentration but different values where obtained for different yeast extract initial concentrations, indicating a probable different metabolism pathway. This aspect was not analysed further in this study. Finally, for a yeast extract initial concentration of 2 gL -1 , the following reaction scheme can represent the growth of Sa. algeriensis and the thiolutine production. The growth on amino acids (R1) is carried out with a weak respiratory activity (Q R =0.689) compared to those obtained for the growth on glucose (R3), Q =1. The assimilated amino acids would be R integrated directly into the level of the proteinic synthesis without participating to the central metabolism of the micro-organism. The reaction extends presented on figure 6 confirms the assumptions made on the occurence of the different reactions: the growth is carried out initially on the free amino acids and the unknown substrate brought by the yeast extract assimilated to glucose (R1). When free amino acids are exhausted, the growth takes on glucose (R3). The production of thiolutine (R4) is uncoupled from the growth. It is carried out during the phase of maintenance on glucose (R2). A kinetic model for each reaction will be established from these reaction extends to express the reaction rate as a function of the substrate concentrations.

R1

CONCLUSIONS

This study shows how it has been possible to acquire some knowledge on Sa. algeriensis metabolism from extra cellular compound measurements. The "reaction engineering" approach enables to perform a qualitative and a quantitative analysis on reliable data with respect to elemental mass balances. So, a four reaction scheme has been suggested as a macroscopic view of the growth of Sa. algeriensis and the associated production of thiolutine. It was shown that Sa. algeriensis was able to grow on both free amino acids and glucose as carbon source. The amino acids of weak molar weight are consumed preferentially to glucose. They seem to be directly integrated to the proteinic synthesis and enriched the medium with ammonia. This ammonia is then assimilated during the growth on glucose. It has been established that the thiolutine production was the result of a secondary e behaviour of this bacterium of a re kind. They give tracks for an exploration of the metabolism to an intracellular level and pose the ases for the production of dithiolopyrrolone at an industrial scale. itouni, A., Lamari, L., Boudjella, H., Badji, B., Sabaou ,N., Gaouar, A., Mathieu, F., Lebrihi, A. & Labeda, D.P. (2004c) Saccharothrix algeriensis sp. nov., isolated from a saharan soil. Int. J. Syst. Evol. Microbiol., 54, 1377-1381.

Fig. 1 :

 1 Fig.1: amino-acid conversion during Sa. algeriensis growth on a semi-synthetic medium (glucose 15 gL -1 , YE 2g.L -1 , Uracil 20 mM) at 30°C, pH 7, 0.5 vvm and 250 rpm.

Fig. 2 :

 2 Fig.2: glucose (raw data (), reconciliated data ()) and amino-acids (raw data (), reconciliated data ()) during Sa. algeriensis growth on a semi-synthetic medium (YE 2g.L -1 , Uracil 20 mM) at 30°C, pH 7, 0.5 vvm and 250 rpm.

Fig. 4 :

 4 Fig. 4 : Batch culture of Sa. algeriensis on a semi-synthetic medium (YE 2gL -1 , Uracil 20 mM) at 30°C, pH 7, 0.5 vvm and 150 rpm. a) Biomass (), α-amino N () et glucose () b) Thiolutine (),αamino N () et ammonium ().

Fig. 5 :

 5 Fig.5: calculated values versus experimental values for biomass (a) and thiolutine (b) for a batch culture of Sa. algeriensis on semi synthetic media (5 g.L -1 glucose and 2 g.L -1 YE) pH 7, 30°C, 250 rpm.

  reaction extends for a batch culture of Sa. algeriensis on a semi-synthetic medium (5 g.L -1 de ), pH 7, 30°C, 250 rpm. F glucose, 2 g.L -1 YE

  metabolism, decoupled from growth. In our experimental conditions it seems to be initiated as soon as a nitrogen source depletion occurs. These results open an interesting way to deepen the knowledge on th ra . (2000) Anticancer property of dithiolopyrrolones. U.S. Patent 6,020,360.

Table 1 :

 1 Concentrations in free amino-acids and other elements of the yeast extract (288620 BactoTM Yeast Extract, Technical)

	Amino-acids	[amino acids] (µmol) per g of YE	Atomic element	[atomic element] (µg) per g of YE
	ASP	142	P	9751
	THR	144	S	36877
	SER	200	W	222
	GLU	119	Zn	286
	GLN	4	Co	5
	PRO	42	Pb	4
	CYS	30	B	41
	GLY	204	Mn	3
	ALA	415	Fe	117
	VAL	236	Mg	802
	MET	51	Al	276
	ISO	164	V	66
	LEU	311	Be	2
	TYR	32	Ca	935
	PHE	140	Zr	1
	LYS	156	Ba	1
	HIS	47	Na	1526
	TRP	15	K	73542
	ARG	109
	ASN	108

Table 2 :

 2 variability of the stoechiometric coefficients with the glucose and the yeast extract initial concentrations

	[glucose] 0 (gL -1 )	[Yeast Extract] 0 (gL -1 )	Other compound	x 1	x 3	t 4
	15	2	Uracil 20 mM	1.2 0.56 0.16
	15	2	benzoic acid 5mM Uracil 20 mM	1.2 0.61 0.16
	15	2	Uracil 20 mM	1.2 0.54 0.16
	15	2	-	1.2 0.52 0.16
	3	2	-	1.2 0.44 0.16
	5	2	-	1.2 0.44 0.16
	8	2	-	1.2	0.5 0.16
	15	6	Uracil 20 mM	0.94 0.32 0.16
	15	4	Uracil 20 mM	0.74 0.26 0.16
			Uracil 20 mM
			Amino-acids mix at
	15	2	the same	0.95 0.53 0.16
			composition as in YE at 2gL -1
	15	2	Uracil 40 mM	1.2 0.54 0.16


M Teles

Dos Santos

G A C Le Roux
email: galoroux@usp.br

V Gerbaud

VIRTUAL LABORATORY FOR THE DESIGN OF STRUCTURED LIPIDS FROM VEGETABLE OILS RESOURCES

Keywords: Lipids, triacylglycerol, vegetable oils, CAMD, CAMbD, group contribution

The use of available vegetable oils natural resources for producing improved molecules/mixtures is a promising area in tailor-made lipids. For this, it is necessary a multiscale modelling of different properties, taking into account the relations between molecular structure and desired properties. This modelling by predictive methods can be coupled with optimization techniques for screening among potential candidates the optimal ones according to the set of end-use properties. The present work shows some aspects of the development of a framework for lipids design using groups-based property prediction models and presents some thermodynamic properties and viscosity of triacylglycerols using Group Contribution Models.

INTRODUCTION

The search for renewable resources in the chemical industry plays is essential for achieving sustainable development. Vegetable oils are renewable lipid-based raw materials. They are already used in food, biofuels, cosmetic, soaps, pharmaceutical and lubricants industry. But new products and applications are sought with tailor-made lipids [START_REF] Osborn | Structured Lipids-Novel Fats with Medical, Nutraceutical, and Food Applications[END_REF] , which are lipids specially designed for give a desired set of properties or functional behaviour. Despite this incentive, the search for new blends of oils designed for a specific application is still based upon empirical experience rather than upon scientific background [START_REF] Ghotra | Lipid shortenings: a review[END_REF] . Vegetable oils industry can better take advantages from market opportunities if a systematic procedure for search the best mixtures/molecules that match the desired properties is available.

This work presents some steps in the development of a computational framework (Virtual Laboratory) for finding the optimal lipid molecules or mixtures for a given application. It is based on the concept of Computer Aided Product Design (CAPD) that can be divided in two sub-problems: Computer-Aided Molecular Design (CAMD) and Computer-Aided Mixture/blend Design (CAM b D) [START_REF] Gani | Chemical product design: challenges and opportunities[END_REF] . CAPD reverse principle is: given a set of desired properties, find the best molecules/mixtures that satisfy the needs. According to Constantinou et al. [START_REF] Constantinou | Computer Aided Product Design: Problem formulations, methodology and applications[END_REF] , the solution of the CAPD reverse problem can find new and stable products minimizing the cost of their development while producing superior products. This approach was successful for applications in others fields: drugs, solvents and refrigerants substitutions, polymers, paints and coatings [START_REF] Constantinou | Computer Aided Product Design: Problem formulations, methodology and applications[END_REF] . We propose to apply it to the design of tailor-made lipids, mainly triacylglycerol molecules/mixtures with improved chemical, physical, medical, nutritious and subjective properties, identifying the best raw materials for achieving it. For a tailor-made lipids virtual laboratory, CAPD concepts need extension to a multi-level modelling to address all the relevant tailor-made lipids length scales, ranging from fatty acids molecular structures to blends of oils solid -liquid properties.

MULTI-LEVEL MODELLING OF STRUCTURED LIPIDS

Structured lipids are modified triacylglycerol molecules by chemically or enzymatically catalyzed reactions. 95% of vegetable oils content is made of hundreds of different triacylglycerols (TAGs); which are molecules formed by 3 fatty acids sterified in a glycerol backbone. For designing specific TAGs molecules or mixtures with improved properties, factors as the chain length of the three fatty acids, their relative position and the number and location of double bonds must be considered at the molecule level. At the mixture level, composition and type of oils used need to be treated as optimisation variables. In Figure 1 Table 1 summarizes the most important natural higher fatty acids. It highlights the large number of possible molecules; a further motivation for a systematic computational procedure able to search the best candidates that fit the desired properties given by the user.

STRUCTURED LIPIDS PROPERTIES FOR INDUSTRIAL APPLICATIONS

Figure 2 shows the exponential increase in the number of possible molecules according to the number of fatty acids allowed, taking into account all types of isomers, including optical ones. Fatty acids type (chain size and double bonds) and position in a TAG structure (Table 1) directly influence thermodynamic and transport properties, solid solubility, crystallization and melting behaviour, solid content, melting profiles, subjective properties and are also related to biochemical reactions in absorption, transport and metabolism of fats in human body [START_REF] Mu | The metabolism of structured triacylglycerols[END_REF] [START_REF] Karupaiah | Effects of stereospecific positioning of fatty acids in triacylglycerol structures in native and randomized fats: a review of their nutritional implications[END_REF] .

In addition to these intramolecular factors, intermolecular issues, in particular polymorphism affects strongly TAGs properties. TAG molecules present basically three different polymorphisms (α, β', β), each one corresponding to a solid crystalline arrangement [START_REF] Sato | Crystallization behaviour of fats and lipids -a review[END_REF] . Differences in molecular packing affect important properties related to crystallization and melting behaviour, such as melting points and enthalpy of fusion (molecule level), melting ranges and solid content (mixture level) as well as subjective properties related to fat crystal networks in edible products. In fat-based food industry, the quality assessment for some products is attached to a particular desired polymorphism. For example, in margarines, β crystals should be avoided (they cause graininess), in other hand β polymorphism is desired in chocolate production [START_REF] Osborn | Structured Lipids-Novel Fats with Medical, Nutraceutical, and Food Applications[END_REF] . The polymorphisms dependent thermal properties play a very important role in quality of edible products and directly influence melting behaviour and related perception of texture [START_REF] Engelen | Oral physiology and texture perception of semisolids[END_REF] . Many subjective and functional properties can also be related to these properties: like texture, creaminess, flavour release, graininess and spreadability.

PROPERTY MODELS

Group Contribution Models (GC)

Group Contribution Methods are predictive properties models [START_REF] Marrero | Group-contribution based estimation of pure component properties[END_REF] well suited for CAPD problems provided that they are correlated to data relevant to the problem. Basically, for a given property, a set of known experimental data for different molecules is used to determine correlations where the parameters are the contribution of each chemical group present in the molecule. The equation (1) represents the general equation for property estimations methods using up to third-order groups [START_REF] Marrero | Group-contribution based estimation of pure component properties[END_REF] :

      j k k k j j i i i E O z D M w C N ) X ( f (1)
where N, M and O are the number of times a first, second or third order group, respectively, appears in the molecule, C, D and E are the contributions for the property and w and z are binary (0,1) variables used to include or not second and third order groups (respectively).

In this work, the GC method of [START_REF] Marrero | Group-contribution based estimation of pure component properties[END_REF] is used for computing the normal boiling point (Tb), critical properties (T c ,V c ,P c ), standard Gibbs energy at 298 K (ΔG f ), standard enthalpy of formation at 298 K (ΔH f ) and standard enthalpy of vaporization at 298 K (ΔH v ). The groups from [START_REF] Marrero | Group-contribution based estimation of pure component properties[END_REF] used to represent any type of triacylglycerols are given in Table 2. The type and number of these groups in TAGs vary according to TAG size and chain unsaturations and are automatically managed in the virtual laboratory.

Table 2: Groups used to model Triacylglycerols molecules.

Group Order

CH 3 First CH 2 First CH First CH=CH First CH 2 COO First CH 2 -CHm=CHn (m,n in 0…2) Second COO-CHn-CHm-OOC (n,m in 1,2) Second

Group Interaction Contribution Models (GIC)

To take into account isomerism, well acquainted in TAGs, Group Interaction Contribution Models (GIC) are also used in the virtual laboratory development, since the interactions parameters for a given group is available. If there are no available parameters for a chain interaction, models by [START_REF] Wesdorp | Liquid -multiple solid phase equilibria in fats[END_REF] can also cope with isomers and symmetry in molecules.

An improved GIC model was presented by [START_REF] Zéberg-Mikkelsen | Predicting the melting points and the enthalpies of fusion of saturated triglycerides by a group Contribution method[END_REF] , and it is used for computing temperature and enthalpy of fusion. This model is based on the following equations:

   IJ IJK , T m T K T  ( 2 )      IJ IJK , H m H K H ( 3 )
Where K is a parameter related to the type of TAG according to the fatty acid repartition (III, IIJ, IJI or IJK) and T ij and H ij are group interaction parameters for temperature and enthalpy of fusion respectively.

RESULTS

Table 3 shows some results for thermodynamic properties calculated using the GC method presented in Equation 1. The virtual laboratory is able to cope with more than 15,600 different TAGs. Knowledge of these fundamental properties through the virtual laboratory environment avoids experimental works in the earlier stages of finding new molecules or mixture for a given application, especially when there is a large possibility of combining the fatty acids on the glycerol backbone. Besides, GC methods enable to evaluate molecules for which thermodynamic experimental data is not available. At some latter step within the CAMD framework, experimental data acquisition will then be made on the sole valid candidate molecules.

Thermal Properties

Table 4 shows the ability of GIC models to cope with isomers. Data for three pairs of isomers in each polymorphisms are showed to illustrate the computing results using the model of [START_REF] Zéberg-Mikkelsen | Predicting the melting points and the enthalpies of fusion of saturated triglycerides by a group Contribution method[END_REF] , which show and absolute average deviation of 0.36 % for melting points and 2.17 % for enthalpy of fusion.

Viscosity

Fats show a complex rheological behaviour due to a mixture of solid crystalline and liquid oil. They have a pseudoplastic behaviour: they look like a solid but can be made to flow by an external force. Under small stresses elastic and viscoelastic deformations occur; under larger stress, yielding and flow occur [START_REF] Walstra | Crystallization Processes in Fats and Lipid Systems[END_REF] . Prediction of pure molecules viscosity gives an understanding of relations between fatty acids structure and viscosity, highlighting the CAMD/CAMbD search procedure.

Rheological properties influence many functional and subjective properties related to fat-based edible products. They can indeed impact flavour perception [START_REF] Bayarri | The role of rheological behaviour in flavour perception in model oil/water emulsions[END_REF] .

To predict the viscosity for pure triacylglycerol molecules, the present work uses the group contribution method of [START_REF] Joback | Estimation of Pure Component Properties from Group-contributions[END_REF] and some results are presented in Figure 3. The viscosity of 5 TAGs vs temperatures is used to illustrate some structure-property relationships. Comparing the molecules, it can be observed that double bounds have a major influence in viscosity. The linoleic acid -based TAG (LinLinLin), despite its bigger chain (18 C) shows the lowest viscosity, due to its higher degree of unsaturation. For TAGs with the same size of fatty acids chain, the higher the number of double bounds, the lower the viscosity.

Although these relations are known, the virtual laboratory gives quantitative data for this structure-property relationship. It also gives a predictive capacity for molecules where no experimental data exist or are difficult to measure.

In Figure 3, the lower temperature used (T = 400 K) was kept above the melting point of β modification for all TAGs, to ensure that all phases are liquid, as β polymorphism has the highest melting points among the polymorphisms. In general, natural unsaturated fatty acids contain cis double bounds. These bonds produces a bend in the straight carbon chain [START_REF] Ghotra | Lipid shortenings: a review[END_REF] that leads to changes in the way molecules are packing together, reducing temperature of fusion and viscosity. Generation of trans double bounds are often related to industrial process of vegetable oils.

CONCLUSIONS

Property prediction is an essential step in the implementation of a Computer Aided Molecular and Mixture Design framework in lipids systems. Group contribution methods can give a predictive capacity for properties needed for modelling structured lipids. They allow insights of how molecular structure affects thermodynamic properties and allow a further investigation for subjective and functional properties of fat-based products. For properties related to crystalline structure, interaction between groups and polymorphisms must be taken into account to a better accuracy in results.

A predictive capacity for different types of properties is also useful in the design, operation and optimization of unit processes, as distillation, heat exchangers, piping and reactors in vegetable oils industry, in a Computer Aided Process Design context.

Besides its use in edible products production, the development of the present framework has also a potential usefulness for other industries that use lipid-based raw materials, like cosmetic, pharmaceuticals, lubricants and biofuels.

The move towards the molecular level for product design is a major trend and challenge for Process System Engineering (PSE) in the future [START_REF] Grosmann | Research Challenges in Process Systems Engineering[END_REF] . Development of predictive capabilities for properties of compounds and mixtures and generation of alternatives in order to apply optimization methodologies developed previously for decisions in process design are needed. We propose that the use of renewable materials as vegetable oils is well suited for this structured molecular approach by using Computer Aided Molecular and Mixture Designs techniques according to this PSE Process and Product Design topic.

Further implementation of mixture property models to be coupled with optimization techniques able to find the best candidate molecular structures/mixtures for structured lipids design development.

Figure 1 :

 1 Figure 1: General multi-scale problem representation.

Figure 2 :

 2 Figure 2: Number of triacylglycerols that can be formed from fatty acids.

Figure 3 :

 3 Figure 3: Viscosity of pure triacylglycerols vs. temperature in liquid state.

Lipids: fatty acyls, glycerolipids, glycerophospholipids, sphingolipids, sterol lipids, prenol lipids, saccharolipids, polyketides -Properties: thermodynamic, transport. -Models: GC, GIC, others correlations. -Lipids: TAGs mixtures, vegetable oils, blends of vegetable oils -Properties: thermodynamic, transport, phase related (Solid Fat Content, clear point, softening point, viscosity, density) and subjective properties (aroma, texture etc). -Models: Solid-Liquid Equilibria, mixture rules, correlations.

  illustrates a general description of the problem.

	-

Table 1 : Natural fatty acids. Carbon Number: Double Bounds Structure (H3C-(R)-COOH) Trivial Names Code

 1

	10:0	-[CH 2 ] 8 -	Capric	Dec
	12:0	-[CH 2 ] l0 -	Lauric	Lau
	14:0	-[CH 2 ] 12 -	Myristic	My r
	16:0	-[CH 2 ] 14 -	Palmitic	Pam
	16:1	-[CH 2 ] 5 CH=CH[CH 2 ] 7 -	Palmitoleic	DPam
	18:0	-[CH 2 ] 16 -	Stearic	Ste
	18:1(9)	-[CH2]7CH=CH[CH2]7-	Oleic	Ol e
	18:1(11)	-[CH 2 ] 5 CH=CH[CH 2 ] 9 -	Vaccenic	Vac
	18:2(9,12)	-[CH2]3(CH2CH=CH)2[CH2]7-	Linoleic	Lin
	18:3(9,12,15)	-(CH 2 CH=CH) 3 [CH 2 ] 7 -	(9,12,15)-Linolenic	aLnn
	18:3(6,9,12)	-[CH 2 ] 3 (CH 2 CH=CH) 3 [CH 2 ] 4 -	(6,9,12)-Linolenic	gLnn
	18:3(9,11,13)	-[CH 2 ] 3 (CH=CH) 3 [CH 2 ] 7 -	Eleostearic	eSte
	20:00	-[CH2]18-	Arachidic	Ac	h
	20:2(8,11)	-[CH 2 ] 6 (CH 2 CH=CH) 2 [CH 2 ] 6 -	---	Δ2Ach
	20:3(5,8,11)	-[CH 2 ] 6 (CH 2 CH=CH) 3 [CH 2 ] 3 -	---	Δ3Ach
	20:4(5,8,11,14)	-[CH 2 ] 3 (CH 2 CH=CH) 4 [CH 2 ] 3 -	Arachidonic	Δ4Ach
	22:0	-[CH 2 ] 20 -	Behenic	Beh
	24:0	-[CH 2 ] 22 -	Lignoceric	Lig
	24:1	-[CH 2 ] 7 CH=CH[CH 2 ] 13 -	Nervonic	Ner
	26:00	[CH 2 ] 24 -	Cerotic	Crt
	28:0	-[CH 2 ] 26 -	Montanic	Mon

Table 3 : Thermodynamic properties of triacylglycerols (GC Method).

 3

	TAG molecule Tc (K) Pc (bar)	Vc (cm3/mol)	ΔGf (kJ/mol)	ΔHf (kJ/mol)	ΔHv (kJ/mol) Tb (K)
	StePamMyr 1017.47	7.54	2951.99	-579.84 -2067.56	287.36	830.48
	LauCrtPam	1039.12	7.27	3289.67	-531.45 -2192.53	316.82	852.22
	PamLigPam 1045.91	7.19	3402.23	-515.33 -2234.19	326.64	859.02
	BehCrtBeh	1088.56	6.81	4190.15	-402.43 -2525.80	395.38	901.55
	AchSteLau	1024.92	7.44	3064.55	-563.71 -2109.22	297.18	837.97
	AchVacPam 1040.31	7.27	3275.82	-452.91 -2078.31	316.34	852.34
	PamLinPam 1027.44	7.46	3036.85	-406.63 -1880.78	296.21	838.21
	OleOleOle	1042.67	7.29	3248.12	-295.83 -1849.88	315.37	852.57
	MyrOlePam 1018.78	7.55	2938.14	-501.30 -1953.34	286.88	830.61
	LigOleLau	1040.31	7.27	3275.82	-452.91 -2078.31	316.34	852.34

Table 4 : Fusion enthalpy and temperature of triacylglycerols (GIC Method).

 4

	TAG molecule	α	β'	β'
	ΔH fus (kJ/mol)
	StePamMyr	89.59 116.10 141.07
	PamSteMyr	94.42 135.71 130.52
	LauPamLau	79.56 125.57 142.15
	PamLauLau	65.10 90.00 115.17
	MyrSteLau	71.48 100.23 108.46
	SteMyrLau	81.61 93.21 127.86
	T fus (K)
	StePamMyr	315.01 330.53 332.26
	PamSteMyr	312.48 331.96 331.69
	LauPamLau	292.64 315.36 319.48
	PamLauLau	295.36 314.87 318.77
	MyrSteLau	299.21 320.05 322.80
	SteMyrLau	304.41 322.71 324.49


Armand Jordana

Avadesh Meduri

Etienne Arlaud

Justin Carpentier

Ludovic Righetti

Risk-Sensitive Extended Kalman Filter

à 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.

Risk-Sensitive Extended Kalman Filter

Armand Jordana, Avadesh Meduri, Etienne Arlaud, Justin Carpentier and Ludovic Righetti Abstract-In robotics, designing robust algorithms in the face of estimation uncertainty is a challenging task. Indeed, controllers often do not consider the estimation uncertainty and only rely on the most likely estimated state. Consequently, sudden changes in the environment or the robot's dynamics can lead to catastrophic behaviors. In this work, we present a risk-sensitive Extended Kalman Filter that allows doing outputfeedback Model Predictive Control (MPC) safely. This filter adapts its estimation to the control objective. By taking a pessimistic estimate concerning the value function resulting from the MPC controller, the filter provides increased robustness to the controller in phases of uncertainty as compared to a standard Extended Kalman Filter (EKF). Moreover, the filter has the same complexity as an EKF, so that it can be used for real-time model-predictive control. The paper evaluates the risk-sensitive behavior of the proposed filter when used in a nonlinear model-predictive control loop on a planar drone and industrial manipulator in simulation, as well as on an external force estimation task on a real quadruped robot. These experiments demonstrate the abilities of the approach to improve performance in the face of uncertainties significantly.

I. INTRODUCTION

Adapting the decisions robots make based on their perception of the world is key to deploying robots outside factories. More precisely, controllers should adapt to the degree of certainty or confidence of the robot's belief of the world. For instance, it is important that a quadruped chooses conservative footholds and slows body movements when its confidence in the location of the ground decreases. Robust output feedback Model Predictive Control (MPC) studies methods that can adapt robot decisions based on the confidence of the perception module. However, the general nonlinear problem is very difficult, and practical algorithms remain often limited to linear systems [START_REF] Mayne | Model predictive control: Recent developments and future promise[END_REF].

The common practice in robotics is to decouple estimation and control (i.e., assume that the certainty equivalence principle holds) [START_REF] Kuindersma | Optimization-based locomotion planning, estimation, and control design for the atlas humanoid robot[END_REF]- [START_REF] Dantec | Whole-body model predictive control for biped locomotion on a torque-controlled humanoid robot[END_REF]. This approach is often chosen due to the availability of separate and tractable control and estimation algorithms that can be deployed on the robot. The estimation module is often a variation of a Gaussian filter, such as an Extended Kalman Filter (EKF) [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF], which computes both This work was in part supported by the National Science Foundation grants 1932187, 1925079, 2026479, 2222815 and the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute), the European Union through the AGIMUS project (GA no.101070165) and the Louis Vuitton ENS Chair on Artificial Intelligence.

Armand Jordana, Avadesh Meduri and Ludovic Righetti are with the Tandon School of Engineering, New York University, Brooklyn, NY. (email: aj2988@nyu.edu, am9789@nyu.edu, lr114@nyu.edu).

Justin Carpentier and Etienne Arlaud are with Inria, Département d'informatique de l'ENS, École normale supérieure, CNRS, PSL Research University, Paris, France (e-mail: justin.carpentier@inria.fr).

the mean and uncertainty of the state estimates from sensor information. In control, an increasingly popular approach is MPC, which consists in solving an optimal control problem numerically, e.g., using Differential Dynamic Programming (DDP) [START_REF] Neunert | Whole-body nonlinear model predictive control through contacts for quadrupeds[END_REF], [START_REF] Erez | An integrated system for real-time model predictive control of humanoid robots[END_REF]- [START_REF] Kleff | Highfrequency nonlinear model predictive control of a manipulator[END_REF], at each time step or at a fixed frequency. The controller can then adapt its behavior online based on the current robot and environment states. During deployment, the estimation module is used to compute the mean of the state estimate, which is then passed on to the controller to compute the optimal behavior [2]- [START_REF] Dantec | Whole-body model predictive control for biped locomotion on a torque-controlled humanoid robot[END_REF]. Unfortunately, relying on the most likely outcome can lead to catastrophic behavior. For instance, on a load-carrying task with a quadruped where the load's mass is unknown, the notion of mean might not be appropriate as this could lead the quadruped to apply insufficient force on the ground and then fall.

Some approaches try to address this issue by adding robustness or safety bounds in either the estimation or control block while keeping them independent. For instance, Robust Extended Kalman filtering [START_REF] Einicke | Robust extended kalman filtering[END_REF] adds robustness to inaccuracies of the EKF or the model. However, the control objective is disregarded, and therefore the controller cannot be robust to estimation uncertainties. Robust MPC has been studied and applied to robots. [START_REF] Villa | Model predictive control of biped walking with bounded uncertainties[END_REF] used tube-based MPC to control a biped. [START_REF] Gazar | Stochastic and robust mpc for bipedal locomotion: A comparative study on robustness and performance[END_REF] used linear stochastic MPC to account for uncertainties in bipedal walking. However, this line of work assumes the state to be known. In contrast to such approaches, we aim to link estimation and control by adding into the estimation module a notion of control performance to improve robustness to the estimation uncertainty.

In this work, we leverage our previous theoretical results on dynamic game control with imperfect state observation [START_REF] Jordana | Stagewise newton method for dynamic game control with imperfect state observation[END_REF] to introduce the Risk-Sensitive Extended Kalman Filter (RS-EKF), a novel filter that enables online risk-sensitive output feedback MPC. The RS-EKF computes state estimates robust to measurement uncertainty while taking into account the value function provided by the controller, i.e., the estimator tailors risk reduction to the control objectives. This, in turn enables automatic modification of robot decisions to be cautious in times of high environmental perturbation. Furthermore, RS-EKF has a similar computational cost to an EKF, allowing real-time deployment. To demonstrate the ability of the filter, we use it together with a DDP-based online non-linear controller to perform risk-sensitive outputfeedback MPC on various simulated robots, such as a quadrotor subjected to arbitrary changes in its mass, and a KuKa robot facing unforeseen environmental disturbances. Finally, we test the filter on a real quadruped robot Solo12 [START_REF] Grimminger | An open torque-controlled modular robot architecture for legged locomotion research[END_REF] to perform an external force estimation and balancing task. These experiments demonstrate that the robots are more robust to perturbations with the RS-EKF algorithm than a classical EKF. To the best of our knowledge, this is the first time that a non-linear risk-sensitive output-feedback MPC controller has been deployed on a robot.

II. BACKGROUND

A. Dynamic game output feedback MPC

To design a controller sensitive to the risk related to estimation uncertainty, Whittle [START_REF] Whittle | Risk-sensitive linear/quadratic/gaussian control[END_REF] introduced a zero-sum game that aims at solving jointly the estimation and control problem. Given a history of measurements y 1:t , a history of control inputs, u 0:t-1 and a prior on the initial state x0 , we aim to find a control sequence u t:T -1 that minimizes a given cost ℓ while an opposing player aims to find the disturbances (w 0:T , γ 1:t ) that maximize this cost ℓ minus a weighted norm of the disturbances. Such a problem is formally written as:

min u t:T -1 max w 0:T max γ1:t ℓ T (x T ) + T -1 i=0 ℓ i (x i , u i ) (1)
- 1 2µ   ω T 0 P -1 ω 0 + t j=1 γ T j R -1 j γ j + T j=1 w T j Q -1 j w j   subject to x 0 = x0 + w 0 , (2a)
x j+1 = f j (x j , u j ) + w j+1 , 0 ≤ j < T, (2b)
y j = h j (x j ) + γ j , 1 ≤ j ≤ t. (2c)
where µ > 0. x j is the state, ω j the process disturbance, y j the measurement,γ j the measurement disturbance, T the time horizon, t the current time. The transition model f j , the measurement model h j , and the cost ℓ j are assumed to be C 2 . The measurement uncertainty R j , the process uncertainty Q j and the initial state uncertainty P are positive-definite matrices. This weighted sum of the disturbances can be seen as the estimation of maximum a posteriori probability (MAP) under Gaussian assumption. Hence, R j , Q j and P can be thought of as the covariance matrix of Gaussian noise. Interestingly, this problem encompasses both formulations of control and estimation. If t = 0, in the limit where µ tends to zero, we find the generic OCP formulation [START_REF] Campi | Nonlinear discrete-time risk-sensitive optimal control[END_REF] which directly minimizes the cost function assuming standard deterministic dynamics. And, if t = T and if we consider all costs ℓ i to be null, then, (1) is equivalent to maximizing the estimation maximum a posteriori probability (MAP). Here, the parameter µ is referred to as the risk-sensitive parameter and regulates how adversarial the problem is.

Whittle [START_REF] Whittle | Risk-sensitive linear/quadratic/gaussian control[END_REF] provided a solution to this min-max problem in the case of linear dynamics and quadratic costs. This solution can be iteratively obtained with two recursions, one on past disturbances and one on future disturbances. Recently, [START_REF] Jordana | Stagewise newton method for dynamic game control with imperfect state observation[END_REF] showed how this solution could be used to implement an efficient Newton's method that iteratively searches for a saddle point of the more general Problem [START_REF] Mayne | Model predictive control: Recent developments and future promise[END_REF]. Exploiting the sparsity of a problem, the proposed Newton step has a linear complexity in the time horizon (note that naive optimization has at least a cubic complexity [START_REF] Copp | Simultaneous nonlinear model predictive control and state estimation[END_REF]). The solution can then be interpreted as a risk-sensitive Kalman smoother coupled to minimax DDP [START_REF] Morimoto | Minimax differential dynamic programming: An application to robust biped walking[END_REF].

In this work, we use these insights to derive a computationally efficient risk-sensitive extended Kalman filter that can be used for output feedback MPC. This is done by simplifying Problem (1) to match assumptions common to EKF and DDP while keeping the adversarial min-max formulation.

B. Extended Kalman Filter

We now discuss the structure of the EKF necessary to derive our filter. The EKF is usually derived by computing the probability a posteriori of the state given measurements, using the linearized dynamics and a Gaussian noise assumption [START_REF] Thrun | Probabilistic robotics[END_REF]. However, the EKF can also be derived from an optimization point of view [START_REF] Bell | The iterated kalman filter update as a gauss-newton method[END_REF]. More precisely, the EKF can be seen as a Gauss-Newton step around a well-chosen point on the log-likelihood of the maximum a posteriori probability (MAP), i.e. log(p(x t , x t-1 |y t )) [START_REF] Bell | The iterated kalman filter update as a gauss-newton method[END_REF]. Assuming disturbances follow Gaussian distributions, γ t ∼ N (0, R t ), ω t ∼ N (0, Q t ), the MAP can be written as:

max xt,xt-1 -(y t -h t (x t )) T R -1 t (y t -h t (x t )) -(x t -f t-1 (x t-1 , u t-1 )) T Q -1 t (x t -f t-1 (x t-1 , u t-1 )) -(x t-1 -xt-1 ) T P -1 t-1 (x t-1 -xt-1 ) (3)
where xt-1 is the prior knowledge on the past state and P t-1 its associated covariance matrix. As shown in [START_REF] Bell | The iterated kalman filter update as a gauss-newton method[END_REF], a Gauss-Newton step around xt-1 and xt = f t-1 (x t-1 , u t-1 ) on ( 3) leads to well-known recursion [START_REF] Thrun | Probabilistic robotics[END_REF]:

Pt = Q t + F t-1 P t-1 F T t-1 (4)
K t = Pt H T t (R t + H t Pt H T t ) -1 (5) P t = (I -K t H t ) Pt (6) μt = K t (y t -h t (x t )) (7) xt = xt + μt (8)
where

F t-1 = ∂ x f t-1 (x t-1 , u t-1 ) and H t = ∂ x h t (x t ).
Here, xt is the most likely estimate and P t is the covariance uncertainty. In practice, at the next time step, xt is used as the prior knowledge on the state. We notice the similar structure of the costs of Problem (1) and Eq. ( 3), except that the EKF only uses one measurement and does not include the control cost ℓ j . Hence, Eq. ( 1) can be seen as a maximization of the estimation log-likelihood up to some cost terms. We will leverage this similarity to derive a risk-sensitive version of the EKF. More precisely, we will add cost-dependent terms in the maximization (3) in order for the filter to adapt to the control objective.

C. Nonlinear MPC

Once the state is estimated, it can then be used in the control module. In this work, we focus on controllers defined as the solution to a minimization problem. At each time step, a stagewise cost defined over a horizon H is minimized over future control inputs u t , u t+1 , . . . , u H-1 :

L t (u t , . . . , u H-1 ) = ℓ t+H (x T +H ) + t+H-1 j=t ℓ j (x j , u j ) (9)
where the state sequence is implicitly defined by the dynamics x j+1 = f j (x j , u j ) and where x t is assumed to be the state estimated by the filter. At each time step, only the term u t is applied to the real system. In the next time step, the state estimate is updated given the new measurement and the optimal control problem is solved again. There are various ways to solve efficiently this minimization problem. In robotics, a popular algorithm is DDP [START_REF] Murray | Differential dynamic programming and newton's method for discrete optimal control problems[END_REF] which reassembles the Newton method but with linear complexity in the time horizon. Additionally, DDP provides a quadratic approximation of the value function which we exploit in our derivation of the risk-sensitive filter.

III. RISK SENSITIVE FILTER

We now introduce RS-EKF, which builds on the dynamic game defined in Equation [START_REF] Mayne | Model predictive control: Recent developments and future promise[END_REF]. First, we modify the game to account for typical assumptions made for MPC while keeping the adversarial part that provides the risk-sensitive behavior. Then, we show how the solution can be computed with a Gauss-Newton step similar to the EKF, leading to an algorithm of similar complexity. Our formulation leads to a modified update in the filter of which the standard EKF can be seen as a limit case.

First, as for the EKF, we consider a history of measurements of length one. Furthermore, we disregard future uncertainties and assume deterministic dynamic equations for the future as is done in classical MPC formulations. Indeed, we expect that the high-frequency re-planning will compensate for model discrepancies. In the end, the problem is supposed to be adversarial only with respect to the uncertainties related to the estimation. This can be written as:

min u t:t+H-1 max wt max wt-1 max γt L t (u t , . . . , u H-1 ) (10)
- 1 2µ γ T t R -1 t γ t + w T t Q -1 t w t + w T t-1 P -1 t-1 w t-1 s.t. x t-1 = xt-1 + w t-1 , (11a) x t = f t-1 (x t-1 , u t-1 ) + w t , (11b)
y t = h t (x t ) + γ t . (11c) x j+1 = f j (x j , u j ), t < j < T. (11d)
As presented in [START_REF] Jordana | Stagewise newton method for dynamic game control with imperfect state observation[END_REF], one of the key features of the dynamic game is that some of the constraints can be removed with an appropriate change of variable. Indeed, we can use the equality constraints of Equations (11a), (11b) and (11c) to replace the disturbance maximization into a maximization over x t-1 , x t :

min u t:t+H-1 max xt-1,xt L t (u t , . . . , u H-1 ) (12)
- 1 2µ (y t -h t (x t )) T R -1 t (y t -h t (x t )) - 1 2µ (x t -f t-1 (x t-1 , u t-1 )) T Q -1 t (x t -f t-1 (x t-1 , u t-1 )) - 1 2µ (x t-1 -xt-1 ) T P -1 t-1 (x t-1 -xt-1 ) subject to x j+1 = f j (x j , u j ), t < j < T,
By definition of the MAP [START_REF] Thrun | Probabilistic robotics[END_REF], this can be written:

min u t:t+H-1 max xt-1,xt L t (u t , . . . , u H-1 ) - 1 µ log(p(x t , x t-1 |y t )) subject to x j+1 = f j (x j , u j ), t < j < T. (13)
Problem ( 13) is intractable in the general case. However, by taking the concave-convex assumption, the minimization and maximization can be interchanged according to the minimax theorem. Consequently, the problem is equivalent to:

max xt-1,xt log(p(x t-1 , x t |y t )) + µV t (x t ), (14)
where V t is the value function of the OCP:

V t (x t ) = min u t:t+H-1 L t (u t , . . . , u H-1 ) (15)
Note that in the simplification from Eq. ( 1) to Eq. ( 10), it is not necessary to disregard future uncertainties as the value function could be the one resulting from minimax DDP [START_REF] Morimoto | Minimax differential dynamic programming: An application to robust biped walking[END_REF].

If µ = 0, we will obtain the unbiased estimate of Kalman filtering and the estimate will be independent of the control objective. Otherwise, if µ > 0, the term µV (x t ) will bias the estimate towards regions with higher value function, which in turn will force the controller to be more conservative.

We now take a Gauss-Newton step on the objective of Eq. ( 14) around the prior: xt-1 and xt = f t-1 (x t-1 , u t-1 ). V t (x t ) is independent of x t-1 therefore, as shown in the appendix, the maximization over x t-1 can be simplified to:

max xt - 1 2 (x t -xt ) T P -1 t (x t -xt ) (16)
+ µ 1 2 (x t -xt ) T V xx t (x t -xt ) + µ(x t -xt ) T v x t
where xt and P t are defined as in Eq. ( 8) and [START_REF] Dantec | Whole-body model predictive control for biped locomotion on a torque-controlled humanoid robot[END_REF]. where V xx t (respectively v x t ) is the hessian (respectively the gradient) of the value function. Those are typically provided by optimal control algorithms such as DDP. In the end, the solution on the maximization over x t is:

xRS t = xt + (I -µP t V xx t ) -1 (μ t + µP t v x t ) (17)
Interestingly, if µ = 0, we recover the EKF. This was to be expected as, when µ tends to zero, the solution of problem ( 1) is exactly the solution of the neutral case where estimation and control are solved independently [START_REF] Whittle | Risk-sensitive linear/quadratic/gaussian control[END_REF]. Otherwise, the estimate is shifted towards regions with higher cost values. Importantly, the magnitude of the shift depends on P t the covariance matrix of the estimation. Note that µ cannot be arbitrarily large as (I -µP t+1 V xx t+1 ) needs to be positive definite. Larger values of µ would make the minmax problem defined in Eq. (1) ill-posed. More details on this limit value can be found in [START_REF] Whittle | Risk-sensitive linear/quadratic/gaussian control[END_REF]. In the end, the estimate is shifted towards P t v x t , i.e. towards a region with a larger cost function, and the magnitude of this shift is increased in the direction corresponding to large eigenvalues of P t V xx t . We obtained the solution to the maximization problem [START_REF] Villa | Model predictive control of biped walking with bounded uncertainties[END_REF]. Therefore, the cost function can now be minimized with respect to the control inputs by taking xRS t as an initial condition of the optimal control problem, which can be solved for example with DDP.

Algorithm 1 summarizes the estimation procedure. It can then be used to do output-feedback MPC efficiently. At each time step, given a measurement, past control input, and a quadratic approximation of the value function, a risksensitive estimate can be computed. This estimate is then used to minimize the cost function [START_REF] Koenemann | Whole-body model-predictive control applied to the hrp-2 humanoid[END_REF] for MPC and the first control input is applied to the real system. Lastly, the quadratic approximation of the value function at t + 1 is saved as it will be used at the next estimation step. IV. EXPERIMENTS This section presents simulation and real robot experiments to illustrate the benefits of the proposed algorithm and demonstrate its applicability to real problems. We study three test problems where we deploy the RS-EKF inside a MPC loop: a planar quadrotor with a load estimation task where we demonstrate qualitatively that the RS-EKF can bring conservatism appropriately in phases of high uncertainty, a push-recovery experiment on a 7-dof industrial manipulator on which we perform a quantitative study and lastly, an external force estimation task on a real quadruped robot to showcase the viability of the method on real systems. In each experiment, we use the DDP implementation provided by Crocoddyl [START_REF] Mastalli | Crocoddyl: An efficient and versatile framework for multi-contact optimal control[END_REF] to solve the optimal control problem given the filter estimate. All the code is available online 1 .

Algorithm 1: Risk Sensitive EKF Input: xt-1 , u t-1 , y t , P t-1 , Q t , R t , V t , v t / * Predict * / 1 Pt ← Q t + F t-1 P t-1 F T t-1 2 xt ← f (x t-1 , u t-1 ) / * Classical Update * / 3 K t ← Pt H T t (R t + H t Pt H T T ) -1 4 P t ← (I -K t H t ) Pt 5 μt ← K t (y t -h t (x t )) / *

A. Planar quadrotor

In this first scenario, we consider a planar quadrotor executing a load-carrying task. The goal is to move the quadrotor from position (p x , p y ) = (0, 0) to position (1, 0) while carrying a load during the first half of the itinerary. The robot mass is 2 kg and the mass's load which is unknown a priori is 3 kg. The system dynamics is:

mp x = -(u 1 + u 2 ) sin(θ), mp y = (u 1 + u 2 ) cos(θ) -mg, md θ = r(u 1 -u 2 ), ( 18
)
where m is the mass of the robot, d the distance between the rotors, θ the orientation of the quadrotor. u 1 and u 2 are the control inputs representing the force applied at each rotor.

1 https://github.com/machines-in-motion/ risk-sensitive-EKF

In this experiment, we want to estimate online the mass parameter that changes in the middle of the flying phase. As it is standard in parameter identification [START_REF] Voss | Nonlinear dynamical system identification from uncertain and indirect measurements[END_REF], we augment the system's state with the unknown parameter and let it be estimated recursively by the filter (EKF or RS-EKF). The state of the system is thus: x = p x p y θ ṗx ṗy θ m T and it is assumed that ṁ = 0 up to some random Gaussian noise. The dynamics are integrated with an Euler scheme and a time step of 0.05. We consider that P 0 = 10 -4 I 7 , R = 10 -4 I 3 , and Q is a 7 × 7 diagonal matrix where all terms are equal to 10 -4 except the last one that we set to 2 to represent the uncertainty in the changes of the load. Lastly, we set µ = 4 × 10 -3 . A stationary cost function of the following form is considered:

ℓ(x, u) = α 1 ∥p x -p des x ∥ 2 + ∥p y -p des y ∥ 2 + α 2 ∥θ∥ 2 + α 3 ∥ ṗx ∥ 2 + ∥ ṗy ∥ 2 + ∥ θ∥ 2 + α 4 ∥u -ū∥ 2 (19)
where ū = mg 2 , mg

2

T and where: α 1 = 100, α 2 = 10, α 3 = 0.01 and α 4 = 0.1. We consider a horizon of 20 nodes and re-plan at each new measurement, i.e. every 0.05s. Furthermore, we only measure: y = p x p y θ T to illustrate the estimator capabilities. We simulate 4s with both output-feedback MPC controllers: the first one relying on the standard EKF estimate and the second one relying on the RS-EKF estimate.

Figure 1 shows the real mass variation and the estimates of both methods. It can be seen that the RS-EKF is more reactive when the load is added or dropped. The RS-EKF estimate spikes in the phases of uncertainty, which adds some robustness. Those spikes can be explained by the fact that the uncertainty increases on the components of the state that are important in the cost function. In other words, some of the eigenvalue of P t V xx t become larger in the phases of uncertainty, which augments the shift on the estimate as shown in Equation [START_REF] Whittle | Risk-sensitive linear/quadratic/gaussian control[END_REF]. It can be seen that the controller relying on the risk-sensitive estimate is more reactive. To evaluate both controllers, we compute the average Mean Square Error (MSE) relative to the reference trajectory. The MSE of RS-EKF is 0.0011, and the one of EKF is 0.0024. Hence, in that situation, the risk sensitivity brings a 54% improvement in tracking. Furthermore, the average cost along the trajectory is equal to 0.0569 for the RS-EKF-based controller, while it is equal to 0.0880 for the EKF-based controller, yielding a 35% improvement. This illustrates how a filter informed of the cost objective can improve the controller's performance. As it can be seen in both Figures 1 and2, RS-EKF introduces a steady state error. This is because the estimation uncertainty never goes to zero due to the nonzero diagonal term of the matrices Q and R. Intuitively, it makes sense that the controller that plans for the worst is slightly sub-optimal in the phase with no environment perturbation uncertainty. The risk-sensitive filter demonstrates its benefits when there are perturbations in the environment as this corresponds to phases where the most likely estimate might be far from reality. In the end, this example illustrates the advantages of having a risk-sensitive controller which is conservative only in phases with large environmental perturbations.

B. Kuka robot

In this example, we consider the 7-DoF torque-controlled KUKA LWR iiwa R820 14. The 14-dimensional state is composed of the joint positions and velocities. The control input is a 7-dimensional vector of the torque applied on each joint. The continuous dynamics and its analytical derivatives are provided by Pinocchio [START_REF] Carpentier | The pinocchio c++ library: A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives[END_REF]. The goal is to track a reference trajectory with the end effector. To do so, we consider the following cost:

ℓ k (x k , u k ) = 10 -2 ∥x k -x∥ 2 2 + 10 -4 ∥u k -ū(x k )∥ 2 2 + 10 2 ∥p target k -p(x k )∥ 2 2 (20) ℓ T (x T ) = 10 2 ∥p target T -p(x k )∥ 2 2 + 10 -2 ∥x k -x∥ 2 2 ,
x, the initial state, is used for the state regularization and is the concatenation of the initial configuration 0.1 0.7 0 -1 -0.5 1.5 0 T and a 7-dimensional zero vector corresponding to the velocity. ū(x k ) is the gravity compensation term given by the rigid body dynamics. p(x k ) is the end-effector position obtained through forward kinematics. p target k is defined such that the end effector follows trajectories forming a circle in the xy plane. We use a horizon of 20 collocation point with an integration step of 0.05s and re-plan at 500 Hz.

In this experiment, we aim to showcase that the risksensitive filter can bring conservatism on phases with large perturbations from the environment which we simulate with large forces applied on the end effector. For the measurement model, we assume that all the states are observed with high accuracy, therefore, we set R = P 0 = 10 -6 I 14 . However, to model the disturbances in the dynamics, we set Q = 10 -1 I 14 . Finally, we consider µ = 7.5 × 10 4 .

P EE z [m]
Unexpected push Fig. 3. End effector trajectory on a tracking task for both the EKF-MPC and the RS-EKF-MPC. An unexpected force is applied between 1s and 2s.

Figure 3 depicts the end effector trajectory for both controllers and their respective estimates. From time 1s to 2s, an external force of norm 80 is applied on the end-effector on the x and z direction. Both controllers are pushed away from the reference trajectory. However, the risk-sensitive estimate overestimates the distance between the reference and the end-effector. Consequently, the robot is more aggressive in its response, and the end-effector remains closer to the reference. This illustrates how taking a pessimistic estimate with respect to the cost improves the running cost. Note that both estimates are originally state estimates but are mapped to end effector space through the forward kinematics to draw Figure 3. The fact that those estimates are pessimistic with respect to the task originally defined in end effector space illustrates well how the method can handle nonlinear dynamics and cost functions. To validate the consistency of the filter, 10, 000 experiments are performed with random external forces. The timing and direction of the forces are uniformly sampled while the duration is fixed to 1s and the norm is fixed to 80. Figure 4 shows the median end effector error trajectory. In average, we obtain a 32% improvement in the MSE. Furthermore, the mean cost is 22% lower with the risk-sensitive filter.

C. Load estimation on a quadruped robot

In this experiment, we deploy the RS-EKF on a real 12-degree-of-freedom, torque-controlled quadruped robot -Solo12 [START_REF] Grimminger | An open torque-controlled modular robot architecture for legged locomotion research[END_REF]. We demonstrate the superior performance of the RS-EKF in estimating wrenches applied on Solo12 while it is standing. We generate the standing behavior on Solo12 using a non-linear MPC scheme. At each control cycle, we minimize a cost function using a centroidal model to compute the optimal forces and trajectory that keep the robot's base at a desired height and orientation. Additionally, we use an augmented state to estimate external forces-torques applied to the robot [START_REF] Rotella | Humanoid momentum estimation using sensed contact wrenches[END_REF]:

ċ = 1 m l (21)
l = mg + Mc i=1 F i + F ext (22) k = Mc i=1 (p i -c) × F i + τ ext ( 23
)
Ḟext = 0, τext = 0 ( 24
)
where m denotes the mass, M c the number of end effectors in contact. Each p i represents a contact location. Here, the state is x = c l k F ext τ ext T which includes the Center of Mass (c), linear momentum (l), Angular Momentum (k) and External Force-torques (F ext , τ ext ). The measurement is y = c l k T up to some noise. A motion capture system measures the base position, velocity, and orientation and an IMU gives the orientation velocity. Joint encodings are provided and their velocities are derived with finite differences. Then, given, q, q, we can compute c, l, k with Pinocchio [START_REF] Carpentier | The pinocchio c++ library: A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives[END_REF], [START_REF] Carpentier | Pinocchio: fast forward and inverse dynamics for poly-articulated systems[END_REF] and can then be used as a measurement by the centroidal-based filter.

The control input,

u = F 1 . . . F Mc T , is a M c × 3
dimensional vector, representing the force applied at each end effector. For this experiment, the robot is standing, therefore, M c = 4. The cost function for the OCP is:

ℓ t (x, u) = (x -x ⋆ ) T H x (x -x ⋆ ) + (u -u ⋆ ) T H u (u -u ⋆ ) + 10 5 Mc i=1 ℓ barrier (u 3i ) ℓ T (x) = (x -x ⋆ ) T H x (x -x ⋆ ) (25)
where H x = BlockDiag(10 2 I 3 , 10I 6 ) and H u is a diagonal matrix where the diagonal terms are made of M c times the following sequence (10 -4 , 10 -4 , 10 -6 ). Lastly,

ℓ barrier (u) =    u 2 if u < 0 0 if 0 < u < 10 (u -10) 2 if u > 10 (26)
Here, x ⋆ is designed to keep a constant height above the ground and to keep the base horizontal, u ⋆ is the gravity compensation. The reference desired angular momentum for the OCP is adapted to bring the base back to a horizontal position. We do this by computing k * = 1 T log 3 (R t R T des ), where R t is the current base rotation matrix, R des the rotation matrix corresponding to the quaternion q des = [0, 0, 0, 1] and T the horizon length. The log 3 is a mapping from SE(3) to se(3). The nominal angular momentum aims to bring the base back to the desired orientation over the time horizon [START_REF] Meduri | Biconmp: A nonlinear model predictive control framework for whole body motion planning[END_REF]. The ℓ barrier is a quadratic barrier function that creates a soft constraint on the maximum forces the robot can apply on the ground.

We solve this OCP at 100 Hz using Croccodyl [START_REF] Mastalli | Crocoddyl: An efficient and versatile framework for multi-contact optimal control[END_REF] and then track the desired forces using a task space impedance dynamics QP [START_REF] Herzog | Momentum control with hierarchical inverse dynamics on a torque-controlled humanoid[END_REF] that we solve at 1 kHz using [START_REF] Bambade | Prox-qp: Yet another quadratic programming solver for robotics and beyond[END_REF].

min f,τ,a 1 2 ∥f -F ∥ 2 (27) subject to M a + g = J T f + S T τ + S T f friction Ja = -J q,
where J is the robot contact Jacobian, M is the mass matrix, S is the selection matrix that projects on the actuated joints, and g is the robot gravity vector at the current time step. The static friction in each joint was estimated independently and is approximately equal to 0.07. However, to keep a continuous model, we consider f friction = -0.07 2 π arctan (2S q). The two constraints ensure dynamics satisfaction and model that the end effectors do not move. We update our state estimate (including the external force torque) using the RS-EKF at 200 Hz with µ = 6. We also use the EKF in place of the RS-EKF to compare the performance of a risk-sensitive filter while keeping the update frequency the same. For both filter, we consider the following parameters:

P 0 = Q = BlockDiag(10 -3 I 6 , 10 -4 I 3 , 10 -1 I 3 , 10 -2 I 3 ) R = BlockDiag(10 -4 I 3 , 10 -2 I 3 , 10 -4 I 3 ).
The first experiment ran on the robot is shown in Fig. 5. Here, the base of Solo12 is pulled up (in the z direction) until an external estimate of 20N is computed by both the filters (Shown in bottom sub-figure 5 at time 1.4s with a vertical line). The base is then released to let Solo12 recover and bring its base back to the nominal desired height. The top sub-figure in Fig 5 shows an overlayed pictorial comparison of Solo12 when the two filters (RS-EKF and EKF) are used on the robot. Each frame corresponds to the state of Solo12 at the same time. The opaque Solo12 image shows the state with RS-EKF, and the blue one shows the EKF. This experiment shows that the RS-EKF helps the OCP to react quicker and bring the base to the nominal location sooner. This happens because the RS-EKF underestimates the base height in z as compared to EKF, which makes the OCP generate higher ground reaction forces to bring the base up sooner. In the end, both filter estimate similar external forces. As it can be seen in Fig. 5, the external vertical force does not converge exactly to zero. We find experimentally that this estimated force is due to friction. Lastly, as it can be seen, the cost of the RS-EKF-based controller is lower after the robot is dropped. The average cost of the RS-EKF is 0.065 while the one of the EKF is 0.130, which corresponds to a 50% improvement. This demonstrates how a filter aware of the cost objective can improve the overall performance.

In order to get a more systematic comparison between the filters and get rid of the human error, we perform two additional experiments where the filters are initialized with exactly the same priors. First, we initialize both the filters with a wrong prior on the external vertical force of 20 N, while in reality, no force is applied on the robot. This experiment creates an identical situation as the previous experiment while also ensuring the very same initial conditions for the robot. The results are shown in Fig 6, where the RS-EKF still performs better. Also, the performance is similar to the first experiment. In that experiment, we obtain a 62.9% improvement in the average cost. Finally, we replicate the previous experiment but now, initialize the filters with a wrong prior of -10 N on the external force, while, in fact, there is no force on the robot. The RS-EKF reacts sooner than EKF once again. It brings the base of Solo12 back to the desired location sooner than EKF, as can be seen in Fig 7. In that experiment, we obtain a 58.9% improvement in the average cost. V. CONCLUSION To conclude, we have introduced a risk-sensitive variation of the EKF based on the zero-sum game introduced by Whittle [START_REF] Whittle | Risk-sensitive linear/quadratic/gaussian control[END_REF]. The filter biases estimates towards high regions of the control cost which result in more robust controllers. Furthermore, the complexity of this filter is similar to the EKF. Lastly, we have demonstrated on several robotics problems, both in simulation and on real hardware, the benefits of this filter for output-feedback MPC. Importantly, we have shown that the proposed filter makes the controller more conservative in phases of high uncertainty leading to better overall control performance. Future work will address scenarios involving more dynamic tasks and more complex measurement uncertainty.

Fig. 1 .

 1 Fig. 1. Mass estimation for both EKF and RS-EKF.

Figure 2

 2 Figure 2 depicts the trajectory in space for both methods.It can be seen that the controller relying on the risk-sensitive estimate is more reactive. To evaluate both controllers, we compute the average Mean Square Error (MSE) relative to the reference trajectory. The MSE of RS-EKF is 0.0011, and the one of EKF is 0.0024. Hence, in that situation, the risk sensitivity brings a 54% improvement in tracking. Furthermore, the average cost along the trajectory is equal to 0.0569 for the RS-EKF-based controller, while it is equal to 0.0880 for the EKF-based controller, yielding a 35%

Fig. 2 .

 2 Fig. 2. Quadrotor trajectory for both the EKF-MPC and the RS-EKF-MPC.

Fig. 4 .

 4 Fig. 4. Median MSE on the tracking over 10, 000 experiments with random external disturbances. The envelope represents the 25th and 75th percentiles.

Fig. 5 .

 5 Fig. 5. Comparison of both methods after an external force of 20 N is applied by pulling the robot vertically. The vertical line indicates the moment when the robot is dropped. The top sub-figure overlays the trajectory of both the EKF in blue and RS-EKF in solid. The RS-EKF-based controller is more reactive to the perturbation and returns to the reference sooner.

Fig. 6 .

 6 Fig. 6. Comparison of the RS-EKF and EKF when initialized with a wrong prior of 20 N on the estimated vertical external force.

Fig. 7 .

 7 Fig. 7. Comparison of the RS-EKF and EKF when initialized with a wrong prior of -10 N on the estimated vertical external force.

APPENDIX

By taking a quadratic approximation of the value function, the Gauss-Newton step can be written as:

where

It can then be found that x t-1 = Q-1 q, where:

xt-1 by using the Woodbury lemma [START_REF] Thrun | Probabilistic robotics[END_REF], It can be shown that:

where Pt is defined as in (4). Finally, using

Pt , we can show that:

where P t , μt are defined as in ( 6), [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF].