Camille Cany
email: camille.cany@cea.fr

Christine Mansilla

Pascal Da Costa

Gilles Mathonnière

Thierry Duquesnoy

Anne Baschwitz

Nuclear and non-dispatchable renewables: two compatible supply options? The case of the French power mix

Keywords: Nuclear power flexibility, non-dispatchable renewable energy sources, Power system economics, back-up, low-carbon source synergy

Introduction

The current international context is characterised by emerging intentions to switch to lowcarbon energy mixes, with country-specific energy transition pathways. The production of heat and electricity is the first contributor to greenhouse gases worldwide as it emitted a quarter of the total emissions in 2010 (French Ministry for the Environment, Sustainable Development and Energy, 2015). Thus electricity production appears to be a key parameter in working towards lower carbon contents. As stated in the SET Plan Integrated Roadmap of the European Commission: "The decarbonisation of electricity production is the centre-piece of the Energy Roadmap 2050. All scenarios studied in the Roadmap show that electricity will have to play a much greater role than now" (European Commission, 2014a). The electric power mix is a core issue of the energy transition: significant decarbonisation of the energy system will be driven by both decarbonising the power sector and enhancing the role of electricity, particularly through sector coupling, like power-to-heat and power-to-mobility, either directly with electricity, or synthetic gas as final energy.

To promote decarbonisation of the power system, the general 3X20 European directive proposes renewable penetration goals (European Union, 2009). The recent European agreement appears proactive; by announcing a binding target of at least 27% of renewable energy on a European level by 2030, it will promote such a development (European Commission, 2014b). In France, 27% of the electricity is to be produced by renewable resources by 2020 (European Union, 2009), and this share will grow up to 2050. In 2013, the share of renewable power already reached 18% of the domestic production, namely approximately 100 TWh [START_REF] Rte | Bilan prévisionnel de l'équilibre offre-demande d'électricité en France[END_REF]. The major part of the growth will come from nondispatchable renewables, which challenges the possibility of maintaining the reliability target level of the power system [START_REF] Gross | The costs and impacts of intermittency: an assessment of the evidence on the costs and impacts of intermittent generation on the British electricity network[END_REF][START_REF] Hart | The Potential of Intermittent Renewables to Meet Electric Power Demand: Current Methods and Emerging Analytical Techniques[END_REF]. Wind and solar are expected to contribute to about 10% of the French electricity production in 2020 [START_REF] Ancre | Scénarios de l'ANCRE pour la transition énergétique[END_REF][START_REF] Rte | Bilan prévisionnel de l'équilibre offre-demande d'électricité en France[END_REF], and according to proactive scenarios, they could contribute to over 50% of the total electricity production by 2050 (ADEME, 2013).

To ensure the reliability target level of a power system, some power plants have to certify that they are available to supply power in case of power deviation from the expected value. This is what is commonly called 'back-up power'. The addition of intermittent power plants in a power system triggers new needs for back-up power, both in the short term (the operational back-up meets balancing requirements) and in the long term (the capacity back-up meets adequacy requirements) [START_REF] Luickx | Considerations on the backup of wind power: Operational backup[END_REF].

Intermittent power plants are characterised by power variability, some uncertainty and nondispatchability with, what is more, a current priority dispatch. In the last few years, the findings first were concerned with defining and quantifying the specificities of intermittent system-dependent production profiles [START_REF] Hart | The Potential of Intermittent Renewables to Meet Electric Power Demand: Current Methods and Emerging Analytical Techniques[END_REF][START_REF] Keppler | Nuclear energy and renewables: system effects in lowcarbon electricity systems[END_REF][START_REF] Luickx | Considerations on the backup of wind power: Operational backup[END_REF][START_REF] Perez-Arriaga | impacts of intermittent renewables on electricity generation system operation[END_REF][START_REF] Ruiz Gomez | Intégration de la production éolienne aux réseaux électriques : approches techniques et économiques[END_REF][START_REF] Wagner | Feature of an electricity supply system based on variable input[END_REF][START_REF] Wan | Analysis of Wind Power Ramping Behavior in ERCOT[END_REF].

They identified average load factors, attempted to quantify additional power gradients and amplitudes induced for the residual load pattern, and examined the issue of power surplus.

A majority of studies on the power system focuses on the impact of non-dispatchable renewable power plants on the power system, and the related needs for more operational and capacity back-up. In particular, the capacity credit of the non-dispatchable sources added to the system must be assessed in order to quantify the need for capacity back-up. Both back-up needs depend on key parameters pointed out by findings shared by [START_REF] Davis | Determining the impact of wind on system costs via the temporal patterns of load and wind generation[END_REF][START_REF] Doherty | Establishing the role that wind generation may have in future generation portfolios[END_REF][START_REF] Gross | The costs and impacts of intermittency: an assessment of the evidence on the costs and impacts of intermittent generation on the British electricity network[END_REF][START_REF] Hand | Renewable Electricity Futures Study[END_REF][START_REF] Hoogwijk | Exploring the impact on cost and electricity production of high penetration levels of intermittent electricity in OECD Europe and the USA, results for wind energy[END_REF][START_REF] Keane | Capacity Value of Wind Power[END_REF][START_REF] Perez-Arriaga | impacts of intermittent renewables on electricity generation system operation[END_REF]. These parameters are inherent to a power system (i.e. relative to the reliability target level, to the features of the power system before introducing intermittent sources and to the features of the intermittent fleet added to the system) on the one hand, and highly linked to approximations used to evaluate these needs on the other hand.

According to [START_REF] Keane | Capacity Value of Wind Power[END_REF], calculations must rely on time series of data for the electric demand that must coincide with the production of non-dispatchable renewable sources, covering at least several years with an hourly time frame, along with a complete inventory of dispatchable sources, associated default rates and maintenance schedules.

Costs related to intermittency are highly sensitive to these parameters, so their quantification should be interpreted with caution, all the more as different terminologies are used from one country to another. These costs are the results of low capacity credits of intermittent sources.

Authors (Connolly et al., 2010a;[START_REF] Gross | The costs and impacts of intermittency: an assessment of the evidence on the costs and impacts of intermittent generation on the British electricity network[END_REF][START_REF] Keppler | Nuclear energy and renewables: system effects in lowcarbon electricity systems[END_REF][START_REF] Perez-Arriaga | impacts of intermittent renewables on electricity generation system operation[END_REF][START_REF] Skea | Intermittent renewable generation and the cost of maintaining power system reliability[END_REF][START_REF] Ueckerdt | System LCOE: What are the costs of variable renewables?[END_REF] agree on the three components which define them:

-Balancing costs that include the change of load factors for the installed capacities.

They result from new operational needs linked to variability and uncertainty of these sources;

-Adequacy costs that are closely related to the assumed capacity factor of the power plants. They result from new needs for capacity back-up; -Grid costs that result from new needs for network reinforcement. [START_REF] Ueckerdt | System LCOE: What are the costs of variable renewables?[END_REF] propose to include all these cost components in the concept of 'system LCOE'1 to make it possible to perform a thorough comparison of technologies.

The French case is very specific. The French power system is currently characterised by a high nuclear penetration: it supplied 82% of the French domestic consumption and 73% of the total demand (including exportations) in 2013 [START_REF] Rte | Bilan prévisionnel de l'équilibre offre-demande d'électricité en France[END_REF]. Nuclear power will remain a significant contributor to the French power system in the medium term, as a low-carbon power source. The nuclear share is to be reduced to 50% of the power production from 2025 onwards, and the renewable share should reach 40% by 2030 (French Government, 2015).

Besides, the choice of the nuclear fleet replacement policy is at the core of the French power debate as half of the fleet will be older than 40 years by 2025.

The energy transition in the French context needs to be pursued by taking into account available technologies and by implementing potential synergies to drive low-carbon power sources with complementary features. With high non-dispatchable renewable penetration, all back-up options should be considered, given their characteristics. Besides peaking unit production -such as gas turbines -which is usually the preferred option, storage, demand (or supply) curtailment, interconnections and even baseload power modulation should be examined as part of the solution [START_REF] Hand | Renewable Electricity Futures Study[END_REF]. In this study, we chose to investigate this question with respect to nuclear flexibility as a way of managing intermittency, which is hardly ever considered in the literature. Our proposed methodology is explained in the first part of this article. It involves designing scenarios in terms of nuclear and non-dispatchable renewable penetration levels, and developing residual load duration curves for each of them.

Herein, we examine and discuss the impact of non-dispatchable renewable energy penetration on the nuclear load factors and the LCOE. Besides, nuclear fleet replacement scenarios were developed and analysed as regards the nuclear LCOE for nuclear participating in load following compared to nuclear operating at baseload. Finally, we discuss nuclear modulation when used as a back-up for wind and solar power sources, we compare it to the gas back-up and we elaborate several policy recommendations.

Methodology and data

As earlier stated, the methodology comprises two major components: the scenario design and the construction of residual load duration curves.

Scenario design

Nuclear and non-dispatchable renewable penetration scenarios

With the aim of achieving a low-carbon power mix in the medium to long term, the so-called energy transition needs to take into account national specificities (power mix history and inertia, regional resources and constraints). In this respect, the French Government launched a broad national consultation on the energy transition in 2013 (Conseil National du Débat, 2013). In this framework, several institutions developed energy prospective scenarios for France. They proposed scenarios in terms of renewable and nuclear energy penetration.

The French National Alliance for Energy Research Coordination (ANCRE2 ) proposed three contrasting scenarios that all meet the factor 4, i.e. a four-fold reduction in greenhouse gas emissions by 2050, compared with the 1990 level [START_REF] Ancre | Scénarios de l'ANCRE pour la transition énergétique[END_REF]. The specificity of these scenarios is a proactive approach to scientific and technological innovation. All these scenarios propose a similar wind and solar penetration for a given time horizon, leading to a maximum energy penetration of 30% by 2050. In order to consider a more proactive assumption in terms of non-dispatchable renewable penetration, we also considered the scenario proposed by the French Environment and Energy Management Agency (ADEME3 ), which examines the case of 50% non-dispatchable renewable penetration by 2050 (ADEME, 2013) . We deliberately focused on nuclear and non-dispatchable renewable resources in each case without clarifying the remaining part of the energy production, which is beyond the scope of this study. However, an explanation should be given in what follows. Today, approximately 13% of the French power production is supplied by hydro power; this figure is to remain stable in the years to come. By assuming 16.5% renewables in 2030, the 'med' scenario would imply that at least 15% of the mix should be produced by fossil resources compared to a 5% value today, hence a higher carbon-content of the power. This led us to consider a higher nuclear share in the 'high nuc' scenario, in which the fossil resources would be reduced compared to today. To enrich the comparisons, we also took into account a 'businessas-usual' nuclear value of 75%, which reflects the current nuclear contribution to the French annual electricity production [START_REF] Rte | Bilan prévisionnel de l'équilibre offre-demande d'électricité en France[END_REF].

In what follows, results will be presented either through the penetration level of nondispatchable renewable and nuclear sources (see section 3.1.), or through the scenario name (see section 3.2.). Sometimes, nuclear penetration rates will be translated into installed capacities for readability issues (see section 3.1).

Nuclear fleet replacement scenarios

In France, 58 nuclear units are currently in operation. They were built between 1977 and 1999 and were initially designed to last, at least, 40 years [START_REF] Cea | ELECNUC. Nuclear power plants in the world[END_REF]. In 2025, half of the fleet will be older than 40 years. In such a context, the nuclear fleet replacement policy choice is at the core of the power mix transition debate. The ten-year inspection programs aim at reassessing the nuclear plants safety and delivering the authorisation to maintain the plant in operation [START_REF] Lokhov | The economics of long-term operation of nuclear power plants[END_REF]. The additional costs associated with life extension programs have been assessed and explained in the 2014 French Court of Auditors' report on the nuclear costs in the French context (Cour des comptes, 2014). We address this question by designing and analysing contrasting nuclear fleet replacement scenarios with an academic viewpoint.

For each non-dispatchable and nuclear penetration scenario, the influence of the nuclear fleet composition on the nuclear LCOE is assessed on a given date (see section 3.2). To this end, two plausible scenarios have been developed as regards the French nuclear fleet replacement.

The first is driven by the power plants' age, while the second is designed to be progressive in order to smooth the investment requirements:

-Current reactors are operated 40 years and then shut down ('40 yrs' scenario);

-The current fleet is progressively replaced to avoid massive investment during a short period of time: 2 GWe are shut down each year starting from 2017 ('prog.' scenario).

In each case, the current fleet made of Pressurised Water Reactors (PWR) is supplemented by Evolutionary Power Reactors (EPR™) in order to achieve the target nuclear penetration level of the selected scenario. An example is provided in Figure 1 for a 50% nuclear penetration in terms of produced energy. In this example, the total energy production is assumed to be 540TWh, which is equivalent to the current level (Eco2mix, RTE database, 2015). As a matter of fact, in this study we neither assume a significant increase in the power demand nor a significant reduction. The PWR load factor is assumed to be 0.75. The EPR™ load factor is assumed to evolve from 0.80 between now and 2040, to 0.85 in 2050, thus taking into account improvements in the technology and organisation of the fleet. From these assumptions, we assess the EPRs™ installed capacity according to the selected replacement scenario, as displayed in Figure 1. In the example illustrated in Figure 1 (50% nuclear penetration case), the nuclear contribution to the total energy production of the power mix declines from 75% in 2015 to 50% in 2030,

where it stabilises at this level. Other assumptions will be examined in section 3.

In the 'prog.' scenario, the PWR shut-down is slower than in the '40 yrs' scenario: for a given year, more PWRs are in operation for the same level of nuclear production capacity. This scenario corresponds to the most realistic one in the French context.

Assessing the nuclear power modulation capacity through the construction of residual load duration curves

To assess the impact of non-dispatchable renewable sources on other power plant capacity requirements and annual productions in a power system, two major types of approaches can be distinguished. First, the construction of load net of wind and solar (called 'residual' in what follows) duration curves for different wind and solar penetration levels can be carried out.

This is used to determine the optimal mix of other generation resources in the long term [START_REF] Davis | Determining the impact of wind on system costs via the temporal patterns of load and wind generation[END_REF][START_REF] Keppler | Nuclear energy and renewables: system effects in lowcarbon electricity systems[END_REF]. Secondly, a minimum cost economic dispatch model can be used to refine the general trends obtained from the first approach [START_REF] Barnacle | Modelling generation and infrastructure requirements for transition pathways[END_REF][START_REF] Barton | The evolution of electricity demand and the role for demand side participation, in buildings and transport[END_REF]Connolly et al., 2010b;[START_REF] Davis | Determining the impact of wind on system costs via the temporal patterns of load and wind generation[END_REF][START_REF] Dautremont | Modèle d'Allocations Efficientes des technologies de production de l'électricité (MAEL), Modélisation économique du système de production électrique français[END_REF][START_REF] Keppler | Nuclear energy and renewables: system effects in lowcarbon electricity systems[END_REF]. The use of such a deterministic input/output model makes it possible to balance short term supply and demand, in response to technical, economic and political constraints. Herein, we focus on the first approach to highlight the general trends in the long term.

The developed methodology involves several steps. First, the construction of annual residual load duration curves considered complementarily with screening curves (see sections 2.2.1 and 2.2.2) makes it possible to assess the annual energy productions and optimal installed capacities from a social welfare viewpoint. One output is the need for baseload installed power that will be fulfilled by nuclear power in the French case.

Thereafter, the nuclear modulation capacity is evaluated based on the realistic French nuclear fleet scenarios presented in the previous section. Thanks to the annual residual load duration curves, we are able to evaluate and compare the nuclear Effective Full Power Hours (EFPH) of both options: nuclear power fleet operated as a baseload capacity, or as a combined option:

baseload and cycling capacity. We then perform an economic comparison.

Construction of load duration curves

A load duration curve represents the sorted hourly load of one year, starting with the highest load hour. The residual load duration curve (see Figure 2) is obtained by withdrawing wind and solar production to the French total production on an hourly time step, thus considering that they have a dispatch priority, and then ordering the residual demand from the highest to the lowest value. wind and solar production profiles were extrapolated for renewable penetration rates, as represented in Figure 2. These rates correspond to scenarios developed by ANCRE for 2030 (16.5%) and 2050 (30%), and by ADEME for 2050 (50%) (see Table 1). Even if wind and solar penetrations are aggregated in the displayed rates, separate assumptions were considered and the corresponding production profiles were treated separately in the calculations.

It should also be noted that the considered demand, to which wind and solar production are withdrawn, actually includes the exportations. Thus, exportations are assumed to be exogenous. Additional exportations may be relevant in case of excess capacities for certain hours and if the interconnection capacities allow it.

Assessing the nuclear power modulation capacity

The first step involved in assessing the nuclear power modulation capacity is to define the 'baseload' nuclear power capacity. We can use the screening curve methodology to do this.

The screening curve methodology is used to design prospective optimal mixes from a social welfare viewpoint [START_REF] Keppler | Nuclear energy and renewables: system effects in lowcarbon electricity systems[END_REF]. It consists in plotting the annual cost of the installed capacity according to its utilisation time (see Figure 3). By plotting the curves for each of the available technology, the optimal minimum annual operating time is obtained for each power plant type so as to minimise the total system costs.

The second step involves copying the obtained minimum operating time for each power plant type on the residual load duration curves, so as to assess the corresponding optimal firm power to operate at full power for at least the minimum operating time resulting from the screening curves. Thus, the screening curve methodology seems to provide a very straightforward way for designing an optimal power mix from the viewpoint of the electric system, or social planner.

However, this method is questionable due to its high sensitivity to the retained assumptions.

We chose to illustrate this sensitivity through the nuclear share -our focus in this paper -and two parameters:

-Discount rate: this assumption proves critical for very capex-intensive technologies, such as nuclear;

-Carbon price: it is a key parameter with respect to increasing the economic competitiveness of low-carbon technologies.

According to our calculations, when varying the discount rate of 3 points, the minimum optimal operating time of nuclear power could vary between 3,000 and 8,000 hours according to discount rate and carbon price assumptions. Figure 3 illustrates the sensitivity of the methodology with two cases for a discount rate of 8%: with no carbon price and with a carbon price of $30/tonCO2. In the first case, nuclear could never become competitive; in the second one, a minimum operating time of 5000 hours would be needed for the nuclear to become the lowest-cost technology compared with coal and gas. As a result, the nuclear power that would be installed from an electric system viewpoint would differ depending on the retained assumptions leading to a resulting firm capacity between 40 GW and 55 GW, for instance, for 16.5% renewable penetration.

Given this uncertainty range, we chose an empirical approach based on the mean historical availability factor of nuclear in France (Power Reactor Information System (PRIS) database, 2015). A value of 7,000 hours was selected. Such an approach rather minimises the resulting installed nuclear capacity. This value was used to assess the 'baseload' firm nuclear power (called P baseload hereafter) according to the residual load to supply. The 'baseload' nuclear production (representing an firm capacity of P baseload ) will then be compared with the nuclear production corresponding to the firm nuclear power capacity provided by the selected scenario (i.e. 'med', 'high nuc', 'high ren', see section 2.1). This will be done to assess the available energy for modulation, i.e. the nuclear modulation capacity. This is illustrated in

Results and discussion

By implementing the methodology that was presented we examine: i/ the impact of renewable penetration on the nuclear power load factors, ii/ the impact of the nuclear fleet replacement policy on the economic competitiveness, and iii/ the potential of nuclear flexibility as a backup power option.

Impact of renewable penetration on the nuclear power load factors

At first, we adapted the assessment illustrated in Figure 4 to each scenario and calculated the share of the total annual electricity production (i.e. 540TWh) that would be provided by nuclear, according to the wind and solar penetration, when nuclear is operated at least 7,000

hours per year at full power. The results are displayed in table 3.

Non-dispatchable renewable penetration (%)

Baseload nuclear penetration (%)

16.5 57 30 41 50 17

Table 3: Baseload nuclear contribution to meet the total power production according to the non-dispatchable renewable penetration (% of produced energy)

From this table, for a 16.5% wind and solar share, it can be seen that 57% of the electricity should be supplied by nuclear as 'baseload', according to the retained assumptions. Hence, 50% nuclear penetration (corresponding to the 'med' scenario in 2030), would be underestimated compared with the 'baseload' reference, and thus no nuclear modulation capacity is available in this scenario.

Attributing dispatch priority to non-dispatchable renewable energy actually lowers the operating time of already existing nuclear reactors, in the short term. This is an issue for capex-intensive technologies such as nuclear. We assessed this impact according to the nuclear power installed capacity and the non-dispatchable renewable penetration level through the use of the residual load curves. The renewable and nuclear penetration assumptions were derived from the scenarios presented in section 2.1. To assess the nuclear installed capacities for the three cases, we assumed a baseload nuclear load factor of 0.75 in accordance with historical data (Power Reactor Information System (PRIS) database, 2015). For the cases with 50% wind and solar, the nuclear annual load factor decreases as low as 50% and up to 40% with a 60GWe installed nuclear capacity (75% share of nuclear).

To provide some orders of magnitude, the nuclear annual revenues were calculated as the produced energy multiplied by the mean value of electricity hourly prices in the day-ahead market in 2014, which equalled 32€/MWh (EPEX SPOT SE: Day-Ahead Auction database, 2014). The figures are presented in Table 5. In Table 6, the 'revenue losses' are presented according to the share of renewable penetration. They are estimated as the difference between the maximum foreseeable revenue (corresponding to the baseload nuclear production, assuming the plants were operated at a load factor of 75%), and the actual revenue that is influenced by the renewable penetration (figures in Table 5). Results in These figures highlight that reduced load factors would represent substantial profit losses in the long run. However, effective losses are expected to be higher because the introduction of renewables with near zero marginal costs will tend to reduce the electricity wholesale market price. Such an effect needs to be taken into account for future works. Nonetheless, significant renewable power penetration does not necessarily mean reducing the nuclear share. Lowcarbon power sources should be promoted in a complementary manner to achieve a lowcarbon power mix. Thus, with a view to make the social planner and plant operator viewpoints coincide, the solution could be to find new outlets rather than reducing the nuclear load factor. This is discussed in section 3.3.

Impact of the nuclear fleet replacement policy on the economic competitiveness

The scenarios presented in section 2.1.2 as regards the power fleet replacement were established based on realistic expectations. They result in different annual LCOE according to the composition of the nuclear mix. The selected assumptions are as follows:

-The annualised investment cost of the nuclear mix was computed from the share between PWRs and recently installed ones (EPR™) for each year of the simulation.

-The annual operation cost is considered constant over the years and has the same value for PWRs and EPRs TM .

-Life extension costs have been added for the 'prog' scenario between 2011 and 2033.

The major part of the additional investment needs falls within 2011 and 2025. As mentioned earlier, nuclear plants safety is reassessed thanks to the ten-year inspection programs. By 2025, the third inspections will be completed for the PWRs of 1300 MWe and more than three-quarters of the fourth inspections will be realised for the PWRs of 900 MWe. What is more, the fourth inspections will end in 2033 for the 900MW and1300MW power plants (Cour des comptes, 2014).

-The latest estimation of the investment cost of the EPR TM under construction in France (Flamanville) was taken from [START_REF] Bezat | Le réacteur EPR de Flamanville touché au coeur[END_REF] with a value of 10 billion€ 2013 . Its life time is assumed to be 60 years.

-A learning rate was assumed for the EPR™; the investment cost was decreased from the first-of-a-kind to the tenth power plant (minus 10% for the second one compared to the first one; minus 9% for the third one compared to the second one, repeating it until minus 1% for the eleventh one compared to the tenth one), and then considered constant.

These cost assumptions are displayed in Table 7 and Table 8 and were taken from the French Court of Auditors' in-depth study of the nuclear costs in the French context (Cour des comptes, 2014;[START_REF] Bezat | Le réacteur EPR de Flamanville touché au coeur[END_REF]. [START_REF] Bezat | Le réacteur EPR de Flamanville touché au coeur[END_REF]. The learning rate assumptions are issued from in-house expertise.

Annual costs per

The nuclear annual LCOE was then computed according to the retained assumptions and the nuclear fleet replacement scenario. To this end, the 'Economic Current Cost' methodology was retained for a discount rate of 8% (actual annualised return between 1971 and 2000), i.e.

the annual LCOE was assessed by annualising the investment including interests over the selected amortisation period (equal to 40 years for the PWRs). After this date, the investment burden is zero: this is the characteristic of this method. Figure 5 -Case 1: for 16.5% renewable penetration in 2030, 'high nuc' scenario (65% nuclear), with a nuclear fleet composition derived from the 'prog.' scenario, costs were compared with the 'baseload' nuclear penetration according to a nuclear fleet composition derived from the '40 yrs' scenario.

-Case 2: for 30% renewable penetration in 2030, 'high ren' scenario (50% nuclear), with a nuclear fleet composition derived from the 'prog.' scenario, costs were compared with the 'baseload' nuclear penetration according to a nuclear fleet composition derived from the '40 yrs' scenario.

Relative differences of nuclear annual LCOE for case 1 and 2 are displayed in Table 9.

Case 1 ("high nuc" / "baseload") -15 % Case 2 ("high ren" / "baseload") -23 % Table 9: Relative difference between the annual nuclear levelised production cost for case 1 and case 2 in 2030

For 16.5% and 30% non-dispatchable renewable penetration in 2030, the 'prog' scenario nuclear replacement policy allows us to consider the possibility of nuclear operating at reduced load factors in order to contribute to load following. In such case, the annual nuclear LCOE remains lower of around 20% than the case with nuclear operating as the 'baseload' with a '40 yrs' replacement policy. The economic advantage of the progressive replacement policy is not compensated by reduced load factors. As a result, under an uncertain context, it is better to promote such a policy, even not presuming of future actual renewable penetration.

This makes it possible to take advantage of the amortised power plants, the history of the French power mix.

Besides enabling a lower production cost for the nuclear fleet, smoothing the replacement of the nuclear capacities would be beneficial to the whole power mix. In the case where nuclear power would be partly replaced by renewables, by somewhat delaying the capacity replacement, it would be possible to 'wait' for cheaper technologies and possibly reduce the overall production cost. In such a case, effective capacity credits should be taken into account according to the nature of the renewable source and the level of penetration into the power mix. Thus, to replace a given capacity of dispatchable power, higher amounts of nondispatchable resources are needed due to intermittency issues. If non-dispatchable renewable production profiles are correlated, their capacity credits tend to decrease with higher penetration levels [START_REF] Keane | Capacity Value of Wind Power[END_REF][START_REF] Keppler | Nuclear energy and renewables: system effects in lowcarbon electricity systems[END_REF].

Nuclear modulation options: leverage to promote decarbonised power mixes

A gap can result from the difference between the installed nuclear capacity and what is considered optimal: either from the electric system viewpoint (with all the attached uncertainties of the screening curve methodology pointed out in section 2.2.2), or from the operator viewpoint that target the maximum load factor. In what follow, nuclear flexibility is examined as a back-up option both in terms of quantity and cost: amount of energy that may be available, and potential economic advantage compared to alternative options.

In this section, the installed nuclear power is compared against the 'baseload' value for each scenario, in order to assess the modulation capacity and energy. For the different cases, Table 10 shows the annual energy dedicated to modulation, its share of the total annual energy produced as well as the corresponding EFPH to be fulfilled by nuclear as back-up power. For a given nuclear penetration level, the share of nuclear energy devoted to modulation increases more than linearly with the renewable penetration level. It is worth pointing out that 50% of nuclear energy would be dedicated to modulation for the 'high ren' scenario in 2050. As stated earlier, the increase in the non-dispatchable renewable share in the power system modifies the residual demand pattern, resulting in higher needs for power ramping and higher amplitude variations, which remain to be quantified [START_REF] Keppler | Nuclear energy and renewables: system effects in lowcarbon electricity systems[END_REF][START_REF] Wagner | Feature of an electricity supply system based on variable input[END_REF][START_REF] Wan | Analysis of Wind Power Ramping Behavior in ERCOT[END_REF].

The French nuclear fleet already contributes to load following because of its high share in the mix. Nuclear power has characteristics compatible with power modulation [START_REF] Bruynooghe | Load-following operation mode at muclear power plants (NPPs) and incidence on operation and maintenance (O&M) costs. Compatibility with wind power variability[END_REF][START_REF] Choho | Major innovations in PWR load follow opearations by AREVA[END_REF][START_REF] Lokhov | Technical and economic aspects of load following with nuclear power plants[END_REF] and the French nuclear fleet has potential for more power ramping and higher amplitude variations according to the power operation track record [START_REF] Bruynooghe | Load-following operation mode at muclear power plants (NPPs) and incidence on operation and maintenance (O&M) costs. Compatibility with wind power variability[END_REF]RTE database -Portail clients -Production réalisée par groupe, 2015). Thanks to the size of the fleet, production could be varied greatly by performing small increments in each power plant. This would remain true as long as the installed nuclear power is sufficient in terms of contributing to the total electricity mix. The quantification of the impacts of intermittent resources on the nuclear hourly production profile (i.e. ramping and amplitudes) will be the topic of future works.

Assuming perfect adaptation to the change in power profile, we compared the nuclear load following option to the gas option, which is considered as the reference back-up option in this study. Since the 'prog' scenario corresponds to the most realistic one in the French context, we will focus on this one. We thus compare:

-The cost of operating the nuclear capacity at a reduced load factor (according to the nuclear and renewable penetrations defined by the scenario: 'med', 'high nuc', 'hig ren'),

-And the cost of operating the nuclear power as a baseload (the nuclear capacity is then only the 'baseload', as defined in section 2.2.2), plus an additional gas capacity to cover what would have been provided by the additional nuclear capacity (defined as the difference between the nuclear penetration hypothesis in the considered scenario and what is required for 'baseload').

LCOE were compared. Cost assumptions as regards nuclear plants were provided by the French Court of Auditors (Cour des comptes, 2014), while the OECD assumptions were taken for the gas-fired power plants [START_REF] Wittenstein | Projected costs of generating electricity[END_REF]. The French assumptions were taken for Combined Cycle Gas Turbines (CCGT) while, due to the lack of data, the German assumptions were taken for the Open Cycle Gas Turbines (OCGT). They are presented in Table 11. The CCGT was selected as the reference, because what would have been provided by the additional nuclear capacity represents annual operation time of more than 1,000 hours per year in our scenarios. For such durations, CCGT are more competitive than OCGT, due to higher OCGT operation costs. Indeed, according to [START_REF] Wittenstein | Projected costs of generating electricity[END_REF], the electrical conversion efficiency of a CCGT equals 61% for the French case and is around 40% for OCGT.

-iii/ if the OCGT was considered as the gas option, the previous results would be modified to around €30/MWh and €5/tonne CO2 respectively. During the first semester of 2015, the mean gas price for industrials in France equalled 35€/MWh without taxes (Eurostat database, 2015). This means that today nuclear would be cheaper than OCGT to provide this back-up.

Even if the CCGT option appears cheaper according to our assumptions, it would generate CO 2 emissions in the range of 10-30 Mt per year, depending on the scenario and the time horizon.

To enlarge the scope, it would be worth considering taking advantage of the available energy to produce valuable services and products [START_REF] Ruth | Nuclear-renewable hybrid energy systems: Opportunities, interconnections, and needs[END_REF], instead of charging power modulation, which would mean charging energy losses. The nuclear economic equation remains unsolved in the case of a new build when the load factor becomes too low. The capital share is too high for efficient use with such low load factors. The question appears more complex if we examine the case of already existing reactors. In systems with high shares of renewables that have a dispatch priority, load factors fall within the range where the competitiveness of nuclear begins to decline when up against gas combined cycles. In particular, these systems would benefit more from opting for the appropriate nuclear fleet replacement policy to compensate from the losses. In the end, the decision of maintaining nuclear power plants with low load factors would be highly linked with payments for balancing services. Future research should address this question in detail.

Thus, we recommend considering nuclear power modulation through the power use. In this case, nuclear power plants would be operated as a baseload and electricity would be supplied to the electric system when requested. For the remaining time, instead of reducing the plant load, the output could be used to produce other valuable energy carriers, such as heat or hydrogen. Hydrogen is also a chemical product that has numerous outlets: ammonia production, refining activities and synthetic fuel production, as well as being a means to enhance connections between energy networks (e.g. power and gas, through power-to-gas).

Finally, hydrogen also makes it possible to integrate high shares of low-carbon energy into a given mix by avoiding both the reducing of the baseload power such as nuclear and the restriction of non-dispatchable renewable energy [START_REF] Caumon | Flexible hydrogen production implementation in the French power system: Expected impacts at the French and European levels[END_REF]. Detailed business models nonetheless remain to be developed in this field. Systems comprising nuclear power plants together with a hydrogen process may propose new energy services such as hydrogen and power supply. What is more, implementing High-Temperature Steam Electrolysis (HTSE) processes in the longer term can make it possible to supply power as a baseload (thanks to the nuclear plant) and for peak demand (thanks to reversible HTSE operation) [START_REF] Haratyk | Nuclear-Renewables Energy System for Hydrogen and Electricity Production[END_REF][START_REF] Mougin | Stack performances in high temperature steam electrolysis and co-electrolysis[END_REF]. During high demand peaks, such a new concept would make it possible to use the power from the nuclear power plant, but also from the reversible electrolyser, thus acting as an energy storage solution and increasing the instantaneous power of the nuclear plant. Such systems could therefore meet both baseload and peak demands. This will be the topic of future research.

Conclusion and policy recommendations

The current French electricity mix is characterised by a high share of nuclear power on the one hand, and the seeking of new energy transition policies on the other hand. Developing low-carbon sources with complementary features is then a central issue which has to be addressed. Nuclear power flexibility is discussed as a way to move in this direction.

In other words, if nuclear power capacities are maintained contrary to political guidelines, it is worth opting for a progressive replacement of installed capacities (PWRs) by EPRs™ in order to cope with the intermittence of non-dispatchable and variable renewable sources. The French power mix specificities are an asset to be maintained if we wish to promote a very-low carbon power mix.

Then, complementarities and even synergies should be encouraged. Despite the fact that nuclear could participate in load following, especially in France since large variations in production could only be achieved by small increments in each power plant, making use of the available energy while operating at baseload would be more profitable. Conversely, even if excess non-dispatchable renewable energy could be curtailed, taking advantage of the excess power seems more fruitful.

For the social planner and plant operator viewpoints to coincide, we therefore recommend modulating the use of nuclear power rather than modulating the production of nuclear power.

This would mean that nuclear power flexibility could accomplish two things with one action: available power would be converted into valuable services to the electric system and, at the same time, into valuable industrial products. In this respect, hydrogen seems especially worthwhile, given its multiple outlets and future specific position within the energy system.

In such a synergistic power system, both nuclear power plant flexibility and nuclear plant operation time extension should be examined, also by considering novel uses and other operational back-up options, such as storage or demand-side management tools, as contributors to the capacity back-up whose optimum solution remains to be established.

Figure 1 :

 1 Figure 1: PWR & EPR™ capacity evolution scenarios (for 50% nuclear penetration)

Figure 2 :

 2 Figure 2: Residual load duration curves in the French case

Figure 3 :

 3 Figure 3: Screening curve example (cost data from Salvadores and Keppler (2010)), 8%

  Figure 4.

Figure 4 :

 4 Figure 4: Assessment of the nuclear power modulation capacity ('high nuc' 2030)

  displays the relative difference between the nuclear annual LCOE for the 'prog.' scenario and the '40 yrs' scenario, for baseload operation. The results are presented for a 50% nuclear penetration in terms of produced energy. The considered period is 2030 -2050: two deadlines considered in the defined scenarios.

Figure 5 :

 5 Figure 5: Relative difference of nuclear annual LCOE (€ 2013 /MWh) between the

Table 1 and

 1 Table 2 and list the assumptions selected for non-dispatchable renewable resources and nuclear share in the total French power production according to the ANCRE and ADEME scenarios.

	Time horizon	Wind penetration (%)	Solar penetration (%)	Total (%)	Reference
	2015	5 %	1 %	6 %	ANCRE, 2013
	2020	8 %	2 %	10 % ANCRE, 2013
	2025	10 %	3 %	13 % ANCRE, 2013
	2030	12 %	4.5 %	16.5 %	ANCRE, 2013
	2040	16 %	7 %	23 % ANCRE, 2013
	2050	20 %	10 %	30 % ANCRE, 2013
	2050	30 %	15 %	45 % ADEME, 2013
	2050	35 %	15 %	50 % ADEME, 2013

Table 1: Non-dispatchable renewable energy penetration scenarios (% of annual produced electricity) Nuclear penetration Equivalent installed capacity (Load factor of 75%) Equivalent annual energy production Reference scenario

	75 %	60 GWe	410 TWh	Business-as-usual
	65%	53 GWe	350 TWh	High nuclear in ANCRE (2013)
				Low nuclear in ANCRE
	50%	40 GWe	270 TWh	(2013) or High nuclear in
				ADEME (2013)

Table 2 : Nuclear energy penetration scenarios

 2

	(% of annual produced electricity)

'High non-dispatchable renewable scenario' (high ren): 30% wind and solar in 2030, 50% in 2050, 50% nuclear in both cases.

Table 4

 4 presents both the total annual nuclear energy production and the nuclear EFPH in each case (a 100% load factor would correspond to 8760 EFPH).

				Installed nuclear capacity (GWe)
			41		53		60
			(50% penetration)	(65% penetration)	(75% penetration)
			TWh	EFPH	TWh	EFPH	TWh	EFPH
	%	16.5	270	6590	330	6300	360	6000
	wind and	30	250	6150	290	5490	310	5120
	solar	50	190	4590	210	3890	210	3530

Table 4 : Resulting nuclear annual energy production and EFPH according to the wind, solar and nuclear penetration levels in the French electricity mix

 4

Table 5 : Annual revenue according to the wind, solar and nuclear penetration levels in the French electricity mix and maximum foreseeable revenue (M€)

 5 Table 6 are expressed as a share of the maximum foreseeable revenue.

				Installed nuclear capacity (GWe)
			41		53	60
			(50% penetration)	(65% penetration)	(75% penetration)
	% wind and solar	16.5 30 50	8640 8060 6020		10690 9310 6590	11520 9820 6780
	Maximum foreseeable revenue	8690		11290	13030
			Installed nuclear capacity (GWe)
			41		53	60
		(50% penetration)	(65% penetration)	(75% penetration)
	% wind	16.5	1		5	12
	and solar	30	7		18	25
		50	31		42	48

Table 6 : Annual revenue losses according to the wind, solar and nuclear installed capacities in the French electricity mix: share of the maximum foreseeable revenue (%)

 6

Table 7 : Assumptions for the cost breakdown of a PWR from

 7 (Cour des comptes, 2014)

	GWe:	M€ 2013 /year/GWe
	Investment	135
	Decommissioning	8
	Operation	158
	Spent fuel and waste management	21
	Total annual costs:	M€ 2013 /year
	Investment for maintenance
	From 2011 to 2025	2300
	From 2026 to 2033	1800
	Life extension
	From 2011 to 2025	2000
	From 2026 to 2033	1800
	Annual costs per GWe:	M€ 2013 /year/GWe
	Investment and decommissioning
	First of a kind	789
	Second	710
	Third	646
	Fourth	594
	Fifth	553
	Sixth	520
	Seventh	494
	Eight	474
	Ninth	460
	Tenth	450
	Eleventh	446
	Operation	158
	Spent fuel and waste management	21

Table 8 : Assumptions for the cost breakdown of a EPR™ from

 8

(Cour des comptes, 2014;

Table 10 : Resulting nuclear annual energy production dedicated to modulation, its share of the total annual energy production and EFPH according to the wind, solar and nuclear penetration levels in the French electricity mix

 10

					Installed nuclear capacity (GWe)
				41			53			60
			(50% penetration)	(65% penetration)	(75% penetration)
			TWh	% EFPH TWh	% EFPH TWh	% EFPH
	% wind	16.5	0	0	0	30	8	6000	50	15	4900
	and solar	30	30	11	6200	70	23	4800	80	27	4200
		50	90	50	4800	110	55	3900	120	56	3400

The system of Levelised Cost Of Electricity.

Alliance Nationale de Coordination de la Recherche pour l'Energie.

Agence de l'Environnement et de la Maîtrise de l'Energie.

Acknowledgements

The authors would like to thank J. C. Bocquet, professor emeritus at CentraleSupelec -Université Paris Saclay, and J. G. Devezeaux de Lavergne, Head of the Institute for technoeconomics of energy systems (I-tésé) for their fruitful discussions and supports.

Considering the nuclear and renewable penetration scenarios defined in section 2.1 ('med', 'high nuc', 'high ren') together with the nuclear fleet replacement scenarios, we must now compare the nuclear fleet costs by considering both effects: i/ the renewable and nuclear energy penetration (impact of load factor reduction as identified in section 3.1), and ii/ the selected fleet replacement (impact of the nuclear fleet composition as identified previously for the 'fleet' cost for baseload).

By 2050, both nuclear replacement scenarios lead to a fleet comprising EPRs™ exclusively, hence costs converge. In what follows, we focus on 2030 to point out differences between the '40 yrs' and 'prog.' scenarios.

On the one hand, in section 3.1, it was highlighted that nuclear load following (particularly resulting from renewable dispatch priority) induced reduced average load factors, hence higher production costs. On the other hand, it was shown in Figure 5 that in 2030, the 'prog.' scenario leads to lower nuclear production costs than the '40 yrs' scenario for baseload operation. Similar trends would obviously be observed when comparing the replacement policies for reduced load factors, since the load factor reduction is only induced by the renewable penetration. Therefore, the question arises: to what extent is the nuclear 'prog.' replacement scenario still cheaper than the '40 yrs' scenario, even when operating the nuclear fleet in a load following mode? Answering this question will put a light on the economic advantage of the progressive replacement policy to participate in load following. To do so, we have used the '40 yrs' replacement scenario with nuclear as the 'baseload' to make our comparison.

We compare annual nuclear LCOE in the two following cases: In each scenario, the CCGT load factor calculation serves determining its LCOE. The load factor is estimated by quantifying the annual energy that has to be provided for modulation and associated installed capacity (calculated with a capacity factor value of 0.85 [START_REF] Wittenstein | Projected costs of generating electricity[END_REF]). Results show that gas would operate with a load factor of 58% to participate in load following by 2030 for the 'high nuc' scenario (respectively, 47% by 2050).

With the retained assumptions, the CCGT cycling option is cheaper than nuclear back-up in each scenario. The gas LCOE which leads to equalising the cost of both options is estimated to be approximately 100€ 2013 /MWh. Thus, to reach this gas LCOE, the associated gas price was determined for a given CO 2 price. Conversely, the same approach was followed for a given gas price.

The results fall within the same orders of magnitude for all scenarios:

-i/ for a CO 2 price of €5/tonne CO2 (which corresponds to the current market value (Thomson Reuters, 2015)), nuclear back-up is competitive when the gas price exceeds €50/MWh;

-ii/ for a gas price of €31/MWh (which corresponds to the ANCRE and OECD assumptionfor 2030 (for 2020 respectively) [START_REF] Ancre | Scénarios de l'ANCRE pour la transition énergétique[END_REF][START_REF] Wittenstein | Projected costs of generating electricity[END_REF], the gas option dominates as long as the CO 2 price is below €100/tonne CO2 ;

At first, the impact of renewable energy penetration on the nuclear load factors and its economic consequences were examined. From the operator's viewpoint, the capital expenses must be amortised, hence the maximum load factor should be targeted. According to installed renewable and nuclear power in the mix, it became clear that significant revenue losses for the nuclear plant operator may result from high penetration rates of renewables.

However, the increasing share of non-dispatchable renewable energy means increasing needs for operational back-up. In economic terms, CCGT gas back-up is more competitive than nuclear in each of the studied scenario. Nonetheless, a carbon tax could smooth the gap.

When compared with CCGT, a value as low as €100/tonne CO2 could change the trends, which is a point in favour of such incentives to avoid greenhouse gas emissions. The French Energy Additionally, nuclear back-up is all the more competitive as amortised power plants are used.

In this transition period, the differences in nuclear production costs between two nuclear fleet replacement policies were assessed with an academic approach. It was shown that a progressive shutdown (i.e. with a longer lifetime of PWRs compared with a shutdown after 40 years of operation) appears to be a relevant choice since it provides a lower nuclear power production cost from 2030 to 2050, even at reduced load factors, while contributing to load following and giving deciders more time before choosing the most sustainable technologies.


Pierre-Olivier Lamare
email: pierre-olivier.lamare@imag.fr

Nikolaos Bekiaris-Liberis
email: bekiaris-liberis@berkeley.edu

Control of 2 × 2 Linear Hyperbolic Systems: Backstepping-Based Trajectory Generation and PI-Based Tracking $

Keywords: Distributed parameter systems, Hyperbolic systems, Backstepping, Lyapunov function

We consider the problems of trajectory generation and tracking for general 2 × 2 systems of firstorder linear hyperbolic PDEs with anti-collocated boundary input and output. We solve the trajectory generation problem via backstepping. The reference input, which generates the desired output, incorporates integral operators acting on advanced and delayed versions of the reference output with kernels which were derived by Vazquez, Krstic, and Coron for the backstepping stabilization of 2 × 2 linear hyperbolic systems. We apply our approach to a wave PDE with indefinite in-domain and boundary damping. For tracking the desired trajectory we employ a PI control law on the tracking error of the output. We prove exponential stability of the closed-loop system, under the proposed PI control law, when the parameters of the plant and the controller satisfy certain conditions, by constructing a novel "non-diagonal" Lyapunov functional. We demonstrate that the proposed PI control law compensates in the output the effect of in-domain and boundary disturbances. We illustrate our results with numerical examples.

Introduction

Control of 2 × 2 systems of first-order hyperbolic PDEs is an active area of research since numerous processes can be modeled with this class of PDE systems. Among various applications, 2 × 2 systems model the dynamics of traffic [START_REF] Fan | Comparative model accuracy of a data-fitted generalized Aw-Rascle-Zhang model[END_REF], [START_REF] Goatin | The Aw-Rascle vehicular traffic flow model with phase transitions[END_REF], hydraulic [START_REF] Bastin | On boundary feedback stabilization of non-uniform linear 2 × 2 hyperbolic systems over a bounded interval[END_REF], [START_REF] De Halleux | Boundary feedback control in networks of open channels[END_REF], [START_REF] Santos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF], [START_REF] Santos | Boundary control of open channels with numerical and experimental validations[END_REF], as well as gas pipeline networks [START_REF] Gugat | Gas flow in fan-shaped networks: Classical solutions and feedback stabilization[END_REF], and the dynamics of transmission lines [START_REF] Curro | A reduction procedure for generalized Riemann problems with application to nonlinear transmission lines[END_REF].

Several articles are dedicated to the control and analysis of 2 × 2 linear [START_REF] Bastin | On boundary feedback stabilization of non-uniform linear 2 × 2 hyperbolic systems over a bounded interval[END_REF], [START_REF] Diagne | Lyapunov exponential stability of linear 1-D hyperbolic systems of balance laws[END_REF], [START_REF] Santos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF], [25] [33], [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2 × 2 linear hyperbolic system[END_REF], [START_REF] Vazquez | Collocated output-feedback stabilization of a 2 × 2 quasilinear hyperbolic system using backstepping[END_REF] and nonlinear [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 one dimensional nonlinear hyperbolic systems[END_REF], [START_REF] Coron | Local exponential H 2 stabilization of a 2 × 2 quasilinear hyperbolic system using backstepping[END_REF], [START_REF] Gugat | Existence of classical solutions and feedback stabilization for the flow in gas networks[END_REF], [START_REF] Prieur | Control of systems of conservation laws with boundary errors[END_REF], [START_REF] Prieur | Robust boundary control of systems of conservation laws[END_REF] systems. Results for the control of n × n systems also exist [START_REF] Di Meglio | Backstepping stabilization of an underactuated 3 × 3 linear hyperbolic system of fluid flow transport equations[END_REF], [START_REF] Di Meglio | Stabilization of a system of n + 1 coupled firstorder hyperbolic linear PDEs with a single boundary input[END_REF], [START_REF] Hu | Control of homodirectional and general heterodirectional linear coupled hyperbolic PDEs[END_REF], [START_REF] Li | Classical solutions for quasilinear hyperbolic systems[END_REF]. Algorithms for disturbance rejection in 2 × 2 systems are recently developed [START_REF] Aamo | Disturbance rejection in 2 × 2 linear hyperbolic systems[END_REF], [START_REF] Tang | Sliding mode control to stabilization of linear 2 × 2 hyperbolic systems with boundary input disturbance[END_REF]. The motion planning problem is solved in [START_REF] Dubois | Motion planning and nonlinear simulations for a tank containing a fluid[END_REF], [START_REF] Petit | Dynamics and solutions to some control problems for water-tank systems[END_REF], for a class of 2 × 2 systems and in [START_REF] Fliess | A distributed parameter approach to the control of a tubular reactor: A multi-variable case[END_REF], [START_REF] Petit | Flatness of heavy chain systems[END_REF] for a class of wave PDEs. Perhaps the most relevant results to the present article are the results in [START_REF] Santos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF], dealing with the Lyapunov-based output-feedback control of 2 × 2 linear systems, the results in [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2 × 2 linear hyperbolic system[END_REF], dealing with the backstepping stabilization of 2 × 2 linear systems, and the results in [START_REF] Petit | Dynamics and solutions to some control problems for water-tank systems[END_REF], dealing with the motion planning for a class of 2 × 2 systems.

In this paper, we are concerned with the trajectory generation and tracking problems for general 2 × 2 systems of first-order linear hyperbolic PDEs with anti-collocated boundary input and output. We solve the motion planning problem for this class of systems employing backstepping (Section 2.1). Specifically, we start from a simple transformed system, namely, a cascade of two first-order hyperbolic PDEs, for which the motion planning problem can be trivially solved. We then apply an inverse backstepping transformation to derive the reference trajectory and reference input for the original system. Our approach is different than the one in [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2 × 2 linear hyperbolic system[END_REF], in that we use backstepping for trajectory generation rather than stabilization, and the one in [START_REF] Petit | Dynamics and solutions to some control problems for water-tank systems[END_REF], in that we employ a different conceptual idea to a different class of systems. The idea of the backstepping-based trajectory generation for PDEs, which was conceived in [START_REF] Krstic | Boundary Control of PDEs[END_REF], is applied to a beam PDE in [START_REF] Siranosian | Motion planning and tracking for tip displacement and deflection angle for flexible beams[END_REF] and the Navier-Stokes equations in [START_REF] Cochran | Motion planning and trajectory tracking for the 3-D Poiseuille flow[END_REF], and is recently extended to general n × n linear hyperbolic systems in [START_REF] Hu | Control of homodirectional and general heterodirectional linear coupled hyperbolic PDEs[END_REF]. We apply this methodology to a wave PDE with indefinite in-domain and boundary damping by transforming (see, for example, [START_REF] Bekiaris-Liberis | Compensation of wave actuator dynamics for nonlinear systems[END_REF]) the wave PDE to a 2 × 2 linear hyperbolic system coupled with a first-order ODE (Section 2.2).

We then employ a PI control law for the stabilization of the error system, namely, the system whose state is defined as the difference between the state of the plant and the reference trajectory. We prove exponential stability in the L 2 norm of the closed-loop system by constructing a Lyapunov functional which incorporates cross-terms between the PDE states of the system and the ODE state of the controller, when the parameters of the system and the controller satisfy certain conditions (Section 3.1). Our result differs than the result in [START_REF] Santos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF] in that we employ PI control on an output of the system in the Riemann coordinates and we construct a non-diagonal Lyapunov functional for proving closed-loop stability. We demonstrate that the proposed PI control law is capable of compensating in the output the effect of additive disturbances affecting the boundary or the interior of the PDE domain (Section 3.2). We present several examples, for the illustration of our methodologies, including a simulation example dealing with the generation of a sinusoidal reference trajectory for a wave PDE (Section 4.1) and a simulation example of a system tracking a sinusoidal reference output (Section 4.2).

Trajectory Generation Using Backstepping

General 2 × 2 Linear Hyperbolic Systems

We consider the following system

z 1 t + ε 1 (x)z 1 x = c 1 (x)z 1 + c 2 (x)z 2 (1) z 2 t -ε 2 (x)z 2 x = c 3 (x)z 1 + c 4 (x)z 2 , (2)
under the boundary conditions

z 1 (0,t) = qz 2 (0,t) (3) z 2 (1,t) = S(t) (4) z 2 (0,t) = y(t) , (5)
where t ∈ [0, +∞) is the time variable, x ∈ [0, 1] is the spatial variable, y is the output of the system, and S is the control input. The functions ε 1 , ε 2 belong to C 2 ([0, 1]) and satisfy ε 1 (x), ε 2 (x) > 0, for all x ∈ [0, 1], and the functions c i , i = 1, 2, 3, 4 belong to C 1 ([0, 1]). Defining the change of variables (see, for example, [START_REF] Bastin | On boundary feedback stabilization of non-uniform linear 2 × 2 hyperbolic systems over a bounded interval[END_REF])

χ 1 (x) = exp - x 0 c 1 (s) ε 1 (s) ds (6)
χ 2 (x) = exp x 0 c 4 (s) ε 2 (s) ds (7)
χ(x) = χ 1 (x) χ 2 (x) , (8)
and the new coordinates

u = χ 1 (x)z 1 (9) v = χ 2 (x)z 2 , (10)
system ( 1)-( 5) is transformed into the following system

u t + ε 1 (x)u x = γ 1 (x)v (11) v t -ε 2 (x)v x = γ 2 (x)u , (12)
with

γ 1 (x) = χ(x)c 2 (x) (13)
γ 2 (x) = χ -1 (x)c 3 (x) . (14)
The boundary conditions become

u(0,t) = qv(0,t) (15) v(1,t) = U(t) (16) v(0,t) = y(t) , (17)
where the original control variable satisfies

U = χ 2 (1)S . (18)
We aim at designing a reference control input U r (t) such that the output y(t) follows a given reference trajectory y r (t), for t ≥ 0. For achieving this we need first to construct the reference trajectory (u r (x,t), v r (x,t)) that satisfies [START_REF] Di Meglio | Backstepping stabilization of an underactuated 3 × 3 linear hyperbolic system of fluid flow transport equations[END_REF], [START_REF] Di Meglio | Stabilization of a system of n + 1 coupled firstorder hyperbolic linear PDEs with a single boundary input[END_REF], [START_REF] Dubois | Motion planning and nonlinear simulations for a tank containing a fluid[END_REF], and [START_REF] Fliess | A distributed parameter approach to the control of a tubular reactor: A multi-variable case[END_REF] with y(t) = y r (t). The trajectory generation problem is solvable when the initial data (u 0 , v 0 ) match the reference trajectory, i.e., when u 0 (x) = u r (x, 0) and v 0 (x) = v r (x, 0) (and hence, the initial data belong to the same space with u r (x, 0) and v r (x, 0)).

Theorem 1. Let y r ∈ C 1 (R) be uniformly bounded. The functions u r (x,t) = qy r (t -Φ 1 (x)) +

x 0 f (ξ ) ε 1 (ξ ) y r (t -Φ 1 (x) + Φ 1 (ξ )) dξ +q x 0 L αα (x, ξ )y r (t -Φ 1 (ξ ))dξ + x 0 L αα (x, ξ ) ξ 0 f (ζ ) ε 1 (ζ ) y r (t -Φ 1 (ξ ) + Φ 1 (ζ )) dζ dξ + x 0 L αβ (x, ξ ) y r (t + Φ 2 (ξ )) dξ (19) v r (x,t) = y r (t + Φ 2 (x)) + q x 0 L β α (x, ξ )y r (t -Φ 1 (ξ )) dξ + x 0 L β α (x, ξ ) ξ 0 f (ζ ) ε 1 (ζ ) y r (t -Φ 1 (ξ ) + Φ 1 (ζ )) dζ dξ + x 0 L β β (x, ξ )y r (t + Φ 2 (ξ )) dξ (20)
U r (t) = y r (t + Φ 2 (1)) + q

1 0 L β α (1, ξ )y r (t -Φ 1 (ξ )) dξ + 1 0 L β α (1, ξ ) ξ 0 f (ζ ) ε 1 (ζ ) y r (t -Φ 1 (ξ ) + Φ 1 (ζ )) dζ dξ + 1 0 L β β (1, ξ )y r (t + Φ 2 (ξ )) dξ , (21)
where

Φ 1 (x) = x 0 1 ε 1 (s) ds (22)
Φ 2 (x) = x 0 1 ε 2 (s) ds (23)
f (x) =    ε 2 (0)K uv (x, 0) , if q = 0 0 , if q = 0 , (24)
and L αα , L αβ , L β α , L β β , K uv are the solutions of the following equations

ε 2 (x)L β α x -ε 1 (ξ )L β α ξ = ε 1 (ξ )L β α -γ 2 (x)L αα (25) ε 2 (x)L β β x + ε 2 (ξ )L β β ξ = -ε 2 (ξ )L β β -γ 2 (x)L αβ (26) ε 1 (x)L αα x + ε 1 (ξ )L αα ξ = -ε 1 (ξ )L αα + γ 1 (x)L β α (27) ε 1 (x)L αβ x -ε 2 (ξ )L αβ ξ = ε 2 (ξ )L αβ + γ 1 (x)L β β (28) ε 1 (x)K uu x + ε 1 (ξ )K uu ξ = -ε 1 (ξ )K uu -γ 2 (x)K uv (29) ε 1 (x)K uv x -ε 2 (ξ )K uv ξ = ε 2 (ξ )K uv -γ 1 (x)K uu , (30)
with the boundary conditions

L β α (x, x) = - γ 2 (x) ε 1 (x) + ε 2 (x) (31)
L αα (x, 0) =    h 1 (x) , if q = 0 ε 2 (0) qε 1 (0) L αβ (x, 0) , if q = 0 ( 32
)
L β β (x, 0) =    1 ε 2 (0) x 0 L β α (x, ξ ) f (ξ )dξ , if q = 0 qε 1 (0) ε 2 (0) L β α (x, 0) , if q = 0 (33) L αβ (x, x) = γ 1 (x) ε 1 (x) + ε 2 (x) (34)
K uu (x, 0) = h 2 (x) (35)
K uv (x, x) = γ 1 (x) ε 1 (x) + ε 2 (x) , (36)
where h 1 , h 2 ∈ C 1 ([0, 1]) are arbitrary, are uniformly bounded and solve the boundary value problem [START_REF] Di Meglio | Backstepping stabilization of an underactuated 3 × 3 linear hyperbolic system of fluid flow transport equations[END_REF], [START_REF] Di Meglio | Stabilization of a system of n + 1 coupled firstorder hyperbolic linear PDEs with a single boundary input[END_REF], [START_REF] Dubois | Motion planning and nonlinear simulations for a tank containing a fluid[END_REF], [START_REF] Fan | Comparative model accuracy of a data-fitted generalized Aw-Rascle-Zhang model[END_REF]. In particular, v r (0,t) = y r (t), for t ≥ 0.

Before proving Theorem 1 we make the following observation, which is also helpful in understanding better the proof strategy of Theorem 1.

Remark 1. The approach for the trajectory generation introduced here is inspired from backstepping. Consider the following system

α t + ε 1 (x)α x -f (x)β (0,t) = 0 ( 37
)
β t -ε 2 (x)β x = 0 , (38)
with boundary condition

α(0,t) = qβ (0,t) , (39)
which follows by directly applying the backstepping transformation

α(x,t) = u(x,t) - x 0 K uu (x, ξ ) u (ξ ,t) dξ - x 0 K uv (x, ξ ) v (ξ ,t) dξ (40) β (x,t) = v(x,t) - x 0 K vu (x, ξ ) u (ξ ,t) dξ - x 0 K vv (x, ξ ) v (ξ ,t) dξ , (41)
where the kernels K uu , K uv , K vu , K vv are given in [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2 × 2 linear hyperbolic system[END_REF], to system (11), [START_REF] Di Meglio | Stabilization of a system of n + 1 coupled firstorder hyperbolic linear PDEs with a single boundary input[END_REF], and [START_REF] Dubois | Motion planning and nonlinear simulations for a tank containing a fluid[END_REF]. It is shown that the functions

α(x,t) = qy r (t -Φ 1 (x)) + x 0 f (ξ ) ε 1 (ξ ) y r (t -Φ 1 (x) + Φ 1 (ξ )) dξ (42) β (x,t) = y r (t + Φ 2 (x)) , (43)
where Φ 1 and Φ 2 are defined in ( 22) and ( 23), respectively, satisfy (37)-(39) with

β (1,t) = y r (t + Φ 2 (1)) (44)
and, in particular, β (0,t) = y r (t). Using the inverse backstepping transformations introduced in [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2 × 2 linear hyperbolic system[END_REF] u(x,t) = α(x,t)

+ x 0 L αα (x, ξ )α(ξ ,t)dξ + x 0 L αβ (x, ξ )β (ξ ,t)dξ (45) v(x,t) = β (x,t) + x 0 L β α (x, ξ )α(ξ ,t)dξ + x 0 L β β (x, ξ )β (ξ ,t)dξ , (46)
and relations (42), (43) one can conclude that the functions u r , v r , and U r = v r (1) solve the trajectory generation problem for system [START_REF] Di Meglio | Backstepping stabilization of an underactuated 3 × 3 linear hyperbolic system of fluid flow transport equations[END_REF], ( 12), ( 15)- [START_REF] Fliess | A distributed parameter approach to the control of a tubular reactor: A multi-variable case[END_REF].

Note that the present approach cannot be directly applied to cases where ε 1 (x) or ε 2 (x) vanish for some x ∈ [0, 1]. This is evident, for instance, from [START_REF] Vazquez | Marcum Q-functions and explicit kernels for stabilization of linear hyperbolic systems with constant coefficients[END_REF] which would imply that the kernel

L β β
of the open-loop control law U r may become infinity for all x ∈ [0, 1].

Proof We first consider the case q = 0. Note that since ε 1 ,

ε 2 ∈ C 2 ([0, 1]) with ε 1 (x), ε 2 (x) > 0, for all x ∈ [0, 1] and γ 1 , γ 2 ∈ C 1 ([0, 1]), system (25)-(34) has a unique solution with L αα , L αβ , L β α , L β β ∈ C 1 (T ) where T = {(x, ξ ) : 0 ≤ ξ ≤ x ≤ 1} [7]
. Hence, from ( 19)-( 21) and the uniform boundedness of y r it follows that u r , v r , and U r are bounded for all t ≥ 0 and x ∈ [0, 1].

Taking the time and space derivatives of u r we get

u r t + ε 1 (x)u r x = q x 0 L αα (x, ξ )y r (t -Φ 1 (ξ ))dξ + x 0 L αβ (x, ξ )y r (t + Φ 2 (ξ ))dξ + ε 1 (x) x 0 L αβ x (x, ξ )y r (t + Φ 2 (ξ )) dξ + qε 1 (x) x 0 L αα x (x, ξ )y r (t -Φ 1 (ξ )) dξ + ε 1 (x)L αβ (x, x)y r (t + Φ 2 (x)) + qε 1 (x)L αα (x, x)y r (t -Φ 1 (x)) . ( 47
)
Integrating by parts the first two integrals we get

u r t + ε 1 (x)u r x = q x 0 ε 1 (x)L αα x (x, ξ ) + ε 1 (ξ )L αα ξ (x, ξ ) + ε 1 (ξ )L αα (x, ξ ) y r (t -Φ 1 (ξ ))dξ + x 0 ε 1 (x)L αβ x (x, ξ ) -ε 2 (ξ ) L αβ ξ (x, ξ ) -ε 2 (ξ ) L αβ (x, ξ ) y r (t + Φ 2 (ξ ))dξ + qε 1 (0)L αα (x, 0) -ε 2 (0)L αβ (x, 0) y r (t) + (ε 1 (x) + ε 2 (x)) L αβ (x, x)y r (t + Φ 2 (x)) . (48)
Due to the fact that L αβ and L αα are the solutions of ( 27) and ( 28) with the boundary conditions [START_REF] Tang | Sliding mode control to stabilization of linear 2 × 2 hyperbolic systems with boundary input disturbance[END_REF] and ( 34) one gets, by using [START_REF] Gugat | Existence of classical solutions and feedback stabilization for the flow in gas networks[END_REF], that u r satisfies [START_REF] Di Meglio | Backstepping stabilization of an underactuated 3 × 3 linear hyperbolic system of fluid flow transport equations[END_REF]. The proof that v r satisfies (12) follows analogously. Setting x = 0 in [START_REF] Gugat | Gas flow in fan-shaped networks: Classical solutions and feedback stabilization[END_REF], [START_REF] Gugat | Existence of classical solutions and feedback stabilization for the flow in gas networks[END_REF] and using ( 22), ( 23), we get that u r and v r satisfy [START_REF] Dubois | Motion planning and nonlinear simulations for a tank containing a fluid[END_REF]. Setting x = 1 in (20) it follows that ( 21) satisfies [START_REF] Fan | Comparative model accuracy of a data-fitted generalized Aw-Rascle-Zhang model[END_REF]. Setting in [START_REF] Gugat | Existence of classical solutions and feedback stabilization for the flow in gas networks[END_REF] x = 0 and using [START_REF] Krstic | Boundary Control of PDEs[END_REF] we get v r (0,t) = y r (t).

Let us consider next the case q = 0. First observe that the PDEs ( 25), [START_REF] Petit | Flatness of heavy chain systems[END_REF] with boundary conditions ( 31), [START_REF] Tang | Sliding mode control to stabilization of linear 2 × 2 hyperbolic systems with boundary input disturbance[END_REF], for the kernels L αα and L β α are decoupled, and hence, L αα and L β α are well-defined [START_REF] Coron | Local exponential H 2 stabilization of a 2 × 2 quasilinear hyperbolic system using backstepping[END_REF]. Hence, since f satisfies (24) and K uv , K uu are well-defined [START_REF] Coron | Local exponential H 2 stabilization of a 2 × 2 quasilinear hyperbolic system using backstepping[END_REF], one can conclude that L αβ and L β β are well-defined as well.

Taking the time and space derivatives of u r we get

u r t + ε 1 (x)u r x = f (x)y r (t) + x 0 L αα (x, ξ ) ξ 0 f (ζ ) ε 1 (ζ ) y r (t -Φ 1 (ξ ) + Φ 1 (ζ )) dζ dξ + x 0 L αβ (x, ξ )y r (t + Φ 2 (ξ ))dξ +ε 1 (x)L αα (x, x) x 0 f (ζ ) ε 1 (ζ ) y r (t -Φ 1 (x) + Φ 1 (ζ ))dζ +ε 1 (x)L αβ (x, x)y r (t + Φ 2 (x)) + ε 1 (x) x 0 L αβ x (x, ξ )y r (t + Φ 2 (ξ ))dξ +ε 1 (x) x 0 L αα x (x, ξ ) ξ 0 f (ζ ) ε 1 (ζ ) y r (t -Φ 1 (ξ ) + Φ 1 (ζ ))dζ dξ . ( 49
)
Integrating by parts the first two integrals we get

u r t + ε 1 (x)u r x = x 0 ε 1 (x)L αα x (x, ξ ) + ε 1 (ξ )L αα ξ (x, ξ ) + ε 1 (ξ )L αα (x, ξ ) × ξ 0 f (ζ ) ε 1 (ζ ) y r (t -Φ 1 (ξ ) + Φ 1 (ζ ))dζ dξ + x 0 ε 1 (x)L αβ x (x, ξ ) -ε 2 (ξ )L αβ ξ (x, ξ ) -ε 2 (ξ )L αβ (x, ξ ) y r (t + Φ 2 (ξ ))dξ + (ε 1 (x) + ε 2 (x)) L αβ (x, x)y r (t + Φ 2 (x)) +y r (t) f (x) + x 0 L αα (x, ξ ) f (ξ )dξ -ε 2 (0)L αβ (x, 0) . (50)
Using ( 20), ( 27), [START_REF] Prieur | Control of systems of conservation laws with boundary errors[END_REF], and (34) one can conclude that u r satisfies (11) if f satisfies

f (x) = ε 2 (0)L αβ (x, 0) - x 0 L αα (x, ξ ) f (ξ )dξ . (51)
This fact can been shown as follows. The inverse of the backstepping transformation (40), ( 41) is uniquely defined and has the form (45), (46) (see, for example, [START_REF] Hochstadt | Integral Equations[END_REF]). Hence, substituting (40), ( 41) in ( 45), (46) we get

x 0 (K uu (x, ξ ) -L αα (x, ξ )) u(ξ ,t) + K uv (x, ξ ) -L αβ (x, ξ ) v(ξ ,t)dξ + x 0 ξ 0 L αα (x, ξ )K uu (ξ , ζ ) + L αβ (x, ξ )K vu (ξ , ζ ) u(ζ ,t) + L αα (x, ξ )K uv (ξ , ζ ) + L αβ (x, ξ )K vv (ξ , ζ ) v(ζ ,t) dζ dξ = 0 . ( 52
)
Performing a change in the order of integration in the second integral of (52) and using the fact that (52) holds for all u and v, one obtains

K uv (x, ξ ) = L αβ (x, ξ ) - x ξ L αα (x, s)K uv (s, ξ ) + L αβ (x, s)K vv (s, ξ ) ds . ( 53
)
Setting ξ = 0 in (53), multiplying (53) by ε 2 (0), and using the facts that K vv (x, 0) = 0 for all x ∈ [0, 1] (see relation [START_REF] Siranosian | Motion planning and tracking for tip displacement and deflection angle for flexible beams[END_REF] in [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2 × 2 linear hyperbolic system[END_REF]) and that f is defined by [START_REF] Li | Classical solutions for quasilinear hyperbolic systems[END_REF], we get that f satisfies (51) for q = 0. The rest of the proof is similar to the case q = 0.

Example 1. We consider the following system

z 1 t + ε 1 z 1 x = - 1 τ z 1 (54) z 2 t -ε 2 z 2 x = - 1 τ z 1 , (55)
with boundary conditions

z 1 (0,t) = qz 2 (0,t) (56)
z 2 (1,t) = S(t) , ( 57
)
where τ is a positive parameter. Among various systems that can be modeled by ( 54)-(57) (for instance, the Saint-Venant equations, see [START_REF] Santos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF], [START_REF] De Halleux | Boundary feedback control in networks of open channels[END_REF]), system (54)-(57) can be viewed as a linearized version of the Aw-Rascle-Zhang (ARZ) macroscopic model of traffic flow in the Riemann coordi-nates

z 1 = w -V (s * )s (58)
z 2 = w , (59)
where w and s correspond to the velocity and density of the vehicles at time t and location x, respectively. The variable V (s * ) is the nominal velocity of the cars and s * is the nominal density.

The opposite transport velocities in (54), (55) correspond to traffic flow in a congested mode. The parameter 1 τ is an indicator of the convergence rate of the velocity w of the cars to the nominal velocity V (s). For more details the reader is referred to [START_REF] Fan | Comparative model accuracy of a data-fitted generalized Aw-Rascle-Zhang model[END_REF]. The boundary condition (56) in the original variables is written as

w = V (s * )s 1 -q . ( 60
)
Hence, the boundary condition (56) dictates that there is a static relation, at the entrance of the road, between the density and the velocity similarly to the static relation between the nominal velocity V (s) and the density of the cars in the road. The change of variables ( 9), ( 10), [START_REF] Santos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF], and (

transform system (54)-(57) to

u t + ε 1 u x = 0 (61) v t -ε 2 v x = - 1 τ exp - 1 τε 1 x u (62) u(0,t) = qv(0,t) (63) v(1,t) = U(t) , (64)
where U(t) is given by [START_REF] Goatin | The Aw-Rascle vehicular traffic flow model with phase transitions[END_REF]. Observing that γ 1 = 0, relations ( 25)-( 34) can be solved explicitly as

L αα (x, ξ ) = 0 (65) L αβ (x, ξ ) = 0 (66) L β α (x, ξ ) = 1 τ (ε 1 + ε 2 ) exp - 1 τε 1 ε 1 x + ε 2 ξ ε 1 + ε 2 (67) L β β (x, ξ ) = qε 1 τε 2 (ε 1 + ε 2 ) exp - 1 τε 1 ε 1 x -ε 1 ξ ε 1 + ε 2 . ( 68
)
Therefore, for system (54)-(57), the reference input which generates the desired output y r (t) is S r (t) = y r t + 1

ε 2 + q τ (ε 1 + ε 2 ) 1 0 exp - 1 τε 1 ε 1 + ε 2 ξ ε 1 + ε 2 y r t - ξ ε 1 dξ + qε 1 τε 2 (ε 1 + ε 2 ) 1 0 exp - 1 τε 1 ε 1 -ε 1 ξ ε 1 + ε 2 y r t + ξ ε 2 dξ . ( 69
)
2.2. Application to a Wave PDE with Indefinite In-Domain and Boundary Damping Let us consider system

z tt = ε(x)z xx + h(x)z t + b(x)z x (70) z x (0,t) = -gz t (0,t) (71) z x (1,t) = W (t) , (72)
with g = 1 √ ε(0) , -1 √ ε(0) , h, b ∈ C 1 ([0, 1]), and ε ∈ C 2 ([0, 1]) with ε(x) > 0, for all x ∈ [0, 1].
The objective is z(0,t) to track a reference trajectory, say, ζ (t), which belongs to C 2 (R). Let us define the output of the system as

ψ(t) = z(0,t) . (73)
With the change of variables

z 1 (x,t) = 1 -ε(0)g 1 + ε(0)g z t (x,t) -ε(x)z x (x,t) (74)
z 2 (x,t) = z t (x,t) + ε(x)z x (x,t) (75)
S(t) = ε(1)W (t) + z t (1,t) , (76)
system (70)-( 72) is rewritten as (1)- [START_REF] Coron | A strict Lyapunov function for boundary control of hyperbolic systems of conservation laws[END_REF] where

y(t) = 1 -ε(0)g ψ(t) (77) ε 1 (x) = ε(x) (78) ε 2 (x) = ε(x) (79) q = 1 (80) c 1 (x) = h(x) 2 - b(x) 2 ε(x) + ε (x) 4 ε(x) (81) c 2 (x) = mc 4 (x) (82) c 3 (x) = 1 m c 1 (x) (83) c 4 (x) = h(x) 2 + b(x) 2 ε(x) - ε (x) 4 ε(x) (84) m = 1 -ε(0)g 1 + ε(0)g , (85)
together with the integrator ψ(t) = 1 1-√ ε(0)g z 2 (0,t). Applying Theorem 1 we get the following reference input

W r (t) = 1 2 ε(1) 1 -ε(0)g exp - 1 0 c 4 (s) ε 2 (s) ds ζ (t + Φ 2 (1)) + 1 0 L β α (1, ξ ) ζ (t -Φ 1 (ξ )) dξ + 1 0 L β β (1, ξ ) ζ (t + Φ 2 (ξ )) dξ -1 + ε(0)g exp 1 0 c 1 (s) ε 1 (s) ds ζ (t -Φ 1 (1)) + 1 0 L αα (x, ξ ) ζ (t -Φ 1 (ξ )) dξ + 1 0 L αβ (x, ξ ) ζ (t + Φ 2 (ξ )) dξ . ( 86
)

Trajectory Tracking Using PI Control

Stability Analysis with a Non-Diagonal Lyapunov Functional

For stabilizing the system around the desired trajectory for any initial condition (u(x, 0), v(x, 0)), rather than only for (u(x, 0), v(x, 0)) = (u r (x, 0), v r (x, 0)), we employ a PI-feedback control law. We first write the dynamics of the tracking errors ũ(x,t) = u(x,t)u r (x,t) and ṽ(x,t) = v(x,t)v r (x,t)

as ũt + ε 1 (x) ũx = γ 1 (x) ṽ (87) ṽt -ε 2 (x) ṽx = γ 2 (x) ũ (88) ũ(0,t) = q ṽ(0,t) (89) ṽ(1,t) = Ũ(t) , (90)
where Ũ = U -U r and U r is the reference input generating the desired reference trajectory. We employ the controller Ũ(t) = -k P ṽ(0,t)

-k I η(t) , (91)
with η(t) = ṽ(0,t) .

Theorem 2. Consider system (87)-(90) together with the control law (91), (92). Let the positive constants µ, β , ρ, γ, ν, κ, and θ be such that the matrices

M 1 =   -q 2 -β k 2 P e µ -1 -κγ 2 -β k P k I e µ + γ 2 (e ν k P + 1) -ρ 2 -β k P k I e µ + γ 2 (e ν k P + 1) -ρ 2 -β k 2 I e µ + γe ν k I -γ 2   (93) M 2 (x) =   M 21 (x) M 22 (x) M 23 (x) M 24 (x)   (94) with M 21 (x) = µ - θ ε 1 (x) e -µx + γ 2 2 (θ ρ -γ) γ 2 2 (x) ε 2 2 (x)
e 2νx (95)

M 22 (x) = - γ 1 (x) ε 1 (x) e -µx -β γ 2 (x) ε 2 (x) e µx - γ 2 2 (θ ρ -γ) γ 2 (x) ε 2 (x) ν - θ ε 2 (x)
e 2νx (96)

M 23 (x) = - γ 1 (x) ε 1 (x) e -µx -β γ 2 (x) ε 2 (x) e µx - γ 2 2 (θ ρ -γ) γ 2 (x) ε 2 (x) ν - θ ε 2 (x)
e 2νx (97)

M 24 (x) = β µ - θ ε 2 (x) e µx - γ 2κ
e 2νx ε 2 2 (x) + γ 2 2 (θ ρ -γ) ν - θ ε 2 (x) 2 e 2νx , ( 98
)
are positive semi-definite for all x ∈ [0, 1], and the inequalities

β ρ > γ 2 e (2ν-µ)x 2ε 2 (x) , ∀x ∈ [0, 1] (99) γ > θ ρ , (100)
hold. Then, there exist positive constants λ and Ω such that, for all initial conditions satisfying ũ0 (x), ṽ0 (x), η0 ∈ L 2 (0, 1) × L 2 (0, 1) × R, the following holds for all t ≥ 0

1 0 ũ2 (x,t) + ṽ2 (x,t) dx + η2 (t) ≤ Ωe -λt 1 0 ũ2 (x, 0) + ṽ2 (x, 0) dx + η2 (0) . ( 101
)
Proof In order to analyze the stability of system (87)-(92) we propose the following Lyapunov functional

V (t) = 1 0   ũ(x,t) ṽ(x,t) η(t)   P(x)   ũ(x,t) ṽ(x,t) η(t)   dx = R 1 (t) + R 2 (t) + R 3 (t) + R 4 (t) , (102)
with

P(x) =     e -µx ε 1 (x) 0 0 0 β e µx ε 2 (x) γe νx 2ε 2 (x) 0 γe νx 2ε 2 (x) ρ 2     , (103)
and

R 1 (t) = 1 0 ũ2 (x,t) e -µx ε 1 (x) dx (104) R 2 (t) = β 1 0 ṽ2 (x,t) e µx ε 2 (x) dx (105) R 3 (t) = γ η(t) 1 0 ṽ(x,t) e νx ε 2 (x) dx (106) R 4 (t) = ρ 2 η2 (t) . ( 107
)
Let us introduce the constants

λ = min x∈[0,1] λ min (P(x)) (108) λ = max x∈[0,1] λ max (P(x)) . ( 109
)
Inequality (99) ensures that P(x) is positive definite and symmetric for all x ∈ [0, 1], and hence, using the fact that

ε 1 , ε 2 ∈ C 2 ([0, 1]) with ε 1 (x), ε 2 (x) > 0, for all x ∈ [0, 1], one can conclude that, λ , λ > 0. Therefore, λ 1 0 ũ2 (x,t) + ṽ2 (x,t) dx + η2 (t) ≤ V (t) ≤ λ 1 0 ũ2 (x,t) + ṽ2 (x,t) dx + η2 (t) . (110)
Using ( 104)-(107) we get along the solutions of system (87)-(92) that

Ṙ1 (t) = -2 1 0 ũ(x,t) ũx (x,t)e -µx dx + 2 1 0 ũ(x,t) ṽ(x,t) γ 1 (x) ε 1 (x) e -µx dx = q 2 ṽ2 (0,t) -e -µ ũ2 (1,t) -µ 1 0 ũ2 (x,t)e -µx dx + 2 1 0 ũ(x,t) ṽ(x,t) γ 1 (x) ε 1 (x) e -µx dx (111) Ṙ2 (t) = 2β 1 0 ṽ(x,t) ṽx (x,t)e µx dx + 2β 1 0 ũ(x,t) ṽ(x,t) γ 2 (x) ε 2 (x) e µx dx = β k 2 P e µ ṽ2 (0,t) + 2k P k I e µ ṽ(0,t) η(t) + k 2 I e µ η2 (t) -ṽ2 (0,t) -µβ 1 0 ṽ2 (x,t)e µx dx + 2β 1 0 ũ(x,t) ṽ(x,t) γ 2 (x) ε 2 (x) e µx dx (112) Ṙ3 (t) = γ η(t) 1 0 ṽx (x,t)e νx dx + γ ṽ(0,t) 1 0 ṽ(x,t) e νx ε 2 (x) dx + γ η(t) 1 0 ũ(x,t) γ 2 (x) ε 2 (x) e νx dx ≤ γ η(t) (e ν (-k P ṽ(0,t) -k I η(t)) -ṽ(0,t)) -νγ η(t) 1 0 ṽ(x,t)e νx dx + κγ 2 ṽ2 (0,t) + γ 2κ 1 0 ṽ2 (x,t) e 2νx ε 2 2 (x) dx + γ η(t) 1 0 ũ(x,t) γ 2 (x) ε 2 (x) e νx dx (113)
Ṙ4 (t) = ρ ṽ(0,t) η(t) , (114)
where we used integration by parts in the first terms of (111)-( 113) and Young's inequality in the second term of (113). Using (102), ( 111)-(114) we get

V (t) ≤ - ṽ(0,t) η(t) M 1 ṽ(0,t) η(t) - 1 0   ũ(x,t) ṽ(x,t) η(t)   M(x)   ũ(x,t) ṽ(x,t) η(t)   dx -e -µ ũ2 (1,t) -θV (t) , (115)
where M 1 is given by (93) and

M(x) = A(x) B (x) B(x) C , (116)
with

A(x) = A 1 (x) A 2 (x) A 3 (x) A 4 (x) , (117)
where

A 1 (x) = µ - θ ε 1 (x) e -µx (118)
A 2 (x) = - γ 1 (x) ε 1 (x) e -µx -β γ 2 (x) ε 2 (x) e µx (119)
A 3 (x) = - γ 1 (x) ε 1 (x) e -µx -β γ 2 (x) ε 2 (x) e µx (120)
A 4 (x) = β µ - θ ε 2 (x) e µx - γ 2κ
e 2νx ε 2 2 (x) (121)
B(x) = -γ 2 γ 2 (x) ε 2 (x) e νx γ 2 ν -θ ε 2 (x) e νx (122)
C = γ -θ ρ 2 . (123)
Using the Schur complement of C in M(x) and (100), (123) one has that M(x) ≥ 0 for all x ∈ [0, 1], if and only if

M 2 (x) = A(x) -B (x)C -1 B(x) ≥ 0 . (124)
Thus, if M 1 ≥ 0 and M 2 (x) ≥ 0, for all x ∈ [0, 1], one has

V (t) ≤ -e -µ ũ2 (1,t) -θV (t) , (125)
and hence, V (t) ≤ e -θt V (0), for all t ≥ 0. Combining this relation with (110) the proof is complete.

Remark 2. A control law with an integral action is designed in [START_REF] Santos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF] for 2 × 2 hyperbolic systems.

Stability of the closed-loop system is proved using a diagonal Lyapunov functional. Here the nondiagonal term in the Lyapunov functional is needed for proving stability using a quadratic Lyapunov functional. Indeed, let us assume that the Lyapunov functional is diagonal. We can write it as

V (t) = 1 0 q 1 (x) ũ2 (x,t) + q 2 (x) ṽ2 (x,t) dx + ρ 2 η2 (t) , (126)
where the functions q 1 and q 2 belong to C 1 ([0, 1]) with q 1 (x), q 2 (x) > 0, for all x ∈ [0, 1]. The time derivative of V along the solutions of system (87), (88) with boundary conditions (89)-( 92) is given by

V (t) =   ṽ(0,t) η(t)     D 1 D 2 D 3 D 4     ṽ(0,t) η(t)   + 1 0   ũ(x,t) ṽ(x,t)   E(x)   ũ(x,t) ṽ(x,t)   dx -q 1 (1)ε 1 (1) ũ2 (1,t) , (127)
where

D 1 = q 1 (0)ε 1 (0)q 2 -q 2 (0)ε 2 (0) + q 2 (1)ε 2 (1)k 2 P (128)
D 2 = 1 2 (q 2 (1)ε 2 (1)k P k I + ρ) (129)
D 3 = 1 2 (q 2 (1)ε 2 (1)k P k I + ρ) (130)
D 4 = q 2 (1)ε 2 (1)k 2 I (131) E(x) =   (q 1 (x)ε 1 (x)) x q 1 (x)γ 1 (x) + q 2 (x)γ 2 (x) q 1 (x)γ 1 (x) + q 2 (x)γ 2 (x) -(q 2 (x)ε 2 (x)) x   . (132)
Using ( 127) and (131) one can conclude that when k I = 0 the inequality V ≤ 0 can not be satisfied for any ũ ṽ η .

In [START_REF] Bastin | On boundary feedback stabilization of non-uniform linear 2 × 2 hyperbolic systems over a bounded interval[END_REF], it is proved that if there exist two boundary controllers for 2×2 linear hyperbolic systems of the form (87), (88) such that the functional

V (t) = 1 0 q 1 (x) ũ2 (x,t) + q 2 (x) ṽ2 (x,t) + q 3 (x) ũ(x,t) ṽ(x,t)dx , (133)
along the solutions of the system (87), (88) with the state ( ũ, ṽ) satisfies V < 0 then the cross term q 3 between ũ and ṽ is necessarily identically zero. However, in the case of stabilization of 2 × 2 linear hyperbolic systems of the form (87), (88) with a PI control law that we consider here, the cross term (106) in the Lyapunov functional (102) between the integral state η of the controller and the state of the plant ṽ is necessary (as explained above) for proving stability of the overall closed-loop system consisting of the plant state ( ũ, ṽ) and the integral state η, using the Lyapunov functional defined in (102) (although a cross term between ũ and ṽ is not necessary).

As explained in Remark 2 the non-diagonal term in the Lyapunov functional is crucial for proving stability using a quadratic Lyapunov functional. However, this term adds considerable complexity in verifying analytically that the matrices (93), (94) are positive definite and that (99) holds. In general, the positivity of M 1 and M 2 must be verified numerically. Yet, from the expression of M 1 we see that a necessary condition is that k I is strictly positive. In addition, from (93) it is evident that k P must satisfy |k P | < 1. Note that from (93)-(100) it seems possible that the positivity of M 1 and M 2 may depend on the values of γ 1 , γ 2 , and q. Next, we numerically verify the conditions of Theorem 2 for the system from Example 1.

Example 2 (Example 1 Continued). We set in (61)-( 63)

ε 1 = 3 (134) ε 2 = 6 (135) τ = 5 (136)
q = 0.2 , (137)
and choose U in (90) according to (91) with

k P = 0.1 (138) k I = 1.0583 , (139)
in order to stabilize the zero equilibrium of (61)-(63). We verify numerically that the conditions of Theorem 2 are satisfied with (β , κ, µ, ν, θ , ρ, γ) = (0.7, 0.2, 0.5, 0.2, 0.7, 2, 2) . The verification of the positive definiteness of matrix ( 94) is more delicate due to its dependence on

x. Figure 1 shows the evolution of the eigenvalues of M 2 (x) and the determinant of matrix (103), which remain positive for all x ∈ [0, 1].

Compensation in the Output of In-Domain and Boundary Disturbances

Let us assume that there exist some disturbances d 1 , d 2 ∈ C 1 ([0, 1]) on the right-hand side of (11), [START_REF] Di Meglio | Stabilization of a system of n + 1 coupled firstorder hyperbolic linear PDEs with a single boundary input[END_REF], respectively and some disturbances d 3 , d 4 ∈ R on the right-hand side of ( 15), [START_REF] Fan | Comparative model accuracy of a data-fitted generalized Aw-Rascle-Zhang model[END_REF], respectively. The error system (87)-(90) becomes

ũt + ε 1 (x) ũx = γ 1 (x) ṽ + d 1 (x) (142) ṽt -ε 2 (x) ṽx = γ 2 (x) ũ + d 2 (x) (143) ũ(0,t) = q ṽ(0,t) + d 3 (144) ṽ(1,t) = Ũ(t) + d 4 , (145)
with Ũ(t) and η(t) given by ( 91) and (92) respectively. The equilibrium of the perturbed system (142)-( 144) and ( 90) is the solution of the following ordinary differential equation

Z (x) = F(x)Z(x) + G(x), (146)
where

F(x) = 0 γ 1 (x) ε 1 (x) - γ 2 (x) ε 2 (x) 0 and G(x) = d 1 (x) ε 1 (x) - d 2 (x) ε 2 (x)
, with boundary conditions

Z 1 (0) = d 3 (147) Z 2 (0) = 0 . (148)
The ordinary differential equation ( 146

Let us define

u(x,t) = ũ(x,t) -ũss (x; d 1 , d 2 , d 3 ) (151) v(x,t) = ṽ(x,t) -ṽss (x; d 1 , d 2 , d 3 ) (152) η(t) = η(t) -ηss . (153)
Using ( 146) with Z = [ ũss , ṽss ] together with (142), (143) it is shown that the variables u and v satisfy

u t + ε 1 (x)u x = γ 1 (x)v (154) v t -ε 2 (x)v x = γ 2 (x)u . (155)
Setting x = 0 in (151), (152), and using (144), (147), and (148) we get that u(0,t) = qv(0,t) .

Setting x = 1 in (152) and using (145), (91), and (150) we get v(1,t) = -k P ṽ(0,t)k I η(t) + k I ηss . Using (152) for x = 0 together with (148) and ( 153) we arrive at

v(1,t) = -k P v(0,t) -k I η(t) . (157)
Using (153) and the fact that v(0,t) = ṽ(0,t) ,

(92) becomes η(t) = v(0,t) . (158) relation
Under the assumptions of Theorem 2 the zero equilibrium of (154)-( 157) and (159) is exponentially stable.

Simulations

Trajectory Generation for a Wave PDE

In this subsection, we illustrate our trajectory generation methodology with a wave PDE of the form (70)-(72). We choose the parameters of the system as

ε = 1 (160) h = 1 (161) b = -1 (162) g = 0 . ( 163
)
The reference for the output is chosen as ζ (t) = sin(3t). As in Example 1, this choice of parameters gives c 2 = c 4 = 0, and hence γ 1 = 0. Therefore, using relations (65)-( 68) we obtain L αα (x, ξ ) = 0 (164)

L αβ (x, ξ ) = 0 (165)

L β α (x, ξ ) = - 1 2 exp x + ξ 2 (166) L β β (x, ξ ) = - 1 2 exp x -ξ 2 . (167)
The reference trajectory z r for system (70)-( 72) is given by z r (x,t) = 1 37 (19 sin (3t + 3x) -3 exp(x) cos (3t -3x) + 3 cos (3t + 3x)

+18 exp(x) sin(3t -3x)) ,

which gives the following reference input W r (t) = 57 37 (cos (3t + 3)exp(1) cos (3t -3)) + 9 37 (exp(1) sin (3t -3)sin (3t + 3)) .

Figure 2 shows the evolution of the reference trajectory z r . Figure 3 shows the evolution of the spatial derivative of z r and, in particular, the control effort W r (t) = z r x (1,t) given by (169).

Trajectory Tracking

In this subsection, a simulation study for the system from Examples 1 and 2 is presented. The numerical approximation of the solution is computed with a two-step variant of the Lax-Friedrichs (LxF) method [START_REF] Shampine | Solving Hyperbolic PDEs in MATLAB[END_REF]. The reference for the output is chosen as y r (t) = cos(t). We add disturbances at the right-hand side of (61), (62) given by

Figure 1 :

 1 Figure 1: Evolution of the eigenvalues of (94) as a function of x (square and cross markers), and of the determinant of P(x) in (103) (star marker) for Example 2.

d 1

 1 (x) = 0.5 exp(x) (170) d 2 (x) = cos(2x) ,

Figure 2 :

 2 Figure 2: Solution to the trajectory generation problem for system (70)-(72) with parameters (160)-(163).

Figure 3 :

 3 Figure 3: The spatial derivative of the reference trajectory of Figure 2. Note in particular the reference input W r (t) = z rx (1,t) given by (169).

Figure 4 :

 4 Figure 4: The output v(0,t) of system (61)-(64) with parameters (134)-(137) under the control law (177) with gains (138), (139) (square marker) and with gains (138), k I = 0 (star marker) for the initial conditions (174)-(176). The single line is the reference output y r (t) = cos(t).

Figure 5 :

 5 Figure 5: Evolution of the state v of system (61)-(64) with parameters (134)-(137) under the control law (177) with gains (138), (139) for the initial conditions (174)-(176).

  ) together with the boundary conditions (147), (148) is a well-posed initial value problem for x. The equilibrium depends on d 1 , d 2 and d 3 . Let us denote this equilibrium by ũss (x; d 1 , d 2 , d 3 ), ṽss (x; d 1 , d 2 , d 3 ). From (145) it follows that the equilibrium value of Ũ, namely Ũss , satisfies Ũss = ṽss (1; d 1 , d 2 , d 3 )d 4 .

			(149)
	Using (91), and (148) with Z = [ ũss , ṽss ] , it follows from (149) that the equilibrium value of η,
	namely ηss , satisfies
	ηss = -	ṽss (1; d 1 , d 2 , d 3 ) -d 4 k I	.

Acknowledgments

The authors are deeply grateful to Antoine Girard and Christophe Prieur for many constructive suggestions.

$ This research was partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025)

together with constant additive disturbances on the boundary conditions (63), (64) given by d 3 = 0.5 (172) d 4 = 0.5 .

(173)

The initial conditions for u and v are chosen as the reference initial conditions given by ( 19), [START_REF] Gugat | Existence of classical solutions and feedback stabilization for the flow in gas networks[END_REF] for t = 0, perturbed by spatially-varying errors as

and the initial condition for η is chosen such that U(0

Figure 4 shows that the output of the system v(0,t) follows the desired trajectory under the PI controller given by

with gains (138), (139), and η(t) = v(0,t)cos(t). One can also observe that with only a P controller (i.e., when k I = 0 in (177)) there is a steady-state tracking error. Figure 5 shows the evolution of the state v.

Conclusions

We presented solutions to the trajectory generation and tracking problems for general 2 × 2 systems of first-order linear hyperbolic PDEs. We solved the motion planning problem with backstepping and the trajectory tracking problem with PI control. We proved exponential stability of the closed-loop system by constructing a Lyapunov functional.