Gaétan Lefebvre
email: g.lefebvre@brgm.fr

Mining economist opinions on using multi-agent methodology to simulate metal markets

Keywords: metal markets, agent-based models, prospective simulation

In coming years, mining economists expect the international demand for metals to rise. Market distortions are likely to appear, because producing countries controlling the resources can implement exportation quotas for certain categories of metals. Consequently, consuming countries like France face increasing risks of shortages of some metals. In response, the French government suggests that a set of prospective tools be implemented, including a prospective simulation tool based on a multi-agent system approach. The goal of this paper is to report and discuss critiques by mining economists on the value and abilities of multi-agent systems (MAS) to simulate critical metal markets whenever the approach is appropriate. The critiques are collected by confronting an existing agent-based computational economics (ACE) model of the lithium market with market reality and then by applying the discussions generally to other metals. The motivation is (a) to define the gaps currently existing between MAS and mining (b) to provide indicators on how the ACE model should be enhanced to reduce these gaps and (c) to produce an improved initial metal classification that is suitable for such a modelling exercise.

INTRODUCTION

In coming years, the growing demand for industrialization processes and progress in technologies means that mining economists will expect the international demand for metals to rise. At the same time, the world market for metals is becoming more and more complex, involving a great variety of players and mechanisms radically different from past decades. For instance, the concept of "critical metals" [START_REF] Graedel | Methodology of Metal Criticality Determination[END_REF] has arisen. Here two situations meet. On one hand, some consuming countries (such as France) and their manufacturing industries depend on importing a certain category of metals (antimony, rare earths, lithium, etc.) necessary for new technologies but unavailable in these countries. On the other hand, producing countries can destabilize the market thanks to their dominant position, by implementing exportation quotas for instance. Consequently, the consuming countries are at risk of shortage. To deal with these situations, the French government is suggesting the implementation of a set of prospective tools. This would allow the government to answer the following question: Given supply uncertainty, how long would a metal supply shortage last (should the case arise) in the world market and in France? One element of this set, which we suggest here, is a prospective simulation tool based on a multi-agent system (MAS) approach. More specifically, the tool contains agent-based computational economics (ACE) models i.e. models of artificial markets that are populated by economic interacting agents [START_REF] Tesfatsion | Agent-based computational economics: A constructive approach to economic theory[END_REF].

The goal of this paper is to report and discuss critiques by mining economists regarding an existing ACE model (see Section 2.2) that simulated the world lithium market. The motivation behind collecting these critiques is that these economists will be more willing to accept the ACE approach for the simulation of metal markets whenever the approach is appropriate. The critiques are collected by confronting the model (called "reference model" in this paper) with market reality and beyond the unique case of lithium. The reference model was implemented by MAS modelers with a view to long-term development of methodologies for metal market forecasts but mining economists need short-term operational tools and are not fully familiar with MAS. This paper does not present an alternative model (i.e. a 'better' tested and validated model than the reference model) stemming from these critiques, or plan to improve any agent architecture. This is a more theoretical discussion that aims (a) to define the gaps currently existing between MAS and mining economists (b) to indicate how the ACE model should be enhanced in the future to reduce these gaps and (c) in parallel, to produce an improved initial metal classification that is suitable for such a modelling exercise.

Section 2 of this paper presents the state of the art for metal market modelling and includes a reminder on the reference model that we discuss throughout the paper. Section 3 presents the critiques of the model from the mining economists as is. Section 4 presents the responses to these critiques by the MAS modelers and suggests possible changes in the model, based on these critiques. Section 5 concludes the paper.

STATE OF THE ART

2.1 Preamble: The general structure of metal markets

Classification of metals

On international markets, metals are usually divided into four categories: base metals, ferrous metals, precious metals and minor metals (Table 2 2.1 are the main references for the market to obtain the global reference of prices for base, ferrous and precious metals. By contrast, for the minor metals, no complete trading platform exists for their markets. In fact, the FANYA Metal Exchange, created in China in 2011, should gradually play this role in the future, but it cannot pretend to be a centralized market place yet. In most cases, trading for this category of metals is done directly between the seller (producer) and the buyer (final industrial manufacturer), without any supervisory authority (an "Over the Counter" exchange).

Metal markets

Metal markets should not be compared with other commodity markets. The logistics of the process of mining metallic deposits are very different from the logistics of exploiting commodities like oil and gas deposits, or of growing crops: the geological nature of metallic deposits (where one deposit can contain many metals) leads to a particular model of industrial development (the metal having higher price will be more exploited from this deposit) and thus, specific impacts on the availability of metal-based products on international markets1 . Generally, the following steps take place in the production process of metals: (a) extraction, where various metals can be obtained, some being by-products of others as a function of their prices (b) metallurgy treatment, i.e. all the processes from separating the mineral of interest to obtaining a tradable chemical product (oxide, ingot, powder), (c) transformation into pieces usable by the manufacturing industries, (d) distribution to these manufacturers and (e) distribution of the final metal-containing goods (e.g. cars having batteries containing lithium, lead, etc.) to consumers.

Summary on the MAS reference model

The ACE reference model [START_REF] Andriamasinoro | Proposal of agent simulation methodology for the prospective analysis of mineral commodities markets[END_REF] discussed in this paper simulated a lithium market. The product was lithium carbonate (Li 2 CO 3 ), an important element in the manufacture of end-use batteries. The model then excluded other products like lithium hydroxide (LiOH) or lithium chloride (LiCl).

Model design was motivated by the fact that whereas MAS generates increasing interest in the social or economic modeling of many types of market commodities (such as agricultural commodities [START_REF] Torii | Modeling agents and interactions in agricultural economics[END_REF], fishing commodities [START_REF] Soulié | Modeling fleet response in regulated fisheries: an agent-based approach[END_REF] or energy commodities [START_REF] Van Pruissen | Energy Efficiency Comparison of a Centralized and a Multi-agent Market Based Heating System in a Field Test[END_REF]) the world of geosciences has shown very little interest in MAS to model mining metal commodities. Geoscience, and especially models dealing with supply shortages, has adopted only the global scale as the level of their studies. Furthermore, to determine probable rupture periods at that level, production and consumption have been estimated in an independent manner and the results then compared arithmetically. The prospective studies are based on a simple extrapolation of current market trends. Consequently, there is no reciprocal balancing between supply and demand (absence of a systemic approach). Figure 2.1 illustrates this situation for lithium and rare earths respectively.

Figure 2.1: Prospective evolution of the lithium market till 2020 (left) [START_REF] Tru | Shocking Future Battering the Lithium Industry through 2020[END_REF] and the rare earth market till 2015 (right) [START_REF] Roskill | Rare Earths & Yttrium: Market Outlook to[END_REF] However, efforts to explore the criticality of metals should not consider only the global level, because organizational differences make a uniform analytical approach for all organizational levels (i.e., global, national and local) impractical [START_REF] Graedel | Methodology of Metal Criticality Determination[END_REF]. In addition, the risks of distribution may be underappreciated when discussing resources at a pure global level [START_REF] Kushnir | The time dimension and lithium resource constraints for electric vehicles[END_REF] because metal resources are distributed unevenly among countries. These arguments emphasize what [START_REF] Labys | New Directions in the Modeling and Forecasting of Commodity Markets[END_REF] already concluded a decade ago regarding modeling commodity markets (including metals) and models' roles in economic development: new frontiers of research should realize that commodity market behavior is intertwined with an international economic mechanism that includes globalization and expanded trade as well as interactions with developed and developing macroeconomics, including related financial institutions. This new vision would allow the government of a consuming country, especially in a restriction policy context, to better evaluate the impacts of the individual behavior and constraints from producers on supply shortage periods (if any) of this government's country.

For at least these reasons, an ACE model (the reference model in this paper, summarized below) was developed. Agent-based modeling was an adequate response because it does not rule out coordination failures, instability or crisis [START_REF] Lengnick | Agent-based macroeconomics: A baseline model[END_REF], a crisis, where in the metal market, an important supply shortage may be a future occurrence. The reference model is then, in the best of our knowledge, one of the pioneers in applying a MAS approach to supply shortage issues in metal markets.

These next subsection recalls the reference model. It will however only present the necessary elements to understand this paper. For a more detail about the formalism of that ACE model, we refer the reader to [START_REF] Andriamasinoro | Proposal of agent simulation methodology for the prospective analysis of mineral commodities markets[END_REF].

The data source

The reference model used international trade data from (GTIS, 2012) as data sources.

The GTIS data presents flows (quantity and price) between producing countries and transit countries as well as between transit countries and consuming countries. However, only the quantity parameter had been considered by the reference model.

The prospective period of the simulation began in 2013. The historic period was between 2005 and 2012, a period when the lithium data necessary for that work were available.

Structure of the reference model

A country was modelled as an agent, which is either a producer, a consumer, or a transit (i.e. a country connecting producers and consumers). For various reasons (ad-ministrative, geographical, etc.), a given consumer can be supplied by the same producer via several transits.

A country may be in the following context: Normal, Restriction, Compensation or Waiting. The model also integrated agents called ambassadors which are delegate agents that handle the flow exchanged between countries (one ambassador per country peer group).

Behaviour of the reference model

At the beginning of the simulation, each country is in a Normal context. A Normal context is the market context of a supply without restriction. A Restriction context is a context where a producer country imposes a quota restriction to a consuming country.

In a Normal context, the interaction between a consumer cck and a producer pci via a transit tcj, at each time step of the market simulation, occurs by following the four stages below, in which the first two points concern the demand stage and the last two points concern the supply stage: 1. Each consumer cck having a demand cck.σd asks its ambassador A(tcj↔cck) (for each j) to calculate the quantity d(tcj←cck) to ask for from pci (for each i), the supply of which will next transit via tcj. A demand over time was calculated via two steps: (a) interpolating the time series of the GTIS data related to the demands from cck to tcj between 2005 and 2013, in order to obtain a regression line which would describe and prolong that demand evolution, and (b) removing, from the resulting interpolation, the current available stock that A(tcj↔cck) already has. 2. When tcj has received the demands d(tcj←cck) (for each k) it transfers them to the ambassador A(pci↔tcj) (for each i). The ambassador then calculates, from these demands, the part d(pci←tcj) for which pci will have to respond. This part is here calculated as being a linear regression on the GTIS data between 2005 and 2013 corresponding to d(tcj←cck) (for all k) 3. When the demand arrives at pci, the latter, in response, calculates the total supply s(pci→tcj) it will provide all consumers. In a normal context supply equals demand, whereas in a restriction context, pci imposes a rate restriction pci.ρr (with 0≤ pci.ρr ≤1). In any case, the supply is then sent by pci (for each i) to tcj via the ambassador A(pci↔tcj). 4. When tcj has received the supplies from pci, it calculates and transfers to cck its part, via the ambassador A(tcj↔cck) (for each j, the sum of which is cck.σs). This part is calculated as being a linear regression on the GTIS data between 2005 and 2013 corresponding to s(pci→tcj) (for all i).

At the end of each time step, the available instantaneous stock of cck is the difference between the initial demand cck.σd and the final supply cck.σs.

A Compensation context is a context where a producer country attempts to partly make up the lack resulting from that restriction by another country. In a compensation context, a consuming country receiving a restriction from pci (let i=r) changes its context from Normal to Compensation and asks all pci (i≠r) to make up the lacking stock. Each pci that receives the message and accepts to make up, either immediately switches its context from Normal to Compensation or waits for a delay. In the latter case, it first switches its context from Normal to Waiting before switching from Waiting to Compensation, once this delay expires. This delay may be necessary for diverse reasons specific to pci: inability to immediately respond, speculation, etc.

Simulation of the reference model and results

The selected producing countries (pci) were Chile (cl), China (cn) and the United States (us). The model also added a (virtual) country called the rest of the world (rw). The selected consuming countries (cck) were France (fr), the subject of that study and, again, the rest of the world (rw). Finally, the transit countries (tcj) are Belgium (be), Germany (de), United Kingdom (uk), Italy (it) and the Netherlands (nl). The model also added, again, the rest of the world (rw). All the ambassadors were next naturally created to connect all these countries (rw included). To comply with the available data in GTIS, the simulation time step is 3 months. We note, for example, quarter 2 of year 2019 as 2/2019.

The pattern of the proposed (and currently fictitious) prospective scenario was the following: one assumes that as of 2014, Chile restricts its supply rate by cl.ρs points. Following this situation, China accepts assuring compensation at a rate of cn.ρp points, and does so immediately. The United States also accepts, with a rate of us.ρp points, but only as of 2016. The purpose of the simulation then consisted in varying the values of these rates to find the shortage end date in France and in the rest of the world.

Table 2.2 shows in detail the list of different scenario instances proposed in this paper. An instance is made of the scenario identifier (written in brackets), the value of the restriction from Chile and the value of compensation, respectively from China and USA. The value chosen in this table also allows a policy maker to analyze the sensitivity of the lithium market after a variation in important indicators (e.g. here, the diverse rates).

id -cl. s +cn. p +us. p (a) -0.15 +0.3 +0 (b) -0.4 +0.3 +0 (c) -0.15 +0.1 +0 (d) -0.4 +0.1 +0.5 (e) -0.4 +0.1 +0.1 (f) -0.1 +0.1 +0
Table 2.2: of All the Scenarios, an Instance Being Composed of a Restriction from Chile (cl) followed by a compensation from China (cn) and USA (us)

The left of Figure 2.2 next shows the shortage end dates obtained for the rest of the world in all scenarios. In this figure, the value of 13,000 (in t/quarter), in absolute values, approximately represents the average demand of lithium of rw (according to GTIS). It means for example that in scenario (f), at the peak time of a supply shortage period, there is still a minimal value of around (13,000-5,800) t/quarter of lithium (more than 50%) which is supplied to this consumer.

The right of Figure 2.2 is the "France equivalent" of the left side, with an average demand of around 350 t/quarter (according to GTIS). In this figure, France reaches a full shortage in all the scenarios where the Chile restriction is high (-0.4), i.e. (b), (d), (e) and (f) and with different durations. The reason of this full shortage is that the linear regressions made on the GTIS data result to a behavior where the model first handles the rest of the world (rw). When the stock is close to 0 again for rw, then France is automatically considered.

Summary of the MAS reference model

The MAS reference model is a purely bottom-up model i.e. it does not have any central mechanism to drive the interactions between producers and consumers. Furthermore, the maximum value of (year) 2023 mentioned in Figure 2.2 is not a time limit imposed by the model user, but an emergent result obtained during the simulation: it was observed that none of the periods of shortage exceeded this year. However, the model was designed, not as a substitution for classical statistical or mathematical approaches but rather as a complement (coupling) to them. Statistical approaches have been integrated to describe, not a global phenomenon, but the behavior of certain agents to which the approach is suitable. It particularly concerns the ambassadors. The linear regression tests were used with this objective. In the resulting system, certain agents behave in a linear manner (ambassadors and transit countries) while others (producing and consuming countries) adopt more complex and more discrete behaviors, which depend on their context.

CRITIQUES OF THE MAS MODEL FROM MINING ECONOMISTS

Relying on macro-and microeconomic reasoning and arguments, mining economists highlight various limits of the reference model when it has been confronted to "real world" situations. The corresponding critiques are categorized into four aspects: the data source, the structure, the behavior, and the simulation. They are detailed below and are followed by conclusions and indicators to improve the model performance. MAS responses to these critiques will next be provided in Section 4.

Regarding the data source

The data source GTIS (which originates from customs data) only contains flows between countries. It does not consider internal country consumption and production. In addition, there is a non-traceability issue. It means that if GTIS and more generally any customs data source like EUROSTATS, or USGSthe United States Geological Surveyare (relatively) accurate concerning the quantities of primary products ex-changed between countries, they do not trace the quantities of exchanged metals when they are already included in end-use products (e.g. quantities of lithium in traded cars). Non-traceability issues also means that for some metals (e.g. rare earths), there is a multiplicity of codes corresponding to different products containing the metal in its different forms.

Regarding the simulation

The time step chosen by the reference model is 3 months and the main agents are only countries. As the prospective exercise is on a long term (≃10 to 15 yearscf. Figure 2.2), these choices seem to be correct.

However, it should be kept in mind that when the prospective exercise is on short (≃2 to 5 years) to medium terms (≃5 to 10 years), these choices (3 month-time step, only country agents, etc.) are questionable. Indeed, at those scales, major industrial companies also have roles, particularly via the influences of their quarterly to bi-annual published results on daily and monthly evolutions in prices and exchanges. The example of the copper prices presented in Figure 3.1 indicates that their cycles are becoming tighter and tighter, requiring the simulation time step to be reduced accordingly. As metal market prospective studies may also cover short and medium-term scales, this aspect should be considered in future work when appropriate.

Regarding the structural aspects

The agents' structures lack properties which could influence greatly the outcomes. Indeed, the reasoning and equations only take quantities of "lithium carbonate" as the input variable while many other factors would need to be taken into account:  in the 'ambassadors': price, geopolitical relations between countries, etc.  in the '(producer) countries': costs of extraction and transformation (economic viability of the project), costs of transportation, etc.  in the '(consumer) countries': internal demands, evolution of technologies It is also noted that the model does not contain the secondary circuit at all. That would allow evaluation of the quantities of metals that can be obtained by recycling end-use products.

Regarding the behavioral aspects

On the demand side

The forecasting technique adopted by the demand model is based on extrapolation of historical data, resulting in the linear trends presented in Figure 2.2. Effectively, such an extrapolation technique was frequently used by the literature to analyze long-term metal market forecasts. We can cite, for instance, [START_REF] Roberts | Duration and characteristics of metal price cycles[END_REF] Again, in conformity with comments regarding simulation time step (cf. Section 3.2), extrapolation may no longer be accurate to analyze medium-term trends of commodity prices and stocks. Indeed, there has been an important change in market fundamentals over at least the last two decades. Finance mechanisms have progressively been applied to commodity markets and have taken on non-negligible importance in their evolutions, introducing new key players such as banks and investment funds. In short, first order fundamentals (meaning global supply and demand, exchange rates, etc.) have become insufficient to represent the market structures and dynamics. There are now a number of second or third order parameters to take into account (speculative buyer behavior, shadow banking, resource nationalism, etc.). In addition, demand trends are becoming harder to forecast by relying on past data because, for most minor metals, they depend on the expansion of new promising technologies where they find their end use. All these situations also apply for short-term forecasting, where nonlinearity of economic data has already been admitted for decades [START_REF] Agnon | Nonparametric, nonlinear, short-term forecasting: theory and evidence for nonlinearities in the commodity markets[END_REF].

On the supply side

On supply side, the reference model put the strong hypothesis that when demand increases, producers can always increase their production capacity to respond accordingly. Unfortunately, the mining industry does not always allow such a mechanism. One of the best examples is rare earths: when China diminished its export by 40% in 2010, none of other producers could effectively compensate such a shortage even though other resources had been identified. The reason is that the industrial treatment facilities were not operational at that time. Currently, there are (non-exhaustive) factors that do not allow the above hypothesis to be realistic: availability of the resources, prices, political factors and internal demand of producing countries, which we explain below.

Availability of the resources: the strong hypothesis above primarily implies that the resource of compensating producers can compete in quality and price with the restrict-ing producer's resources. This parameter is far from obvious and raises the question: if other producers were able to supply more lithium before, why did not they do so? Most of the time, technical constraints explain this. These remain high barriers to alternative producers, who cannot modify their production instantaneously.

Prices: if the reaction of Chile restricting its lithium exportations could be a brutal evolution of lithium market prices, then other producers would also be affected by such a change and could also react by reducing their production voluntarily, awaiting better prices. Indeed, extraction cost is one of the main parameters in determining the rate of production of a mine so that only high enough metal prices can be a credible incentive to increase production (and only if technically feasible).

Political factors: related to the previous argument, price distortion raises the question of the possibility of "organized" market failures, such as monopolistic or oligopolistic situations, arrangements between producing countries, and free-rider issues. In such cases, price distortion could be advantageous for some producers and they would not take any measure to resolve the shortage issue quickly, similar to what OPEC members have done in oil price crises, for instance. The reference model tried to take this aspect into account by introducing a waiting step, allowing for a differed response from some countries, but the determinant of the corresponding variable is currently exogenous and set arbitrarily. Finding the determinants of this variable would be interesting.

Internal demand of producing countries: as producing countries usually import part of their supply from restricting producers (e.g. a large part of China and the USA's lithium also comes from Chile), those countries would also be affected in case of a uniform restriction (i.e. a restriction applied to all connected countries). Their eventual extra-production would then be used first for internal needs, rather than for international markets. Besides, since the chosen customs data (GTIS) do not take into account internal consumption (cf. Section 3.1), the supply value is not accurate.

Conclusions from the economists

MAS seems to be a promising approach to analyze metal markets. However, according to the mining economists, it is far from convincing for operational use. The reference multi-agent model described in this paper still needs a lot of refining to be able to describe market situations that are close to reality. In particular, the following questions are posed:  How can the MAS approach reproduce high non-linear variations of metal market supply, demand and prices for the short or medium term, in an endogenous way?  How can MAS integrate the complex characteristics of some metal markets (e.g. those involved in the exploitation of the 17 groups of rare earths)  For which spatial and temporal scales is MAS modelling more appropriate? Can the approach handle multi-scale situations?

Nonetheless, through this work mining economists have provided a first category of metals that are more suitable for MAS modelling: base metals (copper, aluminum, iron, nickel, etc.), which are less critical than minor metals (antimony, indium, lithium, rare earths, etc.) (cf. Table 2.1 for these categories). The main reason for concluding this is data availability. In fact, in the metal market field, there is an inverse relation between data availability and its criticality. This is due at least to the following reasons: (a) the quantity of metals produced (a few or hundreds of tons for minor metals compared to millions or billions of tons for base metals) while some data sources ignore quantity below a threshold, (b) the way metals are traded as explained in Section 2.1.1 and (c) the speculation around their transformation into high value-added products

Regarding the non-traceability issue outlined in Section 3.1, to temporarily bypass this issue, the economists suggest that future ACE models should limit the number of enduses they manage while diversifying and cross-validating different sources of data (GTIS, EUROSTATS, USGS, Study Groups and Producer Associations, etc.). Limiting means avoiding (at least for the moment) simulating a primary product intended for too many industry end uses. This was what has been done in the reference model with lithium carbonate (mostly used for lithium-ion batteries) but the reasoning is valid for any metals. Obviously, this will not solve the non-traceability issue itself. This issue could be progressively resolved only if manufacturing producers would communicate more about the quantities of metals contained in their products. Meanwhile, this temporary solution will allow the progress in evaluating the MAS approach.

Economists' suggestions for future work

One suggestion from the economists to progress in evaluating the approach and the downstream modelling exercise is demonstrating MAS capacity via a showing-bydoing approach [START_REF] Hamill | Agent-based modelling: the next 15 years[END_REF]. This would implement an ACE model which can reproduce the price behavior of some metals during a period containing the last crisis.

The crisis scenario would be the following:  an expected increase in demand of technologies which are dependent on a specific minor metal (cell phones with tantalum, new technologies with rare earths, etc.),  "consumptionist euphoria" and a speculative bubble by the producers  industrials and financial players building strategic stockpiles, reducing the available raw materials and leading to crazy price increases  a sudden decrease in demand and consequently price.

Another suggestion for future models is to consider integrating other flow modeling techniques, including the modeling flow analysis (MFA) approach. This technique would be applied to analyze, for a given time/metal, how the input (instock, importation, production, etc.) evolves to output (consumption, outstock, exportation, etc.) for each country. This is important because in the current reference model, the flows circulating inside a country have been modeled in a 'simple' way: for example we cannot see (a) for a producer like USA, its importation, (b) for a consumer like France, its exportation, and (c) for transit like Belgium, a more precise analysis on which imported quantities are kept inside the country and which are exported to the final destination. MFA may be an interesting complementary approach to identify such lacks.

DISCUSSION

This section summarizes how the MAS modelers respond to the economic critiques from Section 3, with a view to taking them into consideration appropriately in future modeling work.

Justification for using the MAS approach

Before answering the critiques, it is worth recapping the interest of MAS especially for policy makers.

In fact, the debate on the methodological foundations of macroeconomic theory has gained new momentum during the recent worldwide economic crisis [START_REF] Lengnick | Agent-based macroeconomics: A baseline model[END_REF], a crisis where, in the metal market, an important supply shortage may be a fu-ture occurrence. Jean-Claude Trichet, the former President of the European Central Bank states the following [START_REF] Trichet | Reflections on the nature of monetary policy non-standard[END_REF]: "… As a policy-maker during the crisis, I found (…) in the face of the crisis, we felt abandoned by conventional tools (…). We need to deal better with heterogeneity across agents and the interaction among those heterogeneous agents (…) Such approaches are worthy of our attention…" That is why practitioners, governments and central banks are today very conscious that their understanding of the economic mechanisms has to be improved so that the linkages between various countries' (agents') systems at the international level can be demonstrated more explicitly, linkages which affect aggregate behavior [START_REF] Kirman | The Crisis in Economic Theory[END_REF]. This necessity for scale change also applies to the particular case of metal supply shortage risk according to most of the studies measuring metal criticality gathered by [START_REF] Speirs | Materials Availability, Working paper III: Comparison of material criticality studies -methodologies and results[END_REF]. In fact, those studies, although carried out at only aggregate level, are interesting for policy makers because (a) they are currently used operationally as a decision-aid model by industrialists, the European Union (EU), etc. and (b) they contain rich criticality determinants like geopolitical factors (governance, transparency, etc.), geological factors (substance availability, etc.), environmental factors, and economic factors (supply concentration, etc.). However, one limit of these studies is that they tried to assign a unique aggregated criticality score. Criticality differs from one country to another (depending for example on the geological reserves of the respective countries). The above authors then suggest future studies that would make the interaction between variables resulting in this unique score more explicit. Again, MAS is a solution. The idea is to design a 'country version' of these indicators inside an ACE model (cf. Section 4.3 for a first indicator).

In sum, we have here a research aspect (MAS) and an operational aspect (critical metals indicators at an aggregate level) that are currently attracting the attention of policy makers for analyzing metal criticality, especially as both aspects can be (theoretically at this stage) coupled inside one model.

Regarding the comments from the economists

Regarding simulation aspects (Section 3.2), comments related to the scale and complexity clearly demonstrate a lack of communication between the MAS field and the mining field concerning the foundation of MAS and its potential. MAS can simulate complex systems [START_REF] Wooldridge | An Introduction To Multiagent Systems[END_REF] and can handle multiple time scales such as (if we take examples of ACE models) a day scale [START_REF] Vriend | ACE models of endogenous interactions[END_REF], a month scale [START_REF] Lengnick | Agent-based macroeconomics: A baseline model[END_REF], a quarter scale (the reference model in this work), a year scale [START_REF] Cheng | An agent-based commodity trading simulation[END_REF], etc. MAS can also handle multiple space scales (if we look at the above work again): firm, household/country/world and can consider "financial agents" like banks [START_REF] Tesfatsion | Agent-based computational economics: A constructive approach to economic theory[END_REF].

Regarding behavioral aspects, the strong hypothesis (i.e. assuming that compensations of Chile's restriction by China and USA are always possible) noticed by the economists can be explained as follows: the current ACE work is one of pioneers in the mining field and is also a part of a methodological development in metal market forecasts. So (conceptually) assuming that when producer countries restrict, other countries compensate to maintain a market equilibrium is common sense. This could not simply be verified in the reference model because this model has missing parameters (cf. Section 3.2) and the corresponding missing behaviors. That said, these comments to avoid a strong hypothesis must be obviously integrated in future simulations to make the model realistic and usable at an operational level.

Regarding structural aspects and the missing (economic, geopolitical and environmental) parameters in the reference model, there are no particular comments from the MAS modelers given that they are also aware of these lacks. Geopolitical variables will be integrated into our future work by relying on work concerning metal criticalities as outlined in Section 4.1. In fact, that work not only suggested what geopolitical variables to add to countries (e.g. governance indicators) but also indicated the corresponding data sources (World Bank, Transparency International, International Monetary Fund, etc.).

Regarding data issues and particularly the non-traceability issue, the MAS field cannot solve this alone. At this time MAS will follow the recommendations from the economists (cf. Section 3.5): (a) simulating one specific product intended to one specific industry end-use while (b) diversifying and cross-validating different data sources (GTIS, OECD, EUROSTATS, USGS, etc.).

Regarding the economists' expectation of introducing MFA, the comments clearly highlight the wish to reinforce the coupling of classical and MAS approaches in the modeling of metal markets instead of sticking to a solely MAS approach. This means MAS should be considered as a complementary approach rather than a substitute for the classical approach. The reference model had already started to consider this idea (i.e. coupling MAS and statistical approaches -cf. Section 2.2.5). But for future model developments, it could be further enhanced by henceforth methodologically extending the coupling to the other approaches. Concretely, these approaches are (a) material flow analysis, (b) work on metal criticalities mentioned above (which integrate geopolitical variables), (c) spatialization processes via geological maps and country maps (which consider geological variables and transport distances respectively) and (d) metal life-cycle analysis (LCA). Metal LCA [START_REF] Nuss | Life Cycle Assessment of Metals: A Scientific Synthesis[END_REF] will be introduced in future ACE work to evaluate the environmental impact of metal exploitation as recommended by the EU (EU, 2011) but in this case, inside a MAS simulation. This coupling to LCA can be carried out by considering environmental properties (energy end-use, global warming potential, human health implications and ecosystem damage) in agents' properties and how they are manipulated in agents' behaviors (mining, transport, transformation, recycling, etc.). Throughout these actions, the MAS modelers are aware that it is incumbent on MAS simulation, for better acceptance of its usefulness, to justify its added value in enriching existing methodologies regarding metal market analysis, not the opposite.

How to create the showing-by-doing demonstration

This section explains how the MAS modelers plan to actualize the showing-by-doing demonstration requested by the mining economists, i.e. reproducing some metal price behavior covering periods of the last crisis (cf. Section 3.6).

The metal

Simulating rare earth markets and its 17 groups of elements, as expressed by the mining economists in Section 3.5, is too premature at this stage although more interesting (because critical) and conceptually possible via MAS approach. The chosen metal for that demonstration will be copper. Since this is a medium-term simulation, the exercise will be carried out at a month-time step (not the quarter-time step adopted by the reference model) and with the idea of reproducing high price variations (cf. Figure 3.1) and demand in an endogenous way.

Copper is chosen for the following reasons and more:  Copper's economic data (price, quantity) are more available than data on lithium or rare earths, not only for import-export flow but also for production and consummation inside a country (see Section 3.5 for the explanation). This is true for all primary or secondary (recycling) circuits of this metal (via GTIS, EUROSTAT, IMF, etc.) but also for environmental evaluation (the Eco-invent database).  Copper end-use is, at more than 50% of market share, intended for the construction and automotive sectors where data is more available than for some other metal end-uses. For a modelling exercise, these data (sharing more than 50% of enduse market) are not negligible so they diminish the effect of the non-traceability issue.  We could refer to recent studies on the criticality of copper (i.e. including geopolitical aspects) [START_REF] Nassar | Criticality of the Geological Copper Family[END_REF].

Baseline of the future ACE model

Regarding the future copper ACE model, the following baseline indicators are suggested.

The design of links between producers and consumers would be modeled by direct sequences of ask-bid processes between producers and consumers according to their local knowledge of the market. The future work would then keep the fully bottom-up approach adopted by the reference model (cf. Section 2.2.5) while avoiding the extrapolation process limiting that model as argued in Section 3.4.1. The future work would also rely on the work of [START_REF] Cheng | An agent-based commodity trading simulation[END_REF]. That is interesting because (a) it considers additional market agents (hedgers, traders, etc.), (b) during an ongoing simulation, it considers different market external events that drive complicated price movement with only simple agent strategies and (c) it has been applied to the simulation of world crude oil price evolution, theoretically facilitating its transposition to our copper trading market. Regarding the evaluation of the demand in that future work (mostly in the construction and automotive sectors, as mentioned previously), it would rely on the ACE work of [START_REF] Lengnick | Agent-based macroeconomics: A baseline model[END_REF]. That work is a fully bottom-up market model that computed agents' demands according the household needs (i.e. a variable from which construction and automotive needs could also be computed).

Regarding the future of MAS acceptance

It seems to be clear that MAS is a promising tool for prospective simulation of metal markets, among other reasons because it would be of interest for a policy maker (Section 4.1) in market analysis after the recent economic crisis. MAS is also open enough to be theoretically coupled with more classical existing approaches (Section 4.2). However, while the intention behind these MAS responses is to conform as close as possible to mining economists' suggestion, success will be defined by the economists accepting MAS and is far from obvious for at least two reasons: data obstacles and paradigm obstacles.

A data obstacle

As demonstrated in the entire paper, the data available in metal markets are not 'ready' yet for MAS simulation. One reason is because the industrial metals market at an international scale has never considered the necessity of constructing and organizing data for a simulation purpose. In addition, at more detailed levels, data are confidential to varying degrees, published according to criteria variable from one country to another (the metal, the producing company, the country, etc.). Consequently, although rare earths would be the next metal whose market is interesting to simulate for the showing-by-doing MAS demonstration suggested by the economists (i.e. the capability for simulating the last crisis), copper was finally chosen not because it is the most critical metal but because its market has one of the best quantities of available data. It is hoped that once a future copper model is implemented and validated, the resulting work would incite the mining field to progressive collection and validation of raw material data at the international level for a simulation purpose. The mediumterm idea is to concretely succeed in creating a cyclic process between (a) increased successful transpositions of the reference model, (b) better organization and reliability of data (c) better acceptance of the underlying method, (d) progression in the demand of MAS approach for studies from different entities (consultants, scientists, etc.) at an operational level. Such a cyclic process already exists and is recognized at an operational level but only for a global level, not at an international level.

A paradigm obstacle

Even though it is assumed the showing-by-doing test expected by economists would produce interesting results, we detected many sources of reluctance in using the approach. First, we recognize that from the economist's viewpoint, agent-based modelling has an obvious drawback: it makes it impossible to think in aggregate terms. Multi-agent modelling also has a destructive consequence: the dimension of the model explodes [START_REF] Assenza | E Pluribus Unum: Macroeconomic modelling for multi-agent economies[END_REF]. Second, a MAS is composed of autonomous agents and emergent phenomena while mining economists would like to "keep control" of their analysis, i.e. not letting computers reason in lieu of them. Third, how long it will take to design one application scenario is not clear and during a simulation, the time each agent will take (in real time) to behave is not clear either. There is an uncertainty on how many times, at the most, such potential users should wait before obtaining one scenario result (at a scale of a minute? a week?). All of these sources of reluctance are unlikely to disappear in the short term but we should at least constantly highlight them in future MAS research papers concerning progress towards better acceptance of the MAS method in the metals field.

CONCLUSIONS

In coming years, mining economists expect the international demand for metals to rise. Market distortions are likely to appear because producing countries controlling the resources can implement exportation quotas for certain categories of metals, for political or environmental reasons. Consequently, consuming countries like France are facing increasing risks of shortages for some metals. In response to that, the French government suggests the implementation of a set of prospective tools among which is a prospective simulation tool based on a multi-agent system approach. The goal of this paper was to report and discuss critiques carried out by the economists regarding the value of and ability of a multi-agent system to simulate critical metal markets whenever the approach is appropriate. The critiques were collected by comparing an existing agent-based computational economics (ACE) model of the lithium market to market reality and by generalizing the discussions to other metals. The motivation is (a) to define the gaps currently existing between MAS and mining economists (b) to provide indicators on how the ACE model should be enhanced in the future to reduce these gaps and (c) in parallel, to produce an improved initial metal classification that is suitable for such a modelling exercise. We can summarize the results as follows.

First, finance mechanisms have progressively taken on non-negligible importance in the evolution of macroeconomics (and even more so metal markets, a field of macroeconomics), introducing new key players such as banks and investment funds. Practitioners, governments and central banks are today very conscious that knowledge of the economic mechanisms has to be improved so that the network structure of the financial industry can be more explicitly demonstrated at an international level, a structure which next has an important impact on aggregate behavior. Consequently, unlike what is proposed by the reference model, prospective modeling of metal market demand by only extrapolating its historical data will be less and less appropriate. Trends are becoming harder to forecast by relying on past data because, for most of the minor metals, they depend on the expansion of new promising technologies in which they are used, which could eventually be replaced by even newer technologies in a decadelong period. In sum, demand and price variation will be more and more non-linear and the prospective modeling technique should be improved to better identify the determinants of that evolution.

Second, a Multi-Agent System approach and more specifically, for this work, Agent-Based Computational Economics (ACE) models, is conceptually recognized as appropriate for simulating the complex metal market as described above. The recognition comes from scientists but also from policy makers, especially after the recent economic crisis. Nonetheless, we also recognize that there is a lack of communication to mining economists regarding the MAS foundation and potential for simulating such a market. Thus, for progress towards a better acceptance of MAS by the economists: -mining economists expect demonstration by a showing-doing-approach, in particular an ACE model which would be capable of reproducing high price and demand variation of for a metal during a period including the recent crisis, in an endogenous way.

Copper will be chosen for that demonstration.

classical modeling approaches should be more and more highlighted in the model as far as appropriate, instead of a solely MAS structure.

Third, available data in metal markets are not 'ready' yet for MAS simulation. Consequently, although rare earths would be the next metal whose market was interesting to simulate, copper has been chosen not because it is the most critical metal but because its market has one of the best quantities of available data. More generally, this work has allowed mining economists to provide a first category of metals that are more suitable for a MAS modelling: base metals (copper, aluminum, iron, nickel, etc.) rather than minor metals (antimony, indium, lithium, rare earths, etc.). It is hoped that once a future copper model is implemented and validated, the resulting work would incite the mining field to progressively collect and validate more critical metal data at the international level for a simulation purpose.

Finally, we recognize that even though it is assumed that the showing-by-doing test expected by mining economists would produce interesting results, they will still be reluctant to adopt the MAS approach as a tool at the operational level, at least in the short term, among other reasons because (a) they are used to keeping control of their analysis, while MAS provides autonomous agents and emergent mechanisms and (b) MAS currently contains a level of complexity that is not still suitable for them in terms of model and simulation designs. These sources of reluctance will not immediately disappear but we should constantly highlight them in future MAS research papers to aid progress towards better acceptance of the MAS method in any application.

In terms of perspectives for the present work, the next step is to achieve the showingby-doing demonstration. We plan to base that work on ACE models already existing in the literature [START_REF] Cheng | An agent-based commodity trading simulation[END_REF][START_REF] Lengnick | Agent-based macroeconomics: A baseline model[END_REF].

Figure 2

 2 Figure 2.2: Prospective Evolution of the Lithium Stock in the Rest of the World (left) and in France (right), for all Scenarios

Figure 3

 3 Figure 3.1: Variation of the copper price between 2002 and 2015 (BRGM, 2015)

  who analyzed the prices of 14 metals from January 1947 through December 2007, and (Shafiee & Topal, 2010) who analyzed the prices of gold from January 1968 to December 2008.

  Figure 4.1 shows an example of the efficiency of such a model, in particular the ability of one agent trader to progressively comply with the trend indicated by real-world market price.

Figure 4 . 1 :

 41 Figure 4.1: Example of a trader agent behavior (Cheng et al., 2009): thick line is trader's position balance; thin line is real market price; dotted line is estimated market price dynamics

Table 2

 2 .1).

	Category	Chemical	Trading platform
	Base metals	Aluminum, copper, lead, zinc, tin, nickel	London Metal Ex-
	Ferrous metals Iron ore, chromium, manganese, molyb-	change (LME) and
		denum, niobium, vanadium	Shanghai Futures Ex-
	Precious metals Gold, silver, platinum, palladium	change (SHFE)
	Minor metals	Antimony, beryllium, bismuth, cobalt, gal-lium, germanium, hafnium, indium, lithi-um, rare earths, rhenium, selenium, tanta-lum, tellurium, tungsten, zirconium	FANYA Metal Ex-change (at a starting stage)

.1: Metal categories and trading platforms LME and SHFE introduced in Table

  Table 3.1 gives an example of the issue in the data source (EUROSTAT, 2014) for rare earths. In this table, there are similarities between the products having different codes, making traceability difficult.

	HS code	Description in EUROSTAT
	28053010 Intermixtures or interalloys of rare-earth metals, scandium and yttrium
	28053090 Rare-earth metals, scandium and yttrium (Excl. intermixtures of interal-
		loys)
	28461000 Cerium compounds

28469000 Compounds, inorganic or organic, of rare-earth metals, of yttrium or of scandium or of mixtures of these metals (Excl. cerium)

Table 3 .

 3

1: HS codes for raw rare earth product imports and exports

(EUROSTAT, 2014)

Such heterogeneity exists less in other commodities.


Felipe Castillo

Emmanuel Witrant

Christophe Prieur

Vincent Talon

Luc Dugard

Fresh Air Fraction Control in Engines Using Dynamic Boundary Stabilization of LPV Hyperbolic Systems

In this paper, we consider the boundary control of the fresh air mass fraction in a Diesel engine operated with low-pressure exhaust gas recirculation. The air mass fraction transport phenomenon is modeled using a cascade of first-order linear parameter varying (LPV) hyperbolic systems with dynamics associated with their boundary conditions. By means of Lyapunov based techniques, sufficient conditions are derived to guarantee the exponential stability of this class of infinite dimensional systems. We develop a polytopic approach to synthesize a robust boundary control that guarantees the exponential stability for a given convex parameter set. Simulation results illustrate the effectiveness of the proposed boundary control to regulate the mass fraction of fresh air in a Diesel engine.

I. INTRODUCTION

Regulations of Diesel engine emissions have become stricter, and satisfying simultaneously the emissions legislations and the desired engine drivability objectives is a particularly challenging issue. Although significant improvements were made over the past years, there are still many technical issues that need to be addressed in order to meet the future regulation laws on emissions. The introduction of sophisticated alternative combustion modes such as homogeneous Felipe Castillo, Emmanuel Witrant, Christophe Prieur and Luc Dugard are with UJF-Grenoble 1/CNRS, Grenoble Image Parole Signal Automatique (GIPSA-lab), UMR 5216, B.P. 46, F-38402 St Martin d'Hères, France. (felipe.castillo-buenaventura, emmanuel.witrant, christophe.prieur, luc.dugard)@gipsa-lab.fr. This work is partially supported by HYCON2 Network of Excellence (Highly-Complex and Networked Control Systems), grant agreement 257462. Felipe Castillo and Vincent Talon are with Renault SAS, 1 Allee Cornuel, 91510 Lardy, France (e-mail: vincent.talon@renault.com) charge compression ignition (HCCI), low temperature combustion (LTC) and premixed controlled compression ignition (PCCI) offers a great potential to reduce the engine emissions levels [START_REF] Akihama | Mechanism of the smokeless rich Diesel combustion by reducing temperature[END_REF] [2] [START_REF] Ryan | Fuel requirements for HCCI engine operation[END_REF]. However, these new modes require specific fueling strategies and in-cylinder conditions, thus creating the need for more complex, reliable and precise control systems and technologies.

Dual-loop exhaust gas recirculation (EGR) with both high (HP) and low-pressure (LP) recirculations is one of the new strategies that can provide the appropriate conditions for multiple combustion modes [START_REF] Hribernik | The potential of the high and low-pressure exhaust gas recirculation[END_REF]. Ensuring the adequate in-cylinder conditions is still a difficult task, due to the lack of measurements for EGR flow rates and mass fraction. Several air mass fraction/EGR rate control and estimation methods have been proposed in the literature [START_REF] Chauvin | Motion planning for experimental airpath control of a Diesel homogeneous charge-compression ignition engine[END_REF] [START_REF] Grondin | Control of a turbocharged Diesel engine fitted with high pressure and low pressure exhaust gas recirculation systems[END_REF] [34] [START_REF] Yan | Control of a dual loop EGR air-path systems for advanced combustion Diesel engines by singular perturbation methodology[END_REF].

However, most of the actual air fraction control techniques are based on zero-dimensional (0D) engine modeling, which does not allow taking into account the air transport inside the engine admission air-path. Indeed, the mass transport causes a degradation of the overall engine emission performance during strong transient conditions. This is mostly due to the LP-EGR, as the distance that the gas travels in the engine air-path is much longer than the one associated with HP-EGR.

The control of air fraction taking into account the mass transport time is a significantly less explored subject, although in [START_REF] Bresch-Pietri | Prediction-based trajectory tracking of external gas recirculation for turbocharged SI engines[END_REF] a non-physical representation of the delay by means of PDES was used to design a predictor-based controller.

The flow transport strongly depends on the engine operating conditions. It has been shown in [START_REF] Castillo | Simultaneous air fraction and low-pressure EGR mass flow rate estimation for Diesel engines[END_REF] and [START_REF] Wang | Air fraction estimation for multiple combustion mode Diesel engines with dualloop EGR systems[END_REF] that this dependence can be modeled by using a linear parameter varying (LPV) approach. On the other hand, the transport of mass is often modeled by first order hyperbolic partial differential equations (PDEs), as reported in [START_REF] Bastin | Using hyperbolic systems of balance laws for modelling, control and stability analysis of physical networks[END_REF], [START_REF] Castillo | Dynamic boundary stabilization of linear and quasi-linear hyperbolic systems[END_REF] and [START_REF] Santos | Boundary control of open channels with numerical and experimental validations[END_REF]. An LPV hyperbolic model can thus capture the dynamics of the air mass fraction and be used to solve the associated control problem.

Several results are available in the literature for the control of first order hyperbolic systems.

For instance, sufficient conditions for controllability and observability of quasilinear hyperbolic systems have been obtained in [START_REF] Li | Controllability and observability for quasilinear hyperbolic systems[END_REF]. The boundary control of hyperbolic systems has been considered in [START_REF] Coron | A strict Lyapunov function for boundary control of hyperbolic systems of conservation laws[END_REF], [START_REF] Coron | Dissipative boundary conditions for onedimensional nonlinear hyperbolic systems[END_REF], [START_REF] Di Meglio | Backstepping stabilization of an underactuated 3x3 linear hyperbolic system of fluid flow transport equations[END_REF] and [START_REF] Prieur | Stability of switched linear hyperbolic systems by Lyapunov techniques[END_REF], among other references. Most results consider that the boundary control can react fast enough in comparison with the waves travel time. More precisely, no time response limitation is taken into account at the boundary conditions. For many applications (e.g. [5] [18]), the wave travel can be considered as being much slower than the actuator time response. A static relationship can then be established between the control input and the boundary condition. Nevertheless, there are applications where the dynamics associated with the boundary control cannot be neglected (e.g. the use of LP-EGR to control the fresh air fraction in an engine intake manifold). To address this problem, a discretization of the infinitedimensional system has been used in [START_REF] Castillo | Contrôle de température dans un flux de Poiseuille[END_REF] to apply finite-dimensional control tools. In [START_REF] Castillo | Dynamic boundary stabilization of linear and quasi-linear hyperbolic systems[END_REF] some sufficient conditions for the exponential stability of hyperbolic systems with linear time invariant (LTI) dynamic boundary conditions have been given. In [START_REF] Prieur | ISS-Lyapunov functions for time-varying hyperbolic systems of balance laws[END_REF] and [START_REF] Hante | Converse Lyapunov theorems for switched systems in Banach and Hilbert spaces[END_REF], a strict time-varying Lyapunov function that allows establishing the asymptotic stability of time-varying hyperbolic systems is formulated for time-invariant boundary conditions. In [START_REF] Santos | A multi-model approach to Saint-Venant equations: A stability study by LMIs[END_REF], a multi-model approach (similar to an LPV) with a bilinear matrix inequality has been proposed for the stability of a hyperbolic system representing the flow in an open channel. In [START_REF] Krstic | Backstepping boundary control for first-order hyperbolic PDEs and application to systems with actuator and sensor delays[END_REF], backstepping designs for the boundary control of hyperbolic and time delayed systems have been proposed. However, to the best of our knowledge, no boundary control design for LPV hyperbolic systems with LPV boundary conditions has been addressed in the literature. This article, which is an extension of our previous work [START_REF] Castillo | Dynamic boundary stabilization of linear parameter varying hyperbolic systems: Application to a Poiseuille flow[END_REF], focuses on the stabilization of LPV hyperbolic systems with boundary conditions defined by LPV dynamics and its application to the control of the air fraction in the intake manifold of Diesel engines. To demonstrate the asymptotic stability for this class of hyperbolic systems, we use an extension of the strict Lyapunov function formulation presented in [START_REF] Coron | A strict Lyapunov function for boundary control of hyperbolic systems of conservation laws[END_REF]. A polytopic approach is developed to synthesize a robust boundary control that guarantees the exponential stability for a given convex parameter set. We obtain sufficient conditions in terms of the boundary conditions to prove Lyapunov stability. This feedback design strategy is used for the boundary control of the fresh air mass fraction in an engine air-path with LP-EGR. The effectiveness of our approach is evaluated in simulation using a 1D model of the engine admission air-path. This paper is organized as follows: the problem formulation and the main stability results for LPV hyperbolic systems with LPV dynamic behavior at the boundary conditions are presented in Section II. In Section III, an LPV hyperbolic system is used to model the air fraction in the engine air-path and an air fraction boundary control is designed using the results of Section II. In Section IV, simulation results are presented to illustrate the effectiveness of the proposed control strategy.

Notation. By the expressions H 0, H 0, H 0 and H ≺ 0 we mean that the matrix H is positive semi-definite, negative semi-definite, positive definite and negative definite, respectively.

The usual Euclidian norm in R n is denoted by |.| and the associated matrix norm is denoted . . Given g : [0, 1] → R n , we define its L 2 -norm (when is finite) as:

g L 2 = 1 0 |g(x)| 2 dx

II. STABILITY OF PARAMETER-VARYING LINEAR HYPERBOLIC SYSTEMS WITH DYNAMIC BOUNDARY CONDITIONS

Let n be a positive integer, Θ an open non-empty convex set of R n and Z ϕ a non empty convex set of R l . Consider the general class of first order LPV hyperbolic systems of order n defined as follows:

∂ t ξ(x, t) + Λ(ϕ)∂ x ξ(x, t) = 0 ∀ x ∈ [0, 1], t ≥ 0 (1)
where ξ : [0, 1] × [0, +∞) → Θ, ϕ is a varying parameter vector that takes values in the parameter space Z ϕ , Λ(ϕ) : Z ϕ → R n×n is a diagonal and invertible matrix function (called the characteristic matrix) such that Λ(ϕ) = diag(λ 1 (ϕ), λ 2 (ϕ), ..., λ n (ϕ)), ∂ t and ∂ x denote the partial derivatives with respect to time and space, respectively.

Assumption 1: Assume that the following inequalities hold for all ϕ ∈ Z ϕ :

0 < λ 1 (ϕ) < ... < λ n (ϕ) (2)
Consider the following dynamic boundary conditions for (1):

Ẋc = A c (ϕ)X c + B c (ϕ)u (3) ξ(0, t) = C c X c + D c u with u = Kξ(1, t) (4)
where

X c ∈ R nx , A c : Z ϕ → R nx×nx , B c : Z ϕ → R nx×n , C c ∈ R n×nx , D c ∈ R n×n , K ∈ R n×n ,
u ∈ R n and n x ≥ 1. Define the initial condition for (1) and (3) as:

ξ(x, 0) = ξ 0 (x), ∀ x ∈ [0, 1] X c (0) = X 0 c ( 5
)
where ξ 0 (x) ∈ L 2 ((0, 1); R n ) and X 0 c ∈ R nx . Such an initial condition is required for the existence of a unique classical solution of the Cauchy problem (1) -( 5) [START_REF] Diagne | Lyapunov exponential stability of 1-D linear hyperbolic systems of balance laws[END_REF]. In this section, we consider the problem of finding a control gain K and a positive scalar µ such that the system (1) with boundary conditions (3) and initial condition ( 5) is exponentially stable for all ϕ ∈ Z ϕ satisfying the following inequality.

||X c (t)|| 2 + ||ξ(t)|| L 2 (0,1) ≤ b ||X 0 c || 2 + ||ξ 0 || L 2 (0,1) e -µ t (6)
where b and are positive scalars.

Let the polytope Z ϕ be defined as follows:

Z ϕ := {[ϕ 1 , ..., ϕ l ] T ∈ R l | ϕ i ∈ [ϕ i , ϕ i ], ∀ i = 1, ..., l} (7)
for given l ∈ N + and the parameter extremities ϕ i , ϕ i (minimum and maximum, respectively).

We thus consider that all the admissible values of the vector ϕ are constrained in a hyperrectangle in the parameter space Z ϕ . Consider the polytopic linear representation of the parameter varying characteristic matrix for all ϕ ∈ Z ϕ [START_REF] Angelis | System Analysis, Modelling and Control with Polytopic Linear Models[END_REF]:

Λ(ϕ) = Nϕ i=1 α i (ϕ)Λ(w i ) (8)
where w i ∈ Z ϕ are the N ϕ = 2 l vertices of the polytope formed by all extremities (ϕ i and ϕ i ) of each varying parameter ϕ ∈ Z ϕ , 2 l i=1 α i (ϕ)Λ(w i ) : Z ϕ → R n×n and α i (ϕ) is a scheduling function α i : Z ϕ → [0, 1] defined as:

α i (ϕ) = l k=1 |ϕ k -C(w i ) k | l k=1 |ϕ k -ϕ k | (9)
where C(w i ) k is the kth component of the vector C(w i ) defined as:

C(w i ) k = ϕ k , if (w i ) k = ϕ k ϕ k , otherwise (10)
The scheduling functions α i have the following properties [START_REF] Angelis | System Analysis, Modelling and Control with Polytopic Linear Models[END_REF]:

α i (ϕ) ≥ 0, Nϕ i=1 α i (ϕ) = 1 (11)
The polytopic representation (8) can also be considered for the matrices A c (ϕ) and B c (ϕ) of the boundary conditions (3) using the same scheduling function [START_REF] Castillo | Intrapipe restriction non-homentropic boundary resolution method[END_REF]. Define the polytopic linear representation (PLR) of the parameter varying hyperbolic system (1) with boundary conditions

(3) as follows:

∂ t ξ(x, t) + Nϕ i=1 α i (ϕ)Λ(w i )∂ x ξ(x, t) = 0 ∀ ϕ ∈ Z ϕ , ∀ x ∈ [0, 1], t ≥ 0 (12)
with boundary conditions

Ẋc = Nϕ i=1 α i (ϕ)A c (w i )X c + Nϕ i=1 α i (ϕ)B c (w i )u ξ(0, t) = C c X c + D c u (13)
Based on the PLR ( 12) -( 13), the following theorem states a sufficient condition to ensure the exponential stability for system (1) with boundary conditions (3) and initial condition (5) for all ϕ ∈ Z ϕ .

Theorem 1. [Stability analysis] Along with Assumption 1, assume that there exists two diagonal positive definite matrices P 1 ∈ R nx×nx and P 2 ∈ R n×n and a scalar µ > 0 such that the following matrix inequality is satisfied, for all i = 1, ..., N ϕ ,

M i =      A c (w i ) T P 1 + P 1 A c (w i ) + C T c Λ(w i )P 2 C c P 1 B c (w i )K + C T c Λ(w i )P 2 D c K +µΛ(w i )P 1 K T B c (w i ) T P 1 + K T D T c Λ(w i )P 2 C c K T D T c Λ(w i )P 2 D c K -e -µ Λ(w i )P 2      0 ( 14
)
Then there exist two constant scalars a > 0 and b > 0 such that, for all ξ 0 ∈ L 2 ((0, 1); R n ) and X 0 c ∈ R nx , the solution of (1), ( 3) and (5) satisfies, for all t ≥ 0

||X c (t)|| 2 + ||ξ(t)|| L 2 (0,1) ≤ be -at ||X 0 c || 2 + ||ξ 0 || L 2 (0,1) (15)
Proof: Given the diagonal positive definite matrices P 1 and P 2 , consider (as an extension of the Lyapunov function proposed in [START_REF] Coron | A strict Lyapunov function for boundary control of hyperbolic systems of conservation laws[END_REF]) the quadratic strict Lyapunov function candidate defined for all continuously differentiable functions ξ : [0, 1] → Θ as:

V (ξ, X c ) = X T c P 1 X c + 1 0 ξ T P 2 ξ e -µx dx (16)
where µ is a positive scalar. Computing the time derivative V of V along the classical C 1solutions of (1) with boundary conditions (3) and initial condition (5), yields to the following:

V = ẊT c P 1 X c + X T c P 1 Ẋc + 1 0 ξT P 2 ξ + ξ T P 2 ξ e -µx dx (17)
After integration by parts and taking into account the polytopic linear representation ( 12) and ( 13), the following is obtained:

V = 2 l i=1 α i (ϕ) X T c A c (w i ) T P 1 + P 1 A c (w i ) X c + ξ(1) T K T B c (w i ) T P 1 X c + X T c P 1 B c (w i )Kξ(1) -e -µx ξ T Λ(w i )P 2 ξ 1 0 -µ 1 0 ξ T Λ(w i )P 2 ξ e -µx dx (18)
where ξ(1) = ξ(1, t). The previous equation can be written using the boundary conditions [START_REF] Angelis | System Analysis, Modelling and Control with Polytopic Linear Models[END_REF] as follows:

V = 2 l i=1 α i (ϕ) X T c A c (w i ) T P 1 + P 1 A c (w i ) X c + ξ(1) T K T B c (w i ) T P 1 X c + X T c P 1 B c (w i )Kξ(1) -e -µ ξ(1) T Λ(w i )P 2 ξ(1) + X T c C T c Λ(w i )P 2 C c X c + X T c C T c Λ(w i )P 2 D c Kξ(1) + ξ(1) T K T D T c Λ(w i )P 2 C c X c + ξ(1) T K T D T c Λ(w i )P 2 D c Kξ(1) -µ 1 0 ξ T Λ(w i )P 2 ξ e -µx dx = 2 l i=1 α i (ϕ) -µX T c Λ(w i )P 1 X c -µ 1 0 ξ T Λ(w i )P 2 ξ e -µx dx +   X c ξ(1)   T M i   X c ξ(1)   (19)
where the matrix M i is defined as in [START_REF] Chauvin | Motion planning for experimental airpath control of a Diesel homogeneous charge-compression ignition engine[END_REF]. The definition α i ≥ 0 and the matrix inequality 14) imply that the last term of ( 19) is always negative or zero. This gives the following inequality:

M i 0 from (
V ≤ 2 l i=1 α i (ϕ) -µX T c Λ(w i )P 1 X c -µ 1 0 ξ T Λ(w i )P 2 ξ e -µx dx (20)
From (2) it can be proved that there always exists a > 0 such that Λ(ϕ) -I n×n 0 (e.g. could be the smallest eigenvalue of Λ(ϕ) over Z ϕ ). Moreover, the diagonality of P 1 , P 2 and Λ implies that:

V ≤ -µ V (ξ, X c ) (21)
Therefore, the function ( 16) is a Lyapunov function for the hyperbolic system (1) with boundary conditions [START_REF] Angelis | System Analysis, Modelling and Control with Polytopic Linear Models[END_REF].

By integrating the inequality (21) from 0 to t, we obtain the following:

V (t) ≤ V (0)e -µ t (22)
To obtain the final result [START_REF] Chung | Computational Fluid Dynamics[END_REF], we bound the Lyapunov function as follows:

min{λ min (P 1 ), λ min (P 2 )}(||X c (t)|| 2 + ||ξ(t)|| L 2 (0,1) ) ≤ V (t) ≤ max{λ max (P 1 ), λ max (P 2 )}(||X c (t)|| 2 + ||ξ(t)|| L 2 (0,1) ) (23)
where λ min and λ max are the minimum and maximum eigenvalues of the considered matrices, respectively. Then, using [START_REF] Hribernik | The potential of the high and low-pressure exhaust gas recirculation[END_REF] together with [START_REF] Krstic | Backstepping boundary control for first-order hyperbolic PDEs and application to systems with actuator and sensor delays[END_REF] gives:

||X c (t)|| 2 + ||ξ(t)|| L 2 (0,1) ≤ max{λ max (P 1 ), λ max (P 2 )} min{λ min (P 1 ), λ min (P 2 )} ||X 0 c || 2 + ||ξ 0 || L 2 (0,1) e -µ t (24)
which implies that a = µ and b = max{λ max (P 1 ), λ max (P 2 )}/ min{λ min (P 1 ), λ min (P 2 )} in [START_REF] Chung | Computational Fluid Dynamics[END_REF].

Note that the matrix inequality ( 14) considers, through the Lyapunov matrices P 1 and P 2 , the dynamic coupling between the system and its boundary conditions. Inequality [START_REF] Chauvin | Motion planning for experimental airpath control of a Diesel homogeneous charge-compression ignition engine[END_REF] along with [START_REF] Hante | Converse Lyapunov theorems for switched systems in Banach and Hilbert spaces[END_REF] implies that µ is a tuning parameter of the controller design as it explicitly enables to set the convergence speed of the Lyapunov function. Another interesting convergence feature can be deduced from ( 21): a faster convergence is obtained for larger values of . This implies that hyperbolic systems with high convective velocities converge faster, which is physically consistent.

The following corollary gives a sufficient condition for the design of a stabilizing controller for the particular case where C c is a diagonal matrix, D c = 0 and n = n x .

Corollary 1. [Design of a stabilizing controller] Along with Assumption 1, if C c is diagonal and D c = 0 and if there exists a diagonal positive definite matrix Q ∈ R n×n and a scalar µ > 0 such that the following linear matrix inequality is satisfied, for all i ∈ 1, ..., N ϕ ,

  QA c (w i ) T + A c (w i )Q + C c Λ(w i )QC c + µΛ(w i )Q B c (w i )Y Y T B c (w i ) T -e -µ Λ(w i )Q   0 ( 25
)
where Y = KQ, then there exist two constants α > 0 and M > 0 such that, for all ξ 0 ∈ L 2 ((0, 1); R n ) and X 0 c ∈ R n , the solution of (1) with boundary conditions (3) and initial condition (5) satisfies (15) for all t ≥ 0.

Proof: The proof of Corollary 1 can be found in [START_REF] Castillo | Dynamic boundary stabilization of linear parameter varying hyperbolic systems: Application to a Poiseuille flow[END_REF].

This corollary is interesting because, for systems where C C is diagonal and D c = 0, it provides a constructive approach to obtain the boundary control gain K using convex optimization algorithms after determining a suitable value of µ > 0 (e.g. chosen to obtain a good performance versus robustness trade-off for the system considered).

III. FRESH AIR MASS FRACTION CONTROL

For the air fraction boundary control design problem, consider the admission air-path of a Diesel engine with low pressure exhaust gas recirculation (LP-EGR) presented in Figure 1. F egrl depends on the engine operating conditions and is considered here as a known exogenous input since this quantity is measured in production engines. The compressor increases the enthalpy of the gas which results in an increase of the compressor downstream (subscript d) gas pressure and temperature p dc and T dc , respectively (see Figure 1). To increase the gas mass in the engine cylinders and therefore the engine power, the gas coming from the compressor is cooled down by a heat exchanger, which increases the gas density and therefore the mass inside the cylinders. The pressure and temperature downstream of the heat exchanger are denoted as p dhe and T dhe , respectively. Finally, the gas travels from the heat exchanger to the intake manifold (subscript im) where the gas is delivered to the engine with a mass flow rate Q eng and air fraction F im . In this section, we consider the problem of controlling the air fraction in the intake manifold by regulating the LP-EGR mass flow rate. We focus on the generation of the target Q egrl , which then is translated into LP-EGR valve position by a secondary control.

Since the EGR valves potision controllers have been already implemented in modern engines and is an application well known in the engine control community, then we do not detail its functionally in this work.

The admission air-path presented in Figure 1 can be considered as a series of tubes that are linked together by volumes, i.e. by the volumes upstream and downstream of the heat exchanger V uhe and V dhe , respectively. Using this formulation, we can define three tube sections and three control volumes. The three tube sections are defined as: the tube downstream of the compressor (up to the heat exchanger), the heat exchanger with N d tubes in parallel and the tube downstream of the heat exchanger. The control volumes are: the tube upstream of the compressor (denoted as V uc where gas mixture occurs), the heat exchanger upstream volume V uhe and downstream volume V dhe . The control volumes allow us to formulate a zero-dimensional model for the interfaces between the components of the air-path where most of the mass transport phenomenon occurs (the three tubes). In other words, we do not model the evolution of the states along the space dimensions in these volumes. Instead, we describe the dynamics of the flow using a volumeaverage formulation which eases the modeling of the air-path without affecting significantly the accuracy of the mass transport representation.

The one-dimensional modeling of the engine admission air-path is particularly complex as it involves solving the 1D Euler equations for a compressible gas [START_REF] Winterbone | Theory of Engine Manifold Design: Wave Action Methods for IC Engines[END_REF], which is not appropriate for designing real-time feedback strategies. Therefore, the following assumptions are made to simplify the complexity of the air fraction model, making it suitable for Diesel engine control.

Assumption 2: We assume the following:

A-21: the pressure and particle speed dynamics are much faster than the air fraction dynamics;

A-22: no friction and thermal losses are considered (except for the heat exchanger);

A-23: the speed density equation (equation that accurately predicts the amount of air ingested by an engine during the induction stroke [START_REF] Wang | Hybrid robust control for engines running low temperature combustion and conventional Diesel combustion modes[END_REF]) is considered to model the engine intake mass flow rate. No pulsating flow is considered;

A-24: the cross-sectional area of each of the air-path sections (tube between two devices) is considered constant along the space variable x.

The air fraction in a control volume is modeled by considering the classical 0D approach [START_REF] Wang | Air fraction estimation for multiple combustion mode Diesel engines with dualloop EGR systems[END_REF]:

Ḟv = RT v p v V v n i i=1 (Q i F i ) - no o=1 (Q o F v ) (26)
where F v , T v , p v , V v are the control volume air fraction, temperature, pressure and volume, respectively. n i and n o are the amount of mass inputs and outputs and Q i and Q o are the input and output mass flow rates, respectively. Assumption A-21 implies that

Q eng = Q air + Q egrl .
Therefore, the following equations describe the air mass fraction dynamics in the tube upstream of the compressor and in both heat exchanger volumes (according to [START_REF] Di Meglio | Backstepping stabilization of an underactuated 3x3 linear hyperbolic system of fluid flow transport equations[END_REF]):

Ḟuc = RT uc p uc V uc (-(Q egrl + Q air )F uc + Q egrl F egrl + Q air ) (27)
Ḟuhe = RT dc p dc V uhe (-(Q egrl + Q air )F uhe + (Q egrl + Q air )F dc (L dc , t)) (28)
Ḟdhe = RT dhe p dhe V dhe (-(Q egrl + Q air )F dhe + (Q egrl + Q air )F he (L he , t)) (29)
where F uc , T uc and p uc are the air fraction, temperature, pressure and volume upstream of the compressor. F uhe and F dhe are the air fractions upstream and downstream of the heat exchanger, respectively. F he is the air fraction inside the heat exchanger and F dc (L dc , t) and F he (L he , t)

are the output air fractions of the compressor downstream tube and the heat exchanger tubes, respectively. Note that the dynamics of ( 27) -( 29) are defined in terms of the states with the exception of ( 27) that depends on the time-varying additive term Q egrl F egrl + Q air . This issue is addressed in the sequel to obtain a system in the form of the dynamic boundary conditions [START_REF] Angelis | System Analysis, Modelling and Control with Polytopic Linear Models[END_REF].

For the 1D model of the air fraction dynamics in the tube sections, consider the change in mass fraction across a control volume of length ∂x, which can be expressed as [START_REF] Winterbone | Theory of Engine Manifold Design: Wave Action Methods for IC Engines[END_REF]:

∂ t [ρ(x, t)A t (x)F (x, t)] + ∂ x [ρ(x, t)u(x, t)A t (x)F (x, t)] = 0, ∀ x ∈ [0, L], t ≥ 0 ( 30
)
where L is the length of the pipe, A t the tube cross section, ρ the gas density and u the speed of particles. Equation ( 30) can be expressed as (e.g. see [START_REF] Winterbone | Theory of Engine Manifold Design: Wave Action Methods for IC Engines[END_REF] Chapter 4.2.5 for more details):

∂ t [ρ(x, t)F (x, t)] + ∂ x [ρ(x, t)u(x, t)F (x, t)] + ρ(x, t)u(x, t)F (x, t)∂ x [ln(A t (x))] = 0 (31)
Expanding [START_REF] Santos | Boundary control of open channels with numerical and experimental validations[END_REF] and dividing by F (x, t) gives:

∂ t ρ(x, t) + u(x, t)∂ x ρ(x, t) + ρ(x, t)∂ x u(x, t) + ρ(x, t)u(x, t)∂ x (ln(A t (x))) + ρ(x, t) F (x, t) ∂ t F (x, t) + ρ(x, t)u(x, t) F (x, t) ∂ x F (x, t) = 0 (32)
The first four terms of (32) constitute the global continuity of mass implying, that the sum of the four terms equals zero. Therefore, the air fraction continuity equation implies that:

∂ t F (x, t) + u(x, t)∂ x F (x, t) = 0 F (0, t) = F in (t) F (x, 0) = F 0 (x), ∀ x ∈ [0, L], t ≥ 0 (33)
where F in is the air fraction at the input boundary condition, F 0 (x) is a continuous differentiable function describing the initial conditions and u(t, x) is the air fraction propagation speed (independent of the air fraction).

The assumptions A-21, A-22 and A-24 imply that u(x, t) is only time-varying (constant in space), which simplifies significantly the solution of [START_REF] Tannehill | Computational Fluid Mechanics and Heat Transfer[END_REF]. For the heat exchanger, where A-22 does not apply, the average temperature and pressure are defined as follows:

T he = T dc + T dhe 2 , p he = p dc + p dhe 2 (34)
where T he and p he are the heat exchanger temperature and pressure, respectively. Uniform pressure and temperature distribution are assumed inside the heat exchanger to obtain [START_REF] Wang | Hybrid robust control for engines running low temperature combustion and conventional Diesel combustion modes[END_REF], which allows approximating a space average particle speed inside the heat exchanger as presented later in (38).

From ( 33) and ( 34), the dynamics of each tube section can be modeled with the following set of time-varying first order hyperbolic partial differential equations:

∂ t F dc + u dc (t)∂ x F dc = 0, F dc (0, t) = F uc (t), F dc (x, 0) = F dc0 (x) (35)
∂ t F he + u he (t)∂ x F he = 0, F he (0, t) = F uhe (t), F he (x, 0) = F he0 (x) (36)
∂ t F dhe + u dhe (t)∂ x F dhe = 0, F dhe (0, t) = F dhe (t), F dhe (0, x) = F de0 (x) (37)
for all x ∈ [0, 1] and t ≥ 0. The particle speeds for each air-path section u dc , u he and u dhe can be normalized (implying that x ∈ [0, 1]) and calculated using the ideal gas law as follows:

u dc = RT dc (Q air + Q egrl ) p dc A dc L dc , u he = RT he (Q air + Q egrl ) p he A he L he N d , u dhe = RT dhe (Q air + Q egrl ) p dhe A dhe L dhe ( 38
)
where L dc , L he and L dhe are the lengths, A dc , A he and A dhe are the cross-sectional areas of each of the respective tube sections and N d is the number of parallel tubes in the heat exchanger. To use the results of Theorem 1, the system has to be in the form of (1) with boundary conditions [START_REF] Angelis | System Analysis, Modelling and Control with Polytopic Linear Models[END_REF]. System ( 35) -( 37) is a cascade of hyperbolic systems connected by the dynamics defined by ( 27) -( 29), which does not correspond to the system structure (1). This problem is solved by following a similar approach as in [START_REF] De Halleux | Boundary Control of Quasi-Linear Hyperbolic Initial Boundary-Value Problems[END_REF], where n hyperbolic systems in cascade form are combined in one PDE of order n by re-defining the boundary control. However, it is important to highlight that our realization takes into account the dynamic behavior of the boundary conditions, therefore a different stability analysis as to be taken into account. Consider the definition of the air fraction errors as:

ξ dc (x, t) = F dc (x, t) -F imref , ξ he (x, t) = F he (x, t) -F imref , ξ dhe (x, t) = F dc (x, t) -F imref , (39)
and

ξ uc (t) = F uc (t) -F imref , ξ uhe (t) = F uhe (t) -F imref , ξ dhe (t) = F dhe (t) -F imref (40)
where F imref is a given scalar reference for the air fraction. According to (39) and (40), we can rewrite the dynamics of the simplified air fraction air-path model ( 27) -( 29) and ( 35) -( 37) as an LPV hyperbolic PDE of order n as follows:

∂ t      ξ dc ξ he ξ dhe      +      ϕ 2 A dc L dc 0 0 0 ϕ 3 A he L he N d 0 0 0 ϕ 4 A dhe L dhe      ∂ x      ξ dc ξ he ξ dhe      = 0 (41)
with the dynamic boundary conditions in the form of (3):

     ξuc ξuhe ξdhe      =      -ϕ 1 Vuc 0 0 0 -ϕ 2 V uhe 0 0 0 -ϕ 4 V dhe           ξ uc ξ uhe ξ dhe      +      ϕ 1 Vuc 0 0 0 ϕ 2 V uhe 0 0 0 ϕ 4 V dhe           ṽ ξ dc (1)
ξ he (1)      (42)
The boundary control in the form of (4) that allows writing the cascade of hyperbolic systems as a PDE of order n is:

     ṽ ξ dc (1) ξ he (1)      =      0 0 K 13 1 0 0 0 1 0           ξ dc (1)
ξ he (1)

ξ dhe (1)      (43)
where ṽ is a virtual control input that cancels the additive terms of ( 27) and is defined as:

ṽ = F egrl Q egrl Q air + Q egrl -F imref - Q air Q air + Q egrl (44)
The time-varying parameters ϕ i are given by:

ϕ 1 = RT uc (Q air + Q egrl ) p uc , ϕ 2 = RT dc (Q air + Q egrl ) p dc , ϕ 3 = RT he (Q air + Q egrl ) p he , ϕ 4 = RT dhe (Q air + Q egrl ) p dhe (45)
From ( 43) and ( 44), the LP-EGR mass flow rate is defined as:

Q egrl = Q air (F imref -1 + K 13 (ξ dhe (1))) F egrl -F imref -K 13 (ξ dhe (1)) (46)
The LP-EGR mass flow rate target of ( 46) is then transformed into LP-EGR valve position by a secondary controller which typically is composed by an open-loop control obtained from an inverted orifice equation and a closed-loop control composed by a PID and a LP-EGR mass flow rate estimation as the one presented in [START_REF] Castillo | Simultaneous air fraction and low-pressure EGR mass flow rate estimation for Diesel engines[END_REF]. Figure 2 presents the final schematic of the proposed air fraction control in the intake manifold. The system (41) -( 45) with boundary control (43) is in the adequate form to apply the results of Theorem 1. However, with the definition of the varying parameters presented in (45), the polytope (convex set Z ϕ ) obtained by using the polytopic linear representations (PLR) ( 12) and ( 13) is very large and therefore highly conservative. Since all the varying parameters in (45) are strongly dependent on the compressor mass flow rate Q air + Q egrl , specially from the fact that this flow has a direct impact on the mass transport time in the admission path, it is highly convenient to use this quantity to reduce the size of the polytope and therefore the conservatism of the results. The engine speed could be also used to represent the dynamics of the admission path [START_REF] Chauvin | Motion planning for experimental airpath control of a Diesel homogeneous charge-compression ignition engine[END_REF]; however, it is more complicated in this application because other factors have to be taken into account to express the mass transport phenomenon in terms of N eng . Among these factors we find the high-pressure EGR, the swirl and the engine volumetric efficiency.

Let us redefine the varying parameters as:

φ(t) = Q air + Q egrl , (Q air + Q egrl ) 2 (47)
and the original varying parameters ϕ as a linear combination of φ as follows:

ϕ j (t) = ϑ T j ψ(t), j ∈ [1, ..., 4] (48)
where ϑ j ∈ R 3 is an unknown coefficient vector and

ψ(t) = [1, (Q air + Q egrl ), (Q air + Q egrl ) 2 ].
A classical least square method can be used to determine the optimal set of coefficients ϑ j such that |ϕ j (t) -ϑ T j ψ(t)| 2 is minimized and the coefficient identification can be done using engine benchmark measurements over representative engine operating conditions. With this approach, the amount of varying parameters is decreased from four to two, which significantly reduces the conservatism of the control synthesis.

In Figure 3a, the polytope Z ϕ formed by all the extremities of φ1 and φ2 is shown. As proposed demonstrated in [START_REF] Robert | An H infinity LPV design for sampling varying controllers: Experimentation with a T-inverted pendulum[END_REF], the conservatism can be further reduced from the fact that φ2 = φ2 1 , which allows considering only the polytope formed by three vertexes for the control synthesis, as shown in Figure 3b.

The bounds of the varying parameter vector φ can also be found using experimental measurements over representative engine operating conditions.

To summarize, the air fraction is controlled using the results of Section II over the polytope formed by the extremities of φ to find the gain K 13 in (46) such that the air fraction at the engine intake manifold F im converges in finite time to the reference F imref . In the next section, simulation results illustrate the performance of the proposed air fraction control.

IV. SIMULATION RESULTS

To illustrate the effectiveness of the proposed LPV boundary control strategy, an air fraction control (46) that satisfies Theorem 1 is simulated along with the setup shown in Figure 1.

The simulation is performed using a 1D model based on the Euler equations for compressible gases, which is numerically solved with a MacCormack numerical scheme and a Time Variation Diminishing strategy (for more detail on the numerical methods refer to [START_REF] Winterbone | Theory of Engine Manifold Design: Wave Action Methods for IC Engines[END_REF][START_REF] Chung | Computational Fluid Dynamics[END_REF][START_REF] Tannehill | Computational Fluid Mechanics and Heat Transfer[END_REF]). This 1D model has been experimentally validated for inflow, outflow and intra-pipe restrictions [START_REF] Castillo | Intrapipe restriction non-homentropic boundary resolution method[END_REF] [12] [START_REF] Martin | 0-D -1-D Modeling of the air path of ICE Engines for control purposes[END_REF]. Table I gives the physical and numerical parameters used for the simulation of a light-duty 1.6 liter Diesel engine.

The bounds on each parameter as well as the coefficients of (48) are found from measurements on an engine benchmark operated over a wide range of engine operating conditions. More precisely, we use the new motor vehicle emissions group cycle along with two additional engine cycles. The identification of the coefficient vector ϑ i is performed using a traditional least-squaredmethod and the Penrose-Moore pseudo-inverse. Table II presents the identified root-mean-square deviation (RMSD) of ϑ i , where RM SD j = RM SD(ϕ j (t) -ϑ T j ψ(t)). As shown in Table II, the parametrization (48) is representative of the engine varying parameters ϕ and is therefore adequate for the boundary control synthesis.

Consider, as a criterion for the boundary control design, a maximum air fraction convergence time of 1.2 s, (i.e. the time to reach 90% of the final asymptotic value). From ( 21) and the smallest eigenvalue of Λ(ϕ), we find that µ > 0.6 is required to achieve the desired convergence time. Solving ( 14) over the reduced polytope (presented in Figure 3b), leads to the following control gain and Lyapunov matrices P 1 and P 2 : As depicted in Figure 4, the intake manifold air fraction converges to the reference while respecting the convergence time criterion for the three operating conditions. The convergence time depends on the mass flow rate going through the air system (see Figure 4), which is due to the intrinsic time delays associated with transport in the system. In Figure 5 the LP-EGR mass flow rate is presented for the three engine operating conditions: when the mass flow rate increases in the system, Q egrl also does, to maintain the appropriate air fraction.

The convergence of the Lyapunov function ( 16) is illustrated on Figure 6, where it is shown that the actual system convergence is always faster than the exponential decrease predicted by [START_REF] Hante | Converse Lyapunov theorems for switched systems in Banach and Hilbert spaces[END_REF]. To evaluate the efficiency of the controller under strong engine transient conditions, a simulation of an engine cycle with abrupt changes in fuel quantities and air fraction in the exhaust and intake manifold has been performed. Figure 7 shows the engine operating conditions during the proposed maneuver. Figures 8 and9 present the evolution in time of the air fraction in the intake manifold and the LP-EGR mass flow rate, respectively. As depicted, the air fraction in the intake manifold effectively tracks the air fraction target F imref in a smooth trajectory, only separated by the transport time inherent of the system. Figure 9 shows that also the Q egrl effectively compensates the variation of the engine operating conditions, enabling the controller to properly control the air fraction in the intake manifold.

Due to the low calculation load and low complexity of the boundary control (46), this approach is suitable for real-time implementation in a production engine. The current limitation to implement this strategy mainly comes from the lack of reliable measurements or estimations of the fresh air mass fraction in the intake manifold under strong transient conditions (e.g. lack of a reliable sensor for detecting the mass transport). However, many efforts are being done to address this issue by using air fraction observers that take into account the mass transport phenomenon (e.g. [START_REF] Castillo | Boundary observers for linear and quasilinear hyperbolic systems with application to flow control[END_REF]), which would enable the implementation of the proposed boundary control in production engines.

V. CONCLUSIONS AND PERSPECTIVES

This paper proposes a new solution for the boundary control of the fresh air mass fraction in a Diesel engine operated with low-pressure exhaust gas recirculation. A formulation using a cascade of first order LPV hyperbolic systems with dynamics associated with the boundary conditions has been considered to model the air fraction transport phenomenon. Some sufficient conditions for the exponential stability of such systems were obtained using a matrix inequality approach. The exponential stability was demonstrated by means of a strict Lyapunov function formulation, along with a linear polytopic representation for the LPV hyperbolic system. By defining a minimal air fraction convergence time as the performance criterion, an air fraction regulation strategy was designed using the previous stability results. The conservatism of the boundary feedback control was significantly reduced by redefining the system varying parameters.

The complexity of the obtained boundary control is very low, which makes it suitable for real-time implementation. A 1D simulation of the engine admission air-path with the proposed controller was performed to evaluate the effectiveness of the control approach. The simulation results are promising and motivate future steps toward implementation.

The new control design method presented in this paper has many applications in various systems governed by hyperbolic PDE's, such as mine ventilation, traffic control and hydraulic networks.

However, many questions are still open. In particular, a generalization of some sufficient conditions for the exponential stability of hyperbolic systems with non-linear dynamic boundary conditions seems to be a challenging issue. Considering perturbations in system (1) and in the boundary conditions dynamics (3) is also a natural extension of this work.

Fig. 1 :

 1 Fig. 1: Schematic of a Diesel engine admission air-path

Fig. 2 :

 2 Fig. 2: Resulting air fraction control architecture

Fig. 3 :

 3 Fig. 3: Parameter set reduction

K 13 =Figures 4

 134 Figures4 and 5show the simulation results obtained from a change of air fraction reference (F imref ) from 1 to 0.7 at time 0.1 s. Three different air-path operating conditions are simulated, in order to verify the robustness of the controller with respect to parameter variation and the impact of the mass transport phenomenon in the control of the air fraction in the intake manifold.

Fig. 4 :

 4 Fig. 4: Intake manifold air fraction simulation results for a change of air fraction reference

Fig. 5 :Fig. 6 :

 56 Fig. 5: LP-EGR mass flow rate simulation results for a change of air fraction reference

Fig. 7 :

 7 Fig. 7: Time evolution of the engine operating conditions

Fig. 8 :

 8 Fig. 8: Time evolution of the air fraction in the intake manifold

Fig. 9 :

 9 Fig. 9: Time evolution of the LP-EGR mas flow rate

TABLE I :

 I Simulation parameters. dt is the simulation time step and N uhe , N he and N dhe are the discretization sizes of the tube upstream the compressor, heat exchanger and downstream the compressor, respectively

	L dc	1.3 m	A he	6 × 10 -4 m 2
	A dc	0.002 m 2	V uhe	5 × 10 -4 m 3
	Vuc	3 × 10 -4 m 3	V dhe	5 × 10 -4 m 3
	N d	8	L dhe	1.1 m
	L he	0.58 m	A dhe	0.002 m 2
	dt	40 µs	N uhe	50
	N he	20	N dhe	50
	The bounds obtained for φ are the following:
	φ1 ∈ [0.009, 0.074], φ2 ∈ [0.00009, 0.0054]	(49)

TABLE II :

 II The RMSD obtained after the solution of the least square method