Damien Raynaud

Benoit Hingray
email: benoit.hingray@univ-grenoble-alpes.fr

Baptiste François

Jean Dominique Creutin

Hingray Benoît

Energy droughts from variable renewable energy sources in European climates

Keywords: Energy droughts, Variable renewable energy, Energy demand, Hydro-climatic variability, Europe, Wind, Solar, Hydropower

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

Introduction

During the 21 st Conference of Parties (COP21), the 2015 United Nations Climate Change Conference in Paris, , 175 countries agreed on limiting the temperature increase due to global warming to 2°C above preindustrial levels. Such an ambitious goal necessitates a deep transformation of our societies and first and foremost a reduction of the anthropogenic greenhouse gas emissions. Many countries have already started their energy transition. In Europe, increasing the share of renewables in the electricity production mix is one of the main targets for the next decade. The European Climate Foundation even set a 100% renewable mix objective to be met by 2050 [START_REF]Roadmap 2050: a practical guide to a prosperous, low carbon Europe[END_REF]. In Europe, just like in most regions worldwide, this goal is at least physically realistic, since the resource in renewable energy balances the energy demand several times [START_REF] Hoogwijk | Global potential of renewable energy sources: a literature assessment[END_REF].. Variable renewable energies (VRE) and especially those driven by weather, namely wind, solar and hydro-power from river flow, are expected to play a key role in increasing the share of renewables. Indeed, the installed capacity of VRE is quickly growing worldwide [3]. Exploitable roughly everywhere, wind and solar resources are, for instance, already important contributors of the energy mix in Europe. For some countries, they largely contributed to the early achievement of the 27 % share of renewables targeted for 2030 by the European Council [START_REF] Šturc | Renewable Energy: Analysis of the Latest Data on Energy from Renewable[END_REF]. Hydropower from river flow, classically obtained from run-of-the river power plants (further referred to as RoR power plants), is often given less consideration, especially when compared to the key role of hydropower from large water reservoirs [START_REF] Piria | Norway's Key Role in the European Energy Transition[END_REF]. Even if individual RoR plants have fairly no balancing capacity and only small to very small power capacity (e.g. often smaller than a few MW), the overall RoR production is far from negligible in some regions [START_REF]Statistique suisse de l'Electricité 2013[END_REF] [START_REF] Gse | Statistical Report: Renewable Energy Power Plants in Italy, Technical report Gestore Servizi Energetici[END_REF]. It is also expected to increase significantly in the next years worldwide with the construction of new plants, the upgrade of old ones and the deployment of new technologies such as river hydrokinetic turbines [START_REF] Khan | Hydrokinetic energy conversion systems and assessment of horizontal and vertical axis turbines for river and tidal applications: A technology status review[END_REF] [START_REF] Sample | Scotland׳ s hydropower: Current capacity, future potential and the possible impacts of climate change[END_REF].

Wind power, solar power and hydro-power from river flow are by definition very sensitive to weather conditions, and thus present high space/time variability [START_REF] Engeland | Space-time variability of climate and hydro-meteorology and intermittent renewable energy production -a review[END_REF]. Moreover, wind and solar power are highly intermittent. Consequently, the integration of VRE in the power system is often difficult. They also make power systems rather vulnerable to the hydro-climatic variability and hydro-meteorological extreme events. These issues become even more critical for high shares of VRE [START_REF] Von Bremen | Large-scale variability of weather dependent renewable energy sources[END_REF] [START_REF] François | Increasing climate-related-energy penetration by integrating run-of-the river hydropower to wind/solar mix[END_REF].

In the last decade, many studies have focused on the intermittency of VRE production at small time scales [START_REF] Lave | High-frequency irradiance fluctuations and geographic smoothing[END_REF] [START_REF] Sørensen | Power fluctuations from large wind farms[END_REF] [START_REF] Sevlian | Detection and statistics of wind power ramps[END_REF]. These fluctuations result from the strong variability of weather variables and, to a lower extend, from the cut-in and cut-off thresholds of power generators that only function for a specific range of meteorological conditions (e.g. wind-speed based cut-in and cut-off thresholds for wind turbines). Covering timescales ranging from seconds to several hours, these high-frequency variations in power production lead to a number of critical operational issues (e.g. maintaining the system's stability in order to avoid the system collapse). This especially calls for the support of sufficient mechanical inertia in the system, fast responding back-up power and flexible storage facilities as well as demand-side management [START_REF]Renewable Energy Integration in Power Grids[END_REF].

Low VRE production conditions can be also problematic. They result from weather or hydrometeorological configurations with low VRE resources or from situations for which a production curtailment is necessary to avoid damage on power production units (e.g. when wind speed exceeds the cut-off threshold of wind turbines). They can induce critical situations in term of reliability of the power supply potentially and require the use of large back-up sources or energy storage facilities [START_REF] Vliet | Impacts of recent drought and warm years on water resources and electricity supply worldwide[END_REF] [START_REF] Gleick | Impacts of California's ongoing drought: hydroelectricity generation[END_REF]. Low VRE production conditions have been the centre of attention of numerous studies. They either aim to assess the risks related to low production sequences, in terms of probability of occurrence or duration for both current and future climate conditions, or to improve their predictability, with the objective of supporting the operational management of transmission system operators and guarantee system reliability.

Low river flows affect the hydropower production and a number of other water-related systems. Consequently, a large number of studies focused on the characterization of these conditions in terms of occurrence, durations and severity from local to continental scales (e.g. [START_REF] Hannaford | Examining the large-scale spatial coherence of European drought using regional indicators of precipitation and streamflow deficit[END_REF] for the European domain). A lot of effort is also made to improve their predictability from daily to seasonal scale [START_REF] Lima | Climate informed long term seasonal forecasts of hydroenergy inflow for the Brazilian hydropower system[END_REF] [START_REF] Haro-Monteagudo | Drought early warning based on optimal risk forecasts in regulated river systems: Application to the Jucar River Basin (Spain)[END_REF]. The persistence of wind calms and/or the characteristics of lowwind speed periods (for various thresholds) have been similarly analysed for a number of sites worldwide [START_REF] Lavagnini | Wind persistence in 135 Italian localities[END_REF][23] [START_REF] Deaves | The nature and frequency of low wind speed conditions[END_REF]. Recent studies also consider the occurrence of simultaneous low wind conditions across large areas in the context of modern transmission grids [START_REF] Leahy | Persistence of low wind speed conditions and implications for wind power variability[END_REF][26] [START_REF] Handschy | Is it always windy somewhere? Occurrence of low-wind-power events over large areas[END_REF]. Similarly, low solar power periods due to overcast conditions, persistent low level clouds or dust outbreaks have been recently analysed [START_REF] Pelland | Photovoltaic and solar forecasting: state of the art[END_REF][29] [START_REF] Rieger | Impact of the 4 April 2014 Saharan dust outbreak on the photovoltaic power generation in Germany[END_REF].

The impact of low VRE production periods is also likely to be increased if they co-occur with high energy demand [START_REF] Sinden | Characteristics of the UK wind resource: Long-term patterns and relationship to electricity demand[END_REF][32] [START_REF] Steiner | Critical weather situations for renewable energies-Part A: Cyclone detection for wind power[END_REF], a situation which can result from the fluctuation of some large scale climate phenomenon such as the North Atlantic Oscillation [START_REF] François | Influence of winter North-Atlantic Oscillation on Climate-Related-Energy penetration in Europe[END_REF].

In this study, we propose to complement and enlarge these analyses by characterizing and comparing the periods of low VRE production in Europe for the three main renewable power sources, namely wind, solar and hydro power from river flows. We consider the three sources separately and also combined in an energy mix. For the three sources, we use a same analysis framework and focus on what we call "energy droughts" whose concept rests on the analogy with the classical hydro-meteorological droughts, defined as long periods of very low river flows [START_REF] Lavagnini | Wind persistence in 135 Italian localities[END_REF] or periods with no or fairly no precipitation [START_REF] Mckee | The relationship of drought frequency and duration to time scales[END_REF] [START_REF] Hayes | The Lincoln declaration on drought indices: universal meteorological drought index recommended[END_REF]. We here propose two definitions of "energy droughts", considering either 1) Energy Production Droughts as uninterrupted sequences of days with low power production or 2) Energy Supply Droughts as uninterrupted sequences of days with a high production/demand mismatch. We here characterize Energy Supply Droughts within a hypothetical 100% VRE system where electricity generation comes from wind, solar and/or run-of-the-river hydro-power only. We do not base our analysis on site measurements but on high resolution gridded hydrometeorological datasets obtained for a 30-year historical period from satellite observations and from outputs of weather and hydrological models. This allows us to conduct our study on different European regions [START_REF] François | Increasing climate-related-energy penetration by integrating run-of-the river hydropower to wind/solar mix[END_REF] and thus explore how the characteristics of energy droughts also depend on the climatic spatial variability.

The meteorological datasets and the different weather-to-energy conversion models are presented in section 2 together with the two definitions of energy droughts. Section 3 gathers the results of this study and a comparison of the two droughts definitions. The results are discussed in Section 5 and section 6 concludes our study.

Datasets and models 2.1 Hydro-meteorological data

We consider twelve regions homogeneously spread over Europe and Maghreb and having a surface area of about 40000km² (Fig. 1). As discussed in François et al. [START_REF] François | Increasing climate-related-energy penetration by integrating run-of-the river hydropower to wind/solar mix[END_REF], they draw a picture of the large variety of climatic conditions existing in Europe with four main influences: a North-South gradient from Scandinavian to Mediterranean climate and a West-East contrast from Oceanic to continental influences. Several meteorological variables are required to simulate electricity production and demand time series. Temperature and precipitation data are based on on-site measurements and are extracted from the European Climate Assessment & Dataset (ECAD, [START_REF] Haylock | A European daily high-resolution gridded data set of surface temperature and precipitation for 1950-2006[END_REF]) at a daily time step and on 0.25 degrees grid spacing. Solar radiation series comes from the Surface Solar Radiation Dataset -Heliosat (SARAH) which relies on satellite measurement [START_REF] Müller | Surface Solar Radiation Data Set -Heliosat (SARAH) -Edition 1[END_REF]. It has a 0.05° grid and data are also extracted on a daily basis. Finally, we use pseudo-observations of wind speed derived from simulations of the Weather Research and Forecasting Model which has been forced by the ERA-Interim reanalysis [START_REF] Vautard | Regional climate model simulations indicate limited climatic impacts by operational and planned European wind farms[END_REF]. All these data are available on a 30-yr long research period, extending from 1983 to 2012. River discharges are highly regulated in most regions of Europe. In order to free ourselves from this anthropogenic influence and only depend on the hydro-climatic variability, we use unregulated discharge estimates. They are obtained with a distributed version of the GSM-Socont hydrological model [START_REF] Schaefli | A conceptual glaciohydrological model for high mountainous catchments[END_REF] on a 0.25° grid from each test region [START_REF] François | Increasing climate-related-energy penetration by integrating run-of-the river hydropower to wind/solar mix[END_REF]. Using the previous precipitation, temperature, and wind speed series, this model simulates the main hydrological processes which drive river discharge (e.g. snowpack dynamics, evapotranspiration, infiltration, runoff transfer).

Weather-to-energy conversion models

Daily time series of solar, wind and RoR power are simulated from weather and river discharge thanks to "weather to energy" conversion models. These models are based on generic production systems obtained from generic wind turbine, photovoltaic panel or RoR plant. For each energy source, we assume that all grid cells of the considered region have the same level of equipment. The overall regional production is obtained by summing the simulated production values of all grid cells.

For each grid, solar power (PPV) time series is derived from global solar irradiance Ieff (Wm -2 ) and air temperature Ta (°C) following Perpiñan et al. [START_REF] Perpinan | On the calculation of energy produced by a PV grid-connected system[END_REF]. The conversion of wind speed into wind power relies on a two-step process. Firstly, the 10m wind speed series are extrapolated to 80m (height of the generic turbine hub), assuming a logarithmic increase in wind speed with increase in altitude. Daily wind speed is then converted into wind power with the idealized power curve described in Francois et al. [START_REF] François | Increasing climate-related-energy penetration by integrating run-of-the river hydropower to wind/solar mix[END_REF]. To estimate RoR hydro power from river discharge PRoR (kW), we use a classic conversion function (Eq.1) using the water head h (m) defined as the altitude difference between the current grid cell and the lowest point of the region.

𝑃 (𝑡) = 𝜂 𝑔ℎ𝜌 𝑞(𝑡) (1)
where 𝜂 the efficiency of the turbine (%), 𝑞 the flow passing through the turbine (𝑚 . 𝑠 ), g the acceleration of gravity (=9.81 𝑚. 𝑠 ), 𝜌 the water density (=1000 𝑘𝑔. 𝑚 ). The water flowing through the turbine is modelled as a non-linear function of the river discharge as described in Francois et al. [START_REF] François | Increasing climate-related-energy penetration by integrating run-of-the river hydropower to wind/solar mix[END_REF]. The local runoff produced in any given grid is assumed to be ideally harnessed from its location of origin to the catchment's outlet.

The daily energy load required here to analyse the Energy Supply Droughts is finally estimated for each grid cell from daily temperature (Eq.2). Thus, we put aside all non-weather related influences on the demand, and only focus on the effects of temperature fluctuations. In this way, results from different regions are easily comparable whatever the population or their economic development. The daily energy load L (Wh) is estimated as:

𝐿(𝑡) = 𝑎 × [𝑇 -𝑇(𝑡)] + 𝑏 𝑖𝑓 𝑇(𝑡) < 𝑇 𝐿(𝑡) = 1 𝑖𝑓 𝑇 < 𝑇(𝑡) < 𝑇 𝐿(𝑡) = 𝑎 × [𝑇(𝑡) -𝑇 ] + 𝑏 𝑖𝑓 𝑇(𝑡) > 𝑇 , (2)
where THeat and TCool are the heating and cooling thresholds (15°C and 20°C respectively), b (Wh) is the base load for the temperature-independent range and aTheat (Wh °C-1 ), aTCool (Wh °C-1 ) are the weather sensitivity coefficients for heating and cooling respectively. The parameters were estimated using observed data of electricity demand from the European Network of Transmission Systems Operators of Electricity (ENTSOE, https://www.entsoe.eu/home/) (more details in François et al. [START_REF] François | Increasing climate-related-energy penetration by integrating run-of-the river hydropower to wind/solar mix[END_REF]).

Definition of energy droughts 3.1 Energy Production Droughts (EPD)

We define a low production period as a contiguous sequence of days during every day of which the production is below a given low-production threshold. This analysis requires the calculation of a daily Deficiency Index (DI) time series for each region and energy source. The DI time series simply reduces the power production to a binary time series equals to 0, when the daily production (P) is greater than the low-production threshold (P0), or to 1 when lower (Eq.3).

DI(j) = 1 if P(j)  P0 (3) DI(j) = 0 if P(j) > P0
In our case, the threshold is set to a given percentage of the mean daily production estimated over the whole period of analysis, i.e. 1983-2012. The duration of a low production period is defined as the number of consecutive days for which the DI index equals 1 (Fig. 2a).

Our definition is different from the one classically used for the analysis of hydrological droughts. In the latter case, the threshold changes periodically in time and usually corresponds to a given percentile of observed discharge data (e.g. the 90% percentile). It is usually estimated by using a moving seasonal window and all available observations falling within [42][43]. In our study, the threshold is time-invariant. Thus, droughts events are defined per-se and not as events that are unusual when compared to a calendar climatology. This definition allows the comparison of droughts regimes and severity obtained for different energy sources or different regions.

We will consider in turns two thresholds which correspond respectively to 0.2 and 0.5 times the mean daily production of the considered region/VRE source. We will refer to the corresponding low productions as moderate and severe "Energy Production Droughts" respectively (EPD in the following).

EPD are defined with respect to the mean daily production obtained for the considered region. They are thus defined per-se and do not depend on the electricity demand. Hence, they only represent a signature of weather and hydrometeorological variability.

Energy Supply Droughts (ESD)

As mentioned in Section 1, the severity of "energy droughts" can be reinforced (resp.

reduced) if they occur during higher (resp. lower) than usual demand periods.

An impact oriented assessment of energy droughts leads to consider Energy Supply Droughts, where the production/demand (im-)balance is accounted for. For the present analysis, we account for the production/demand (im-)balance that would be obtained for each region individually. We thus consider that each region is "autonomous", i.e. that there is no energy exchanges with the neighbouring regions. We further consider that all the regional production comes from solar, wind and/or RoR hydro power only and that the 30-yr mean regional VRE production equals the 30-yr mean regional electricity demand. This configuration corresponds to the so-called 100 % renewable scenario considered in a number of recent works (e.g. [11][12]).

We finally use the "copper plate" assumption within each region and thus consider that there is no energy transmission limitation within each region.

We here also define what we will refer to as low satisfaction periods, as sequences of days for which the non-satisfied demand is greater that a given fraction k of the mean electricity demand. This analysis also requires the calculation of a Supply Deficiency Index (SDI) time series for each energy source and region. This index is estimated on a daily basis as following (Eq.4):

SDI(j) = 1 if (D(j)-P(j)) > k.Dm SDI(j) = 0 if otherwise (4)
Where D(j) is the energy demand on day j, Dm is the mean daily energy demand over the whole 1983-2012 period. The main characteristics of low satisfaction periods are derived from the time series of the SDI index similarly to what has been described for Low Production Periods in Section3.1 (see Fig. 2b). We consider two thresholds values, k=0.5 and k=0.8, which correspond to what we will call moderate and severe Energy Supply Droughts respectively (ESD in the following).

The DI and SDI indices are very similar. Actually, the DI index corresponds to the SDI index for which the electricity demand would be constant in time and equals to the mean production Dm. In our case, it is easy to show that the thresholds considered for the definition of moderate EPD and moderate ESD refer to a same drought intensity level. The same applies to severe EPD and severe ESD. Thus, moderate (resp. the severe) droughts, can be compared for both definitions, making possible to estimate how the production/demand co-variability structure influences the characteristics of energy droughts.

A visual comparison between the two definitions of energy droughts is presented for one year-long time series of simulated hydro power production and energy demand on Fig. 2.

Fig. 2: "Energy Production Droughts" (a) and "Energy Supply Droughts" (b). Illustration for one year of hydropower power production (blue) and energy demand (black) in England. Production and demand time series have been normalized so that the annual average equals 1 for both series. Energy Droughts sequences correspond to the grey shaded bars. In (a) droughts correspond to a production smaller than 0.2 times the mean annual production. In (b) droughts correspond to a difference between demand and production greater than 0.8 (in this configuration the production falls below the dotted line which is defined from D(j)-0.8Dm). These two configurations correspond to the definition of what we call "severe" droughts.

Results

For the different energy sources, we now describe the number of days with low production and present the characteristics of EPD and ESD. Some figures only gather information for 3 regions (Norway (NO), Germany (GE), Andalucía (AN)) which draw a rather representative picture of energy droughts features for Europe.

Days in drought state

Fig. 3 presents the cumulative distribution functions (cdfs) of daily power production simulated for wind, solar and hydro power in NO, GE and AN. To make the comparison between energy sources easier, the times series have been scaled dividing them by the mean daily production over the whole 1983-2012 period. The two dashed lines represent the 0.5 and 0.2 thresholds used to identify the days in a moderate or severe drought state respectively.

Wind power classically exhibits the larger range of daily production values with heavytailed cdfs (e.g. in GE and AN). The proportion of days in a moderate drought state is rather similar between regions and of about 35%. Conversely, regional differences are much higher for days in a severe drought state with much fewer days in southern Europe (here in AN).

The variability of daily solar power production and the proportion of days in drought state are much higher in Northern Europe than elsewhere. This is obviously a direct consequence of the latitude influence on the space/time daylight variations in Europe. About 25% and 40% days are for instance below the severe and moderate thresholds respectively in NO while only few days are in drought conditions in AN (about 1% for the severe droughts threshold and 10% for the moderate drought one).

For hydropower, the cdf of daily production can have very different shapes from one region to the other. They have similar shapes in Alpine areas, Northern and southern European regions (here NO and AN) with a large proportion of days below the drought thresholds (15% below the 0.2 threshold and about 30% below the 0.5 one). In high elevation or high latitude regions, these days often occur in winter as a result of the typical snow-dominated hydrological regimes of river basins (not shown). The cdfs of hydro power production present less skewness in regions located in plains (here GE), leading to much smaller numbers of days in drought conditions.

Energy Production Droughts

The characteristics of EPD also depend on the energy source and on the region. Fig. 4 gathers, for the 12 European regions, the mean duration (x axis) and the mean annual frequency This representation gives a synthetic but partial picture of drought characteristics. The characteristics of individual drought sequences (e.g. duration and intensity) may obviously differ greatly between events. For illustration, Fig. 5 presents the cdfs of the droughts durations obtained for severe and moderate EPD in the NO, GE and AN regions. In the following analysis, this variability should be kept in mind. For both moderate and severe droughts, the homogeneity of wind power EPD characteristics is noteworthy: all individual droughts have very short durations (1 to 7days -Fig. 5). The main characteristics of moderate wind power EPD are also very similar from one region to the other (Fig. 4). Indeed, whatever the region, moderate droughts are numerous (from 50 to 60 events per year) and last from two to three days in average. The mean duration of severe wind droughts, smaller than 2days, is also very similar between regions. Conversely, their frequency is more variable with fewer droughts sequences in the Mediterranean regions. This result induces a much lower number of days in severe drought conditions for these regions.

Hydropower droughts are rather rare (less than 5 events per year). Individual hydropower droughts also present a large diversity of duration with relatively short episodes to extremely long ones which can last several months (Fig. 5). Conversely to wind EPD, short events (up to a few days) represent a small proportion of hydro-power droughts. Indeed, the hydrological processes and underground storage capacity of river basins work as powerful low-pass filters and ease most of the daily meteorological variability. The characteristics of hydropower droughts depend a lot on the region (Fig. 4). The longest droughts are found in Scandinavian and Mediterranean regions. In NO and AN, hydropower droughts can last up to 5 months. However, these EPD occur in different seasons: in summer in TU, AN and GR, in winter in NO, IT and FI (not shown). They respectively result from long periods with no or few precipitation and, for high elevation/latitude regions, from the storage in snowpack of most winter precipitation. Finally, the characteristics of hydropower droughts also significantly depend on the drought intensity level.

Severe hydropower EPD are roughly twice less numerous than moderate ones. They also last 30% to 50% time less. The regional specificities are roughly the same for both drought intensity levels.

The characteristics of solar droughts are in-between those of wind and hydropower droughts (Fig. 4). The number of solar drought events is medium in all regions. However, there are two different types of droughts with on one hand very short events, and on the other hand 2 to 3 months EPDs (Fig. 5). Short events correspond to cloudy weather situations, which mainly last a few days. Conversely, long lasting EPDs result from short daylight duration in winter. This effect is logically strengthened in Northern regions (here NO) where it is clearly visible for both moderate and severe EPD. In medium latitude countries, it only impacts moderate droughts whereas no effect can be found in Mediterranean regions.

Whatever the energy source, severe droughts are expected to be less numerous and to present lower mean durations than moderate ones. This is roughly the case for all regions.

However, solar EPD in Nordic regions do not respect this rule as a result of the bimodal distribution of drought duration mentioned previously. For these regions, when going from moderate to severe droughts, only short duration events are removed whereas long duration ones remain. It then results in an increase of the mean drought duration.

Energy Supply Droughts

As mentioned in Sec.1, low production periods can be more (resp. less) critical when cooccurring with high (resp. low) energy demand. Fig. 6 presents the differences in number of days in drought conditions, mean duration and mean annual frequency of droughts between EPD and ESD. For wind power droughts the mean duration of EPD and ESD are very similar whatever the region. The differences in droughts frequencies are also limited for most regions. However, the number of drought days and the annual frequency of drought events tend to be slightly smaller for ESD in Northern regions. It likely results from the partial synchronism between wind power and demand which both present a seasonal pattern with higher values in winter (higher wind speed and higher demand from heating systems). However, the production / demand covariability structure can sometimes lead to opposite changes for moderate and severe droughts (e.g. NO).

The characteristics of solar power droughts are strongly impacted by the energy demand.

The number of days in drought situation is higher everywhere, except in the southern regions. This result is consistent with the synchronism existing between low winter production due to short day length and high winter energy consumption in northern and central Europe. It is noteworthy that for most regions, higher numbers of drought days are associated with smaller frequencies of drought events (e.g. all regions for moderate droughts and Scandinavian regions for severe droughts). Naturally, this effect is balanced with higher droughts durations. Finally, for some remaining cases (e.g. RO for severe droughts), higher number of drought days result from more droughts sequences but with no changes in mean duration.

For hydro power, the differences between EPD and ESD are strongly region-dependent.

The number of droughts days slightly decreases in regions where hydrological regimes present higher river discharge in winter (e.g. FR, GE, GA, GR). ESD tend to be shorter than EPD. This results from the seasonal co-variability between precipitation and electricity demand with higher precipitation, discharge and hydropower production in winter, when heating systems are operating. The opposite co-variability structure happens for the other regions (continental, mountainous and Scandinavian regions) where a higher winter energy demand co-occur with low winter discharges. This latter configuration increases the number of days in drought conditions with shorter and more frequent ESD sequences in NO or longer sequences with similar annual frequency in RO for instance.

Whatever the energy source, small differences exist between EPD and ESD in southern regions. It is probably a consequence of a lower sensitivity of the demand to regional temperature changes. Indeed, the weather is relatively mild from autumn to spring, making the demand roughly weather-independent. In other countries, heating systems tend to boost the winter demand. On the other hand, the summer increase in demand from air conditioning is of relatively small magnitude in our temperature-demand model (Fig. 7).

Discussion

In Sec.4, we showed that the characteristics of Energy Droughts depend on both energy source and regional climate. They are also expected to depend on non-climatic factors such as technical components and structural choices for power supply systems. In this section we illustrate how Energy Drought's characteristics would be changed for power system using a mix of energy sources or including some storage facilities to balance the production/demand temporal mismatch.

Energy droughts for an energy mix

Mixing energy sources is known to reduce the variability of the production and classically improves the production/demand match [START_REF] Von Bremen | Large-scale variability of weather dependent renewable energy sources[END_REF]. For the 12 European regions considered here, François et al. [START_REF] François | Increasing climate-related-energy penetration by integrating run-of-the river hydropower to wind/solar mix[END_REF] estimated energy mixes based on hydro, solar and wind power, which maximize the mean daily production/demand match (c.f. share coefficients in Tab.1). We here illustrate how this mix may also influence the statistical characteristics of Energy Droughts. Fig. 8 presents the characteristics of EPD obtained with these regional mixes for both moderate and severe droughts. It should be compared to Fig. 4 which presents the results of single VRE sources. In all cases, a large decrease in drought severity is obtained with the mix. The number of days under drought conditions is drastically decreased, by a factor 2 to 10 depending on the region and on the single VRE source considered initially. This decrease results in a large reduction of both duration and frequency of energy droughts. For some regions, the combination of VRE sources even leads to the complete removal of severe droughts (e.g. AN). Thus, the complementarity that exists between VRE sources can make power systems much less prone to EPD. Fig. 7: Effect of a VRE mix (see Tab.1) on moderate and severe Energy Production Droughts. Figure to be compared with Fig. 4 (caption details in Fig. 4).

Tab.1: Regional energy mix proposed by François et al. [START_REF] François | Increasing climate-related-energy penetration by integrating run-of-the river hydropower to wind/solar mix[END_REF]. Share coefficients are given for solar power (SPV), wind power (SW) and hydropower (SRoR) for the 12 regions.

Energy droughts and storage

Storage systems, such as water reservoirs and batteries, can also ease the integration of VRE; balancing part of the missing production with the energy surplus stored on previous days [START_REF] Rasmussen | Storage and balancing synergies in a fully or highly renewable pan-European power system[END_REF][45] [START_REF] François | Effects of Increased Wind Power Generation on Mid-Norway's Energy Balance under Climate Change: A Market Based Approach[END_REF]. Fig. 9 presents the effect of storage on ESD characteristics in NO, GE and AN. Similar results are obtained for the other regions. We consider 4 simplified and idealized storage systems with different storage capacities corresponding respectively to 1, 7, 30 and 90 days of mean energy load. The energy transfer from and toward the storage system is simulated on a daily basis. Depending on the daily production/demand balance, the storage system either fills up, storing the energy surplus, or runs out, to satisfy more demand (see [START_REF] Francois | Complementarity between solar and hydro power: sensitivity study to climate characteristics in Northern-Italy[END_REF] for details). For simplification, we assume a perfect storage (no power loss).

For wind power, the frequency of drought episodes is efficiently lowered even with small storage capacity. When a large storage is involved, only few drought events remain. This result directly comes from the high frequency variations of wind power and the short durations of wind power droughts.

As mentioned previously, hydro power droughts have rather long duration in most regions. Therefore, small size storage systems do not lead to a significant reduction of droughts:

the number of drought days does not change much. However, small size storage leads to a reduction of droughts frequency, which is counterbalanced with an increase in mean drought duration. A significant reduction of drought severity is obtained for storage capacities larger than 30 days of mean demand. The seasonal based size (S = 90) is even more efficient but still does not lead to similar improvement as for wind power.

Similar comments apply for solar power. Small storage capacity does not change the number of drought days but often leads to longer and less numerous droughts events. Different modifications occur depending on the regional location in latitude. Even with large storage systems, Northern regions still suffer from long and numerous solar power droughts. Conversely, these critical periods are almost all removed in Mediterranean regions, even with medium storage capacity (here AN).

Conclusions

The similarities in stochastic nature and characteristics (occurrence, duration, intensity) between low production periods from VRE sources and low flow periods in hydrology call for the introduction of the concept of Energy Droughts.

In this study, we showed the strong variations of Energy Droughts characteristics between sources and within Europe. Wind droughts are very numerous but present very short duration, whatever the location. Solar droughts are of two types: 1) those that directly derive from the seasonal cycle of daylight which are long lasting especially in northern regions 2) those induced by the weather variability (changes in cloud cover) which usually last a few days maximum.

Hydropower droughts are much more region and time dependent as a result of the combined effects of meteorological droughts (snow accumulation or lack of rainfall / snow melt) and of the low pass filtering behaviour of river basins. Whatever the region, the complementarity that exists between wind, solar and hydro would make VRE mixes much less prone to energy droughts than systems based on one single VRE source.

The characterization of Energy Droughts is expected to enlighten decision makers in the design of future renewable power systems, including the provision of backup and/or storage solutions. As illustrated in this study, the analysis of Energy Droughts should take into account the electricity demand sensitivity to weather conditions. Indeed, it can potentially increase / decrease the socio-economic impact of drought events, depending on the co-variability structure between production and demand.

Our work is based on a number of assumptions, data and modelling choices, which potentially lead to some uncertainty in our results. The hydro-meteorological data used here are not observations but pseudo-observations (e.g. radiation from satellite products, surface winds from WRF). Of course, such data present limitations. However, conversely to observed data, pseudo-observations are available at high space -time resolution and classically cover rather long periods of time, even in ungauged regions. They allow having a more relevant view of unusual hydro-meteorological events. The temporal and spatial organisation (and next co-variations) of these pseudo-observations is also a priori consistent with the weather conditions of the 1983-2012 period. Thus, our analysis is expected to give a reasonable picture of the main features of energy droughts in different European climates for the last decades and of their specificities for the different types of VRE sources. However, further studies should consider the possibility to use the higher quality datasets (e.g. high resolution reanalyses of surface weather) that should be available in the coming years.

The choice of the weather-to-energy model is also another source of concern. Generic production systems have been considered for all regions. However, a variety of technological solutions which have their own weather-to-energy conversion relationships exist to harvest VREs. These relationships are classically highly non-linear and vary from one plant to the other.

For instance, this is the case for wind turbines (cut-in / cut-off thresholds in the wind speedwind power relationship) or for hydropower-plants which present an even greater variety with site-dependent specificities. Consequently our results are expected to significantly depend on these technical features. Similarly, the weather driven part of the energy demand was estimated with a unique model for all regions. However, the sensitivity of the demand to temperature can strongly vary from one region to the other. For instance, the winter sensitivity in FR and NO is actually twice higher than the average relationship estimated for our 12 European regions; most heating systems being based on electricity in these regions. Similarly, the summer sensitivity is roughly 5 times greater in IT and GR than elsewhere as a result of a much higher proportion of households equipped with air conditioning. The weather sensitivity of the demand obviously determines the differences between Energy Production Droughts and Energy Supply Droughts.

Further studies should explore how our results depend on the technical and socio-economic configuration. However, the main conclusions of our analysis are not expected to drastically change.

As mentioned previously, Energy Production Droughts are a signature of weather and hydrometeorological variability. Thus, they are defined per-se. Conversely, Energy Supply Droughts depend on a number of socio-economic and technical factors, including those pertaining to the energy mix and the power system considered. To provide some insight on Energy Supply Droughts, part of our analysis was carried out with the 100% renewable energy scenario. In this context, we assume that the resource is sufficient to provide the whole energy demand. This is likely not the case for all the configurations we considered. This would depend on the available mean resource itself which can be rather low in some regions (e.g. solar in the north) and on the mean demand which can be rather large in highly populated areas. The operational feasibility of such a 100% scenario is also obviously questionable, as discussed in a number of recent publications (see for instance [START_REF] Clack | Evaluation of a proposal for reliable low-cost grid power with 100% wind, water, and solar[END_REF] and [START_REF] Jacobson | The United States can keep the grid stable at low cost with 100% clean, renewable energy in all sectors despite inaccurate claims[END_REF]). Whatever those possible limitations, the 100% VRE scenario offers a convenient and simple analysis framework able to give informative perspective of a configuration to which a number of countries / regions are prone to head toward [START_REF]Roadmap 2050: a practical guide to a prosperous, low carbon Europe[END_REF]. In our analysis, it allows highlighting the interest of mixing of solar/hydro/wind sources, using the space/time complementarity between single sources. It also allows analysing how low resources sequences could interplay with high demand sequences as a result of the co-variability between weather variables. The analysis proposed here thus allows drawing a first but already informative picture of other energy droughts issues. Further works should consider a more realistic representation of the energy mix and power system. We here for instance focused our analysis on specific European regions assumed to be independent from each other. This is obviously not the case owing to the wide European electricity grid which connects regions and makes them less vulnerable to local weather variations. Further works should then consider Energy Droughts in a broader space-time framework to determine the characteristics of Energy Droughts over a continuum of spatial scales from local to European ones. The potential severity of European-wide events is of major interest for building a reliable pan-European power system. On the other hand, Energy droughts should also be studied at a very local scale, for a reliability analysis of minigrids and stand-alone power systems which are appropriate for remote and isolated communities and/or industries.

Energy Droughts are a key security issue of the energy supply challenge in a context where the share of VRE is continuously increasing worldwide. This call for a detailed characterization of their main statistical properties and extreme return levels similarly to what is classically done for meteorological and/or hydrological extremes. This characterization should account for the likely intensification of the hydrological cycle and of hydro-meteorological extremes expected from the ongoing climate change. Energy droughts will likely also challenge the operational management of the energy systems. To support their real-time management for instance, further works should also explore / understand the large-scale meteorological and/or atmospheric configurations prone to Energy Droughts and assess the possibility to forecast such events. We can likely expect that these different issues will gain increasing attention in the coming years.

Acknowledgement

Fig. 1 :

 1 Fig.1: Locations and boundaries of the 12 test regions. The two letters are further used as regional IDs.

Fig. 3 :

 3 Fig.3: Cumulative distribution functions (cdfs) of normalized power production for each VRE source in 3 European regions. Hydro -blue, Wind -orange, Solar -red. See Fig.S1 for results on all 12 regions

(

  y axis) of energy droughts. Additionally, the background colour scale gives the proportion of days under drought conditions discussed in Sec.4.1. With this representation, all points located on a descending diagonal and for a fixed background colour present the same number of drought days. The representation also makes possible to differentiate configurations with numerous but short droughts (points located in the left/upper side of the diagonal) from situations with rare but long lasting ones (points located in the right/bottom side of a given diagonal).

Fig. 4 :

 4 Fig.4: Annual frequency (y-axis) versus mean drought duration (x-axis) of Energy Production Droughts (EPD) for (top) moderate and (bottom) severe droughts. The colour shade gives the proportion of days experiencing energy drought conditions within the 30-year period (wind: orange; solar: red; hydro: blue). All points located on a descending diagonal at a chosen colour present the same number of days in drought conditions. The ID of each region is indicated within each dot.

Fig. 5 :

 5 Fig.5: Reverse cumulative distribution of the duration of Energy Production Droughts for both moderate (top: 50% level) and severe (bottom: 20% level) droughts in NO, GE and AN. Hydro power droughts (blue); wind power droughts (orange), solar power droughts (red). All panels show the percentage of drought events of the 30-yr period that exceed a given duration in days. See Fig.S2.a and FigS2.b for results on all 12 regions.

Fig. 6 :

 6 Fig.6: Difference in statistical characteristics from Energy Supply Droughts to Energy production Droughts (ESD-EPD). x-axis: difference in mean duration of drought sequences; y-axis: difference in annual frequency. Colour scale: difference in the number of days in drought state over the simulation period.

Fig. 6 :

 6 Fig.6: Temperature to energy demand function (dashed broken linear line) with seasonal and annual distributions of daily temperature (DJF-blue, MAM-green, JJA-ref, SON-orange, annual-black). See Fig.S3 for results on all 12 regions.

Fig. 8 :

 8 Fig.8: Effect of storage on Energy Supply Droughts characteristics. Mean number of drought episodes versus mean drought duration. Results are displayed for hydro-power (blue), solar-power (red), windpower (orange). The points associated to different storage capacities but from the same energy source have been connected in a logical order. They correspond to S = 1, 7, 30 and 90 days of mean energy load (S=1 in bold). The colour scale gives the mean annual number of days under drought conditions. See Fig.S4.a and Fig.S4.b for results on all 12 regions.

20

  This work is part of the FP7 project COMPLEX (Knowledge based climate mitigation systems

	for	a	low	carbon	economy;	Project	FP7-	ENV-2012	number:	308601;
	http://www.complex.ac.uk/).


Yohei Hosoe

Dimitri Peaucelle

Tomomichi Hagiwara

Linearization of Expectation-Based Inequality Conditions in Control for Discrete-Time Linear Systems Represented with Random Polytopes

Keywords: Stochastic dynamics, discrete-time linear systems, LMI, random polytopes, robust control

come

Introduction

In this paper, we deal with a class of discrete-time linear systems with stochastic dynamics described by random coefficient matrices. Such systems can be seen as a discrete-time linear case of random dynamical systems in the field of analytical dynamics [START_REF] Arnold | Random Dynamical Systems[END_REF]. The system class is useful for representing the influence of various stochastic phenomena such as packet interarrival times [START_REF] Paxson | Wide area traffic: The failure of Poisson modeling[END_REF] and failure occurrences [START_REF] Finkelstein | Failure Rate Modelling for Reliability and Risk[END_REF] in system dynamics. In addition, it is also considered to be compatible with the statistical approaches for estimation and modeling such as the ensemble Kalman filter [START_REF] Evensen | The ensemble Kalman filter: Theoretical formulation and practical implementation[END_REF], whose direction of research would have a potential of being a trend toward future control technologies.

As a traditional class of systems with stochastic dynamics that have been dealt with in the field of control, This paper was not presented at any IFAC meeting. Corresponding author Y. Hosoe. Tel. +81-75-383-2253. Email addresses: hosoe@kuee.kyoto-u.ac.jp (Yohei Hosoe), peaucelle@laas.fr (Dimitri Peaucelle), hagiwara@kuee.kyoto-u.ac.jp (Tomomichi Hagiwara).

Markov jump systems [START_REF] Costa | Discrete-Time Markov Jump Linear Systems[END_REF] are well known, whose dynamics is determined by a Markov chain. Although most of earlier studies on the systems were assuming that the underlying Markov chain only has finite modes (i.e., given with a countable state space), [START_REF] Costa | Stochastic stability of jump discrete-time linear systems with Markov chain in a general Borel space[END_REF] weakened this part of assumption so that a more general stationary Markov process (i.e., that with an uncountable state space) can be dealt with, through alternatively assuming the boundedness of the supports of random coefficient matrices (depending on the process). The research team involved in the present paper, on the other hand, has studied control of systems with dynamics determined by an i.i.d. stochastic process [START_REF] Hosoe | Robust stability analysis and state feedback synthesis for discrete-time systems characterized by random polytopes[END_REF][START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF], which is the time sequence of random vectors that are independent and identically distributed (i.i.d.). One might consider that our system class is included in that for Markov jump systems, since i.i.d. processes are a special case of stationary Markov processes. However, it is not true because the random coefficient matrices in our studies are allowed to have unbounded supports (e.g., normally distributed coefficients can be dealt with). Hence, the results for Markov jump systems and those in this paper do not cover each other, in general.

In this paper, we propose a technique useful for deriving numerically-tractable inequality conditions for control of systems with the dynamics determined by an i.i.d. process. In particular, to take account of possibilities of the presence of modeling errors even in such stochastic systems, we assume uncertainties in the random coefficient matrices that are represented with random polytopes [START_REF] Hug | Random polytopes[END_REF], as in [START_REF] Hosoe | Robust stability analysis and state feedback synthesis for discrete-time systems characterized by random polytopes[END_REF]. In [START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF], a Lyapunov inequality condition characterizing nominal stability was proposed and proved to be necessary and sufficient. However, the inequality involves decision variables contained in the expectation operation, which we call an expectation-based inequality. Such inequalities are difficult to transform toward advanced control compared with the conventional linear matrix inequality (LMI) conditions for deterministic systems [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF][START_REF] Ebihara | S-Variable Approach to LMI-Based Robust Control[END_REF]. Hence, [START_REF] Hosoe | Robust stability analysis and state feedback synthesis for discrete-time systems characterized by random polytopes[END_REF] discussed robust stability analysis and synthesis for the systems represented with random polytopes only under the approximation treatment of the distribution of the underlying i.i.d. process with its samples. The technique proposed in the present paper resolves this inconvenient situation in a direct fashion.

It is also remarked that the class of the systems with dynamics determined by an i.i.d. process is essentially the same as that of the systems with white parameters studied in De Koning (1992); Van Willigenburg & De Koning (1999) (which is a generalization of the systems with multiplicative noise [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF]). The contributions of our study compared to those earlier results mainly lie in our (re-)developing control theory for the systems based on LMIs through some new treatment of expectations. Our direction of the study naturally enables us, e.g., to achieve multi-objective control and to deal with random polytopes for the stochastic systems, as is the case with deterministic systems. The latter advantage about ease of handling of random polytopes will be clearer through the arguments in this paper.

The contents of this paper are as follows. Section 2 briefly describes the discrete-time linear systems represented with random polytopes. Section 3 shows two key lemmas, and extends the Lyapunov inequality in [START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF] toward robust stability by using the lemmas. In particular, the obtained robust stability condition is reformulated as a standard LMI. The same technique of using the two lemmas is also applied to robustly stabilizing state feedback synthesis in Section 4. Then, Section 5 compares this new approach for robust control with our earlier approach in [START_REF] Hosoe | Robust stability analysis and state feedback synthesis for discrete-time systems characterized by random polytopes[END_REF] through a numerical example.

We use the following notation in this paper. The set of real numbers, that of positive real numbers and that of non-negative integers are denoted by R, R + and N 0 , respectively. The set of n-dimensional real column vectors and that of m × n real matrices are denoted by R n and R m×n , respectively. The set of n × n symmetric matrices and that of n × n positive definite matrices are denoted by S n×n and S n×n + , respectively. The identity matrix of size n is denoted by I n ; the subscript will be dropped when the size is obvious. The Euclidean norm of the vector (•) is denoted by ||(•)||. The vectorization of the matrix (•) in the row direction is denoted by row(•), i.e., row(•) := [row 1 (•), . . . , row m (•)], where m is the number of rows of the matrix and row i (•) denotes the ith row. The Kronecker product is denoted by ⊗. The (block) diagonal matrix is denoted by diag(•). For the real square matrix M , He(M ) := M + M T , where M T denotes the transpose of M . The expectation (i.e., the expected value) of the random variable (•) is denoted by E[(•)]; this notation is also used for the expectation of a random matrix. If s is a random variable obeying the distribution D, then we represent it as s ∼ D.

2 Discrete-Time Linear Systems Represented with Random Polytopes

Let us consider the Z-dimensional discrete-time stochastic process ξ = (ξ k ) k∈N0 satisfying the following assumption.

Assumption 1 ξ k is independent and identically distributed (i.i.d.) with respect to k ∈ N 0 .

This assumption naturally makes ξ stationary and ergodic [START_REF] Klenke | Probability Theory: A Comprehensive Course[END_REF]. The support of ξ k is denoted by Ξ . By definition, Ξ ⊂ R Z , and Ξ corresponds to the set of values that ξ k can take.

With such a process ξ, consider the discrete-time linear system

x k+1 = A(ξ k )x k (1)
where x k ∈ R n , A : Ξ → R n×n , and the initial state x 0 is assumed to be deterministic. This system class is general in the sense that any discrete-time linear finitedimensional system with stochastic dynamics can be described by representation (1) with ξ satisfying Assumption 1 if the random coefficient matrix (denoted by A k ) determining the stochastic dynamics is i.i.d. with respect to the discrete time; in such a case, we can always take A(•) and ξ such that A(ξ k ) = A k .

In this paper, we assume that the random coefficient matrix A(ξ k ) in (1) has an uncertainty and is represented by

A(ξ k ) = A(ξ k ; θ) := L l=1 θ (l) A (l) (ξ k ) (2)
with a deterministic uncertain parameter vector θ (= [θ (1) , . . . , θ (L) ] T ) satisfying

θ ∈ E L := θ ∈ R L : θ (l) ≥ 0 (l =1, . . . , L), L l=1 θ (l) =1
(3) and the given mappings A (l) : Ξ → R n×n (l = 1, . . . , L) satisfying the following assumption.

Assumption 2 For each l = 1, . . . , L, the squares of elements of A (l) (ξ k ) are all Lebesgue integrable, i.e.,

E[A (l) ij (ξ k ) 2 ] < ∞ (∀i, j = 1, . . . , n), (4)
where

A (l) ij (ξ k ) denotes the (i, j)-entry of A (l) (ξ k ).
This assumption naturally makes the squares of elements of A(ξ k ) = A(ξ k ; θ) also Lebesgue integrable for each θ ∈ E L . Although Assumption 1 is stronger than that used in the earlier study [START_REF] Costa | Stochastic stability of jump discrete-time linear systems with Markov chain in a general Borel space[END_REF], in which ξ is assumed to be stationary and Markovian, Assumption 2 is weaker than another corresponding assumption used in the same study about boundedness of (not Ξ but) the support of A(ξ k ). This implies that Assumptions 1 and 2 would lead to such results that cannot be covered with the studies on Markov jump systems, in general. The system (1)-( 3) with ξ satisfying Assumptions 1 and 2 has been already dealt with in our earlier study [START_REF] Hosoe | Robust stability analysis and state feedback synthesis for discrete-time systems characterized by random polytopes[END_REF], and we call it a discrete-time linear system represented with (or characterized by) random polytopes [START_REF] Hug | Random polytopes[END_REF].

The random polytope {A(ξ k ; θ)} θ∈E L is a stochastic counterpart of the usual deterministic matrix polytope and is useful for representing uncertainties in random matrices. In particular, the introduced parameter θ can be reinterpreted as uncertainties in the distributions of random matrices in some cases. For example, it is stated in [START_REF] Hosoe | Robust stability analysis and state feedback synthesis for discrete-time systems characterized by random polytopes[END_REF] that the random matrix

A k =     0.6 0 a 2k -0.4 a 1k 0 0 -0.4 -a 1k a 2k     (5)
involving the random variables a 1k and a 2k obeying respectively the uncertain continuous uniform distribution with minimum 1.1-0.3δ 1 and maximum 1.4-0.3δ

1 (δ 1 ∈ [0, 1]) (denoted by U(1.1-0.3δ 1 , 1.4-0.3δ 1 )
) and the uncertain normal distribution with mean 0.6 and standard deviation 0.3δ 2 (δ 2 ∈ [0, 1]) (denoted by N(0.6, (0.3δ 2 ) 2 )) can be represented in the random polytope form without loss of generality, by associating θ with the uncertainties δ 1 and δ 2 under some appropriate vertices.

As the stability notion for system (1)-( 3), this paper uses the following second-moment exponential stability 1 [START_REF] Kozin | A survey of stability of stochastic systems[END_REF].

Definition 1 For a fixed θ ∈ E L , the system (1)-(3) (satisfying Assumptions 1 and 2) is said to be exponentially stable in the second moment if there exist a = a(θ) ∈ R + and λ = λ(θ) ∈ (0, 1) such that

E[||x k || 2 ] ≤ a||x 0 ||λ k (∀k ∈ N 0 , ∀x 0 ∈ R n ). (6)
If the system is (exponentially) stable for all θ ∈ E L , we say that the system is (exponentially) stable robustly with respect to E L .

3 Expectation-Based Inequality Conditions, Their Linearization, and Application in Robust Stability Problem

This section first reviews nominal stability conditions that are represented with what we call expectationbased inequalities, and then shows and discusses two key lemmas and their roles through extending the conditions toward robust stability. The inequality condition finally obtained in this section can be solved as a standard linear matrix inequality (LMI), and hence, we call the technique for such an extension using the two key lemmas a linearization technique. Our linearization technique is expected to be useful not only for the robust stability problem but also for other robust control problems associated with random polytopes.

First Key Lemma and Its Implication

The following theorem giving nominal stability conditions is shown in [START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF].

Theorem 1 Suppose ξ satisfies Assumption 1 and A (l) (ξ k ) (l = 1, . . . , L) satisfy Assumption 2. For each θ ∈ E L , the following three conditions are equivalent.

1. The system (1)-( 3) is exponentially stable in the second moment. 2. There exists

P = P (θ) ∈ S n×n + such that E[P -A(ξ 0 ) T P A(ξ 0 )] > 0. ( 7
)
3. There exist

P = P (θ) ∈ S n×n + and λ = λ(θ) ∈ (0, 1) such that E[λ 2 P -A(ξ 0 ) T P A(ξ 0 )] ≥ 0. ( 8
)
The expectation-based inequality ( 7) is a stochastic counterpart of the Lyapunov inequality for discrete-time deterministic linear time-invariant systems. In addition, ( 8) is also useful when the convergence rate with respect to the sequence

E[||x k || 2 ] k∈N0
is required, since λ in (8) corresponds to that in (6). For this reason, we will derive a robust stability condition from (8) with a technique using two lemmas shown later. However, the same technique can be used also for ( 7), which leads to another slightly different robust stability condition without λ.

The following is the first key lemma used in our linearization technique.

Lemma 1 For a given random vector ξ 0 (with support Ξ ) and a given mapping R : Ξ → S n0×n0 , the following two conditions are equivalent.

1. The following inequality holds.

E[R(ξ 0 )] ≥ 0 (9)
2. There exists T : Ξ → S n0×n0 such that

E[T (ξ 0 )] ≤ 0, (10)
T (ξ ) + R(ξ ) ≥ 0 (∀ξ ∈ Ξ ). (11)
PROOF. 1⇒2: T (•) = -R(•) readily leads to ( 10) and ( 11).

2⇒1: It follows from ( 11) that for each sample of ξ 0 , we have

T (ξ 0 ) + R(ξ 0 ) ≥ 0. This implies E[R(ξ 0 )] ≥ -E[T (ξ 0 )
] for random ξ 0 , which together with (10) leads to (9).

The newly introduced mapping T in this lemma enables us to make the random matrix R(ξ 0 ) uncontained in the expectation operation as in ( 10) and ( 11); the converse assertion (i.e., the assertion from condition 2 to condition 1) is also important and will be used soon. Although this lemma itself might seem almost trivial, it is actually very powerful in deriving inequality conditions associated with the use of random polytopes in system representation. For example, we can easily obtain the following corollary, which gives a stochastic counterpart of the LMI condition with S-variables (i.e., auxiliary variables) [START_REF] Ebihara | S-Variable Approach to LMI-Based Robust Control[END_REF] developed for deterministic systems.

Corollary 1 For given ξ 0 (with Ξ ), P ∈ S n×n

+ , N 1 : Ξ → R n1×n0 (n 1 ≥ n 0 ), U : R × S n×n → S n1×n1 and λ ∈ (0, 1), E[N 1 (ξ 0 ) T U (λ, P )N 1 (ξ 0 )] ≥ 0 (12)
holds if there exist S ∈ R n1×n2 , T : Ξ → S n0×n0 , J : S n0×n0 → S n1×n1 and N 2 : Ξ → R n2×n1 satisfying [START_REF] Ebihara | S-Variable Approach to LMI-Based Robust Control[END_REF] and

J(T (ξ )) + U (λ, P ) + He(SN 2 (ξ )) ≥ 0 (∀ξ ∈ Ξ ), (13)
N 1 (ξ ) T J(Π )N 1 (ξ ) = Π (∀Π ∈ S n0×n0 , ∀ξ ∈ Ξ ), (14)
N 2 (ξ )N 1 (ξ ) = 0 (∀ξ ∈ Ξ ). ( 15
)
PROOF. By post-and pre-multiplying N 1 (ξ ) and its transpose on (13), we obtain 16) by ( 14) and ( 15). This, together with Lemma 1, completes the proof.

T (ξ ) + N 1 (ξ ) T U (λ, P )N 1 (ξ ) ≥ 0 (∀ξ ∈ Ξ ) (
The significance of this corollary lies in the linearity of the inequality condition ( 10) and ( 13) in the operator N 2 (•), despite the nonlinearity of ( 12) in N 1 (•) (which has a close connection with N 2 (•) through ( 15)). This significance can be demonstrated by showing the following theorem, which gives a robust stability condition for the systems represented with random polytopes.

Theorem 2 Suppose ξ satisfies Assumption 1 and A (l) (ξ k ) (l = 1, . . . , L) satisfy Assumption 2. The system (1)-( 3) is exponentially stable in the second moment robustly with respect to E L if there exist S ∈ R 2n×n , P (l) ∈ S n×n + , T (l) : Ξ → S n×n (l = 1, . . . , L) and λ ∈ (0, 1) satisfying

E[T (l) (ξ 0 )] ≤ 0 (l = 1, . . . , L), (17)
T (l) (ξ ) + λ 2 P (l) 0 0 -P (l) + He S A (l) (ξ ) I ≥ 0 (l = 1, . . . , L; ∀ξ ∈ Ξ ). ( 18
)
PROOF. Fix θ ∈ E L . Fix also l and take

U (λ, •) = λ 2 (•) 0 0 -(•) , N (l) 1 (•) = I -A (l) (•) (19)
and then

J(•) = (•) 0 0 0 , N (l) 2 (•) = A (l) (•) I (20)
so that ( 14) and ( 15) with

N i (•) replaced by N (l) i (•) (i = 1, 2) hold. Then, it follows from (18) that J(T (l) (ξ )) + U (λ, P (l) ) + He(SN (l) 2 (ξ )) ≥ 0 (∀ξ ∈ Ξ ). (21)
Since this is true for each l, multiplying θ (l) on this inequality and ( 17) and taking the sums leads to [START_REF] Ebihara | S-Variable Approach to LMI-Based Robust Control[END_REF] and ( 13) with 14) and ( 15). Hence, it follows from Corollary 1 that

T (•) = T (•; θ) := L i=1 θ (l) T (l) (•), P = P (θ) := L i=1 θ (l) P (l) and N 2 (•) = A(•) I . Further- more, taking N 1 (•) = I -A(•) T T leads to (
E   I -A(ξ 0 ) T λ 2 P 0 0 -P I -A(ξ 0 )   ≥ 0, (22)
which is nothing but (8). This together with Theorem 1 completes the proof.

Although the direct extension of the inequality condition (8) toward robust stability is difficult because of its nonlinearity in A(ξ 0 ), Corollary 1 about the introduction of S-variables enables us to circumvent this issue as in the above proof even for the present stochastic systems case.

Obviously, Lemma 1 has played a crucial role in showing Corollary 1, and hence, the above arguments indeed demonstrate an aspect of the powerfulness of Lemma 1; the lemma (and the associated corollary) is useful for deriving a sufficient inequality condition described with the vertices of the random polytope, as in the present robust stability problem.

If we restrict S in (18) to be [0, G T ] T with G ∈ R n×n , present Theorem 2 immediately reduces to Theorem 2 in our earlier study [START_REF] Hosoe | Robust stability analysis and state feedback synthesis for discrete-time systems characterized by random polytopes[END_REF], which was derived with a completely different approach (i.e., without using Lemma 1). Hence, present Theorem 2 corresponds to a generalization of the earlier result and would enable us to achieve less conservative robust stability analysis. The simultaneous inequality condition consisting of expectation-based and infinitely-many matrix inequalities such as ( 17) and ( 18), however, is generally difficult to directly use for numerical analysis (even when the mappings T (l) are confined to some tractable class). Indeed, this was an obstacle and an analysis method using samples of ξ 0 was additionally discussed in [START_REF] Hosoe | Robust stability analysis and state feedback synthesis for discrete-time systems characterized by random polytopes[END_REF], in which it was required to evaluate the confidence level of the analysis result since the result inevitably depends on the used samples. Hence, if we stop our arguments here, the analysis based on present Theorem 2 would also have the same issue. To directly resolve this inconvenient situation so that the theoretically equivalent level of analysis (i.e., that using Theorem 2) can be performed without such a sample-based method, we provide in the next subsection the second key lemma and discuss a way of further transforming expectationbased inequality conditions. Such arguments will simultaneously lead us to a numerically-tractable inequality condition for robust stability that can be solved without the sample-based method despite no increase of the conservativeness in the corresponding robust stability analysis (compared with that using present Theorem 2).

Second Key Lemma and Numerically-Tractable Robust Stability Condition

The following is the second key lemma in our linearization technique.

Lemma 2 For matrices

M 1 =      M 11 . . . M 1p      (M 1i ∈ R r×m1 ), (23)
M 2 =      M 21 . . . M 2q      (M 2i ∈ R r×m2 ) ( 24
)
and

H ∈ R p×q , M T 1 (H ⊗ I r )M 2 = p i=1 q j=1 h ij M T 1i M 2j (25)
holds, where h ij is the (i, j)-entry of H.

This lemma can be shown with the direct calculation of the matrices, and corresponds to a generalization of the idea implicitly used in [START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF] for solving the expectation-based inequality conditions ( 7) and (8) (involving the decision variable P contained in the expectation operation) as standard matrix inequalities.

Before proceeding to the description of our linearization technique for expectation-based inequalities, we demonstrate how to use Lemmas 1 and 2 through deriving a numerically-tractable LMI condition for robust stability with the lemmas. Consider the (additional

) constraint 2 S 2 + S T 2 -P (l) > 0 (l = 1, . . . , L) (26)
on P (l) and S =:

[S T 1 , S T 2 ] T (S 1 , S 2 ∈ R n×n ).
Then, under this constraint, the Schur complement technique [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF] ensures that the inequality in ( 18) is equivalent to

T (l) (ξ ) + λ 2 P (l) + He(S 1 A (l) (ξ )) -(S T 1 + S 2 A (l) (ξ )) T • (S 2 + S T 2 -P (l) ) -1 (S T 1 + S 2 A (l) (ξ )) ≥ 0. ( 27
)
Hence, Lemma 1 leads us to the following corollary.

Corollary 2 For given S = [S T 1 , S T 2 ] T (S 1 , S 2 ∈ R n×n ), P (l) ∈ S n×n + (l = 1, . . . , L) and λ ∈ (0, 1), the following two conditions are equivalent.

1. There exist T (l) : Ξ → S n×n (l = 1, . . . , L) such that ( 17), ( 18) and ( 26) hold. 2. The following inequality and (26) hold.

E[λ 2 P (l) + He(S 1 A (l) (ξ 0 )) -(S T 1 + S 2 A (l) (ξ 0 )) T • (S 2 + S T 2 -P (l) ) -1 (S T 1 + S 2 A (l) (ξ 0 ))] ≥ 0 (l = 1, . . . , L) (28)
This kind of transformation of inequalities using Lemma 1 has never been considered in our earlier studies and will be a key in using Lemma 2 next.

We see that each term, after appropriate expansion, in the expectation operation in ( 28) can be represented in the form

M 1 (ξ 0 ) T HM 2 (ξ 0 ) for M i (ξ 0 ) = A (l) (ξ 0 ) or I n . For example, S 1 (S 2 +S T 2 -P (l) ) -1 S 2 A (l) (ξ 0 ) can be rep- resented as in that form with M 1 (ξ 0 ) = I n , M 2 (ξ 0 ) = A (l) (ξ 0 ) and H = S 1 (S 2 + S T 2 -P (l) ) -1 S 2 .
Let us temporarily focus on this form of expectations. Then, by using Lemma 2 under r = 1, this expectation can be transformed as

E[M 1 (ξ 0 ) T HM 2 (ξ 0 )] = n i=1 n j=1 h ij E[row i (M 1 (ξ 0 )) T row j (M 2 (ξ 0 ))], (29)
where h ij is the (i, j)-entry of H (∈ R n×n ). The expectation operation in the right-hand side of this equation does not contain decision variables. Here, if we can take X(l)

Ai ∈ R nl ×n and X(l) Ii ∈ R nl ×n (for some nl ) satisfying

E[row i (A (l) (ξ 0 )) T row j (A (l) (ξ 0 ))] = X(l)T Ai X(l) Aj , (30) E[row i (I n ) T row j (A (l) (ξ 0 ))] = X(l)T Ii X(l) Aj , (31)
E[row i (I n ) T row j (I n )] = X(l)T Ii X(l) Ij ( 32
)
for each l = 1, . . . , L, then Lemma 2 applied to the righthand side of (29) under r = nl leads us to

E[A (l) (ξ 0 ) T HA (l) (ξ 0 )] = X (l)T A (H ⊗ I nl ) X (l) A , (33) E[HA (l) (ξ 0 )] = X (l)T I (H ⊗ I nl ) X (l) A , (34)
(H =) E[H] = X (l)T I (H ⊗ I nl ) X (l) I , (35)
where

X (l) A =      X(l) A1 . . . X(l) An      ∈ R nn l ×n , X (l) I =      X(l) I1 . . . X(l) In      ∈ R nn l ×n . ( 36
)
The right-hand sides of ( 33)-( 35) are tractable from the viewpoint of using LMI techniques. The matrices X(l)

Ai and X(l)

Ii making this transformation possible can be obtained as follows. Take full row rank matrices X(l)

1A ∈ R nl ×(n 2 +1) (l = 1, . . . , L) satisfying X(l)T 1A X(l) 1A =E[X 2(l) 1A (ξ 0 )], ( 37
) X 2(l) 1A (ξ 0 ) :=[1, row(A (l) (ξ 0 ))] T [1, row(A (l) (ξ 0 ))]. (38)
Then, X(l) Ai can be obtained as

X(l) 1A =: X(l) 1 , X(l) A1 , . . . , X(l) An ( X(l) 1 ∈ R nl ×1 , X(l) Ai ∈ R nl ×n (i = 1, . . . , n)), (39) while X(l)
Ii can be constructed as

X(l) Ii = X(l) 1 row i (I n ) ∈ R nl ×n (i = 1, . . . , n), (40)
which can be confirmed by a direct calculation to satisfy (30)-( 32).

With the transformation of expectations ( 33)-( 35), the inequality condition (28) can be equivalently rewritten as

λ 2 P (l) + He( X (l)T I (S 1 ⊗ I nl ) X (l) A ) -(S T 1 ⊗ I nl ) X (l) I + (S 2 ⊗ I nl ) X (l) A T • (S 2 + S T 2 -P (l) ) ⊗ I nl -1 • (S T 1 ⊗ I nl ) X (l) I + (S 2 ⊗ I nl ) X (l) A ≥ 0 (l = 1, . . . , L). ( 41
)
Since (S 2 + S T 2 -P (l) ) ⊗ I nl > 0, the Schur complement technique again ensures that (41) (i.e., ( 28)) is equivalent to

λ 2 P (l) 0 0 -P (l) ⊗ I nl + He S 1 ⊗ X(l)T 1 S 2 ⊗ I nl X (l) A I nn l ≥ 0 (l = 1, . . . , L), (42)
where we used X (l)

I = I n ⊗ X(l)
1 , which can be confirmed from ( 40) and ( 36). Hence, this, together with Theorem 2 and Corollary 2, leads us to the following theorem.

Theorem 3 Suppose ξ satisfies Assumption 1 and A (l) (ξ k ) (l = 1, . . . , L) satisfy Assumption 2. The system (1)-( 3) is exponentially stable in the second moment robustly with respect to E L if there exist S i ∈ R n×n (i = 1, 2), P (l) ∈ S n×n + (l = 1, . . . , L) and λ ∈ (0, 1) satisfying ( 26) and ( 42), where X(l) 1 and X (l) A are given as ( 39) and ( 36) with full row rank matrices X(l) 1A ∈ R nl ×(n 2 +1) satisfying ( 37) and ( 38).

Once X(l)

1 and X (l)

A are calculated with the vertices A (l) (ξ 0 ) (l = 1, . . . , L), the inequality condition ( 26) and ( 42) can be seen as a standard LMI (for each fixed λ) with the decision variables S i (i = 1, 2) and P (l) (l = 1, . . . , L). Hence, we can analyze robust stability by minimizing λ with respect to the inequality condition through the bisection with respect to λ. This inequality condition is numerically more tractable than that in Theorem 2 in the sense that the sample-based method and the associated evaluation of the confidence level are not required. That is, the obtained λ readily corresponds to an upper bound of the convergence rate with respect to

E[||x k || 2 ] k∈N0 for each θ ∈ E L .
Remark 1 When L = 1, S 1 = 0 and S 2 = P (= P (1) ), the inequality condition (42) essentially reduces to that for nominal stability derived in our earlier study [START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF]. A transformation similar to that from (28) to (41) was also discussed in the study. The main differences from the earlier result in terms of this part are that the core idea of the technique used in the transformation was summarized explicitly as Lemma 2 so that the applicability of the technique is improved, and that we newly provided a key idea for dealing with the left-hand sides of ( 34) and ( 35) as well as (33) through the matrices X (l)

A and X

I . This is because the present paper dealing with robust stability came to handling the expectations of the entries of the matrix A(ξ 0 ) as in (34) through the introduction of the S-variable, unlike the earlier study in which it sufficed to deal with the expectations of the products of these entries as in (33); for concisely dealing with these two types of expectations, introducing the matrices X (l)

A and X (l) I is convenient, and is actually quite essential in the robust stabilization synthesis studied in the following section. Furthermore, the advantage of taking the full rank decomposition in (37) in reducing the associated computation cost was not claimed in the earlier study.

Summary of Linearization Technique Using Two Key Lemmas

In the preceding subsections, a numerically-tractable robust stability condition was derived by using Lemmas 1 and 2, which condition itself is new and useful. However, not such a specific condition but the technique used in its derivation is more significant and indeed constitutes a main contribution in this paper. The roles of the two lemmas used in the derivation can be summarized as follows.

Role of Lemma 1: To obtain θ-independent (i.e., uncertainty-free) inequalities associated with the vertices of the random polytope such as (28) from θdependent inequalities such as (8).

Role of Lemma 2: To obtain inequalities with decision variables uncontained in the expectation operation such as (42) from those with decision variables contained in the expectation operation such as (28).

Hence, even when we have only a θ-dependent inequality condition with decision variables contained in the expectation operation (i.e., a θ-dependent expectationbased inequality), these lemmas could lead us to a θindependent inequality condition with decision variables uncontained in the expectation operation from that original inequality condition. This is helpful for further using the conventional LMI techniques, and indeed, the LMI condition ( 26) and ( 42) was finally derived in the present robust stability problem case through the combination of the lemmas with the S-variable approach and the Schur complement technique. The whole technique of deriving LMI conditions in such a way with the two key lemmas is nothing but what we propose in this paper. Note that this technique itself is not specific to the present robust stability problem, and many other robust control problems could also be tackled in a similar fashion.

Remark 2 To the best of the authors' knowledge, no similar arguments arise in the field of Markov jump linear systems. This is because the inequality conditions such as the Lyapunov inequalities for standard Markov jump linear systems can always be written directly as standard matrix inequalities since any underlying expectation can be described merely as a weighted summation of a finite number of matrix terms. Hence, we encounter no obstacles in applying the standard LMI techniques such as the Schur complement technique directly to the matrix inequalities. Indeed, robust stability problems have already been discussed in de Souza ( 2006) for these systems only through the standard treatment.

Robustly Stabilizing State Feedback Synthesis

The Linearization technique discussed in the preceding section can be easily applied to the robustly stabilizing state feedback synthesis problem for the systems represented with random polytopes. Let us consider the Zdimensional process ξ satisfying Assumption 1 and the associated system

x k+1 = A op (ξ k )x k + B op (ξ k )u k , (43)
A op (ξ k ) = A op (ξ k ; θ) := L l=1 θ (l) A (l) op (ξ k ), (44)
B op (ξ k ) = B op (ξ k ; θ) := L l=1 θ (l) B (l) op (ξ k ) (45)
and ( 3), where

x k ∈ R n , u k ∈ R m , A (l)
op : Ξ → R n×n and B (l) op : Ξ → R n×m (l = 1, . . . , L). On the vertices of the coefficient matrices of the above system, we make the following assumption similar to Assumption 2.

Assumption 3 For each l = 1, . . . , L, the squares of elements of A Let us also consider the state feedback

u k = F x k (46)
with the static time-invariant gain F ∈ R m×n . Then, the closed-loop system can be described by ( 1)-(3) with

A (l) (ξ k ) = A (l) op (ξ k ) + B (l) op (ξ k )F (l = 1, . . . , L). (47) Note that if A (l) op (ξ k ) and B (l)
op (ξ k ) satisfy Assumption 3 then the above A (l) (ξ k ) also satisfy Assumption 2 (for each fixed F ). This section shows an inequality condition for designing F such that the closed-loop system is exponentially stable in the second moment robustly with respect to θ ∈ E L .

By Theorem 2, we see that the inequality condition for robust closed-loop stability can be given by ( 17) and (18) with A (l) (ξ k ) in (47). Since this inequality condition is difficult to linearize for the synthesis of F under the general class of S even in the case with deterministic systems, we restrict the class of S-variables to the following form3 with a given matrix A s ∈ R n×n (A s can actually be confined to a Schur stable matrix without loss of generality, although the details are omitted).

S = diag(G, G) A T s I n (G ∈ R n×n ) (48)
Then, G is nonsingular from P (l) > 0 and the lower right block of ( 18), and applying the congruence transformations using G -1 and diag(G -1 , G -1 ) respectively to ( 17) and ( 18), together with the change of variables

K = F V , Q (l) = V T P (l) V and W (l) (•) = V T T (l) (•)V for V = G -1
, leads to the following synthesis-oriented condition: a robust stabilization gain exists if there exist l) : Ξ → S n×n (l = 1, . . . , L), and λ ∈ (0, 1) satisfying

V ∈ R n×n , K ∈ R m×n , Q (l) ∈ S n×n + , W (
E[W (l) (ξ 0 )] ≤ 0 (l = 1, . . . , L), (49)
W (l) (ξ ) + λ 2 Q (l) 0 0 -Q (l) + He A T s I n A (l) op (ξ )V + B (l) op (ξ )K V ≥ 0 (l = 1, . . . , L; ∀ξ ∈ Ξ ). (50)
A robust stabilization gain can be obtained from the solution of this inequality condition as F = KV -1 . If we take A s = 0, the above condition immediately reduces to that in Theorem 3 in [START_REF] Hosoe | Robust stability analysis and state feedback synthesis for discrete-time systems characterized by random polytopes[END_REF]. Hence, the above condition is a generalization of the earlier result.

Such a new condition, however, is also incompatible with numerical synthesis as in the case of robust stability analysis based on Theorem 2. Hence, we further derive an equivalent numerically-tractable inequality condition from the above condition by using Lemmas 1 and 2.

First, under the additional constraint

V + V T -Q (l) > 0 (l = 1, . . . , L), (51)
the inequality condition ( 49) and ( 50) can be equivalently transformed to

E[λ 2 Q (l) + He(A T s (A (l) op (ξ )V + B (l) op (ξ )K)) -(V T A s + A (l) op (ξ )V + B (l) op (ξ )K) T (V + V T -Q (l) ) -1 • (V T A s + A (l) op (ξ )V + B (l) op (ξ )K)] ≥ 0 (l = 1, . . . , L) (52)
by using Lemma 1 and the Schur complement technique. We see that this inequality involves decision variables contained in the expectation operation. Hence, we next consider applying Lemma 2 to take out the decision variables from the expectation operation (recall the arguments in Subsection 3.2). Then, such transformation of the inequality and the application of the Schur complement technique finally lead us to the following theorem.

the confidence level for λ min < 1 (i.e., robust stability) was calculated and confirmed to be sufficiently high.

For the same system, we next design a robustly stabilizing gain F with Theorem 4. Before solving the LMI, we need to calculate the LMI coefficients Ȳ

(l) 1 , Y (l)
A and Y (l) B (l = 1, . . . , L) given by ( 54)-( 57). This calculation can be numerically performed with MATLAB and Symbolic Math Toolbox, and we obtained all the coefficients in 38.38 s under the computation environment described above. The associated ranks ml became 2 for l = 1 and 2, and 4 for l = 3 and 4, which are less than the worst case value n 2 +nm+1 = 13 in the sense of the size of the expectation in (56) for the present situation. With such coefficients and A s = 0, we minimized λ under ( 51) and ( 53) through a bisection with respect to λ. Then, we obtained 0.6837 as the minimal λ, and the solution leading to the gain F = KV -1 = [-1.0678, 2.5312, 1.6770]. In the bisection for the synthesis, we solved the LMI seventeen times, and the average time of the single LMI computation was 0.3079 s. Since the above value of minimal λ is the result obtained under the fixed A s (i.e., under the constraint on the S-variable), the post-synthesis analysis using Theorem 3 without such a constraint might provide us with a less conservative result. Hence, we numerically computed the minimization of λ under ( 26) and ( 42) in a fashion similar to the above synthesis. Then, we obtained 0.6806 as the minimal λ, which is less than the value obtained at the synthesis stage.

As we can see, the new robust stability analysis and synthesis framework successfully provided a better result (i.e., smaller λ) than that in [START_REF] Hosoe | Robust stability analysis and state feedback synthesis for discrete-time systems characterized by random polytopes[END_REF] revisited above, despite no need to use the sample-based method; the value of minimal λ obtained in the new framework is not an estimate, and readily corresponds to an upper bound of minimal λ satisfying (6) for all θ ∈ E 4 . Possible reasons for this improvement could be attributed to a quality issue of the generated samples of ξ 0 in the earlier study and the use of less conservative inequality conditions in the new framework; in the present situation, the latter reason is related essentially only with the restriction on mappings needed to introduce in the earlier study, since we took A s = 0 above (recall the arguments below (50)). In summary, we can say that although the earlier framework has an advantage that the exact model of the distribution of ξ 0 is not necessarily required (since only samples of ξ 0 are used), the new framework provides us with more accurate and less conservative results when the LMI coefficients can be calculated.

Remark 3 In the synthesis using Theorem 4, we have a degree of freedom in the selection of A s . By exploiting this degree of freedom, we can actually further improve the results in the new framework. As we can see from (48), S 1 and S 2 in Theorem 3 correspond to GA T s and G, respectively. Hence, if we take A s = (S -1 2 S 1 ) T with S 1 and S 2 obtained in the post-synthesis analysis, and if we design the gain F with such A s again, we can guarantee that the obtained result at least does not become worse than that in the analysis (since a trivial solution always exists). By using this idea once in the above numerical example, we obtained F = [-1.0494, 2.5652, 1.8049], which led us to 0.6794 as the minimal λ under ( 26) and (42).

Conclusion

This paper proposed a technique of using two key lemmas for linearizing parameter-dependent inequality conditions with decision variables contained in the expectation operation for discrete-time linear systems represented with random polytopes. As a demonstration of usefulness of the technique, we derived with the lemmas new numerically-tractable LMI conditions for robust stability of the systems. The obtained LMI conditions correspond to a solid stochastic counterpart (in other words, a generalization in the stochastic sense) of one of the latest robust stability results for deterministic systems [START_REF] Ebihara | S-Variable Approach to LMI-Based Robust Control[END_REF]. Hence, they enable us to achieve the same level of robust stability analysis and synthesis as such latest deterministic systems results, despite our dealing with systems represented with random polytopes. The improvement of the inequality conditions from the viewpoint of conservativeness was numerically confirmed through comparison with our earlier result. The proposed linearization technique is expected to contribute to deriving LMI conditions not only for the present robust stability problem but also for other robust control problems, and would facilitate the future development of the LMI-based control framework for the systems represented with random polytopes.

  k ) are all Lebesgue integrable.

Note that second-moment exponential stability (also called exponential mean square stability) is shown in[START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF] to be equivalent to second-moment asymptotic stability (i.e., asymptotic mean square stability) and quadratic stability for the present class of systems.

This does not actually increase the conservativeness in the associated analysis because of the positive definiteness of P(l) although the details are omitted.

This form itself does not actually decrease generality of S if As is also regarded as a decision variable.

Acknowledgements

This work was partially supported by JSPS KAKENHI Grant Number 17K14700.

Theorem 4 Suppose ξ satisfies Assumption 1 and A (l) op (ξ k ) and B (l) op (ξ k ) (l = 1, . . . , L) satisfy Assumption 3. There exists a gain F such that the closed-loop system (1)-( 3) and ( 47) is exponentially stable in the second moment robustly with respect to E L if there exist V ∈ R n×n , K ∈ R m×n , Q (l) ∈ S n×n + (l = 1, . . . , L) and λ ∈ (0, 1) satisfying (51) and

for a given matrix A s ∈ R n×n , where

A and Y

are given as

A1

. . .

B1

. . .

with full row rank matrices Ȳ (l) 1AB ∈ R ml ×(n 2 +nm+1) satisfying

In particular, V becomes nonsingular, and F = KV -1 is one such state feedback gain.

As in the case of Theorem 3, the inequality condition in the above theorem can be solved as a standard LMI (for each fixed λ) once

A and Y (l)

B (l = 1, . . . , L) are calculated. Hence, we can design a robustly stabilizing state feedback by using this theorem, without using the samples of ξ 0 . The increase of conservativeness in this inequality condition, compared to that in Theorem 3, could be led only by the fixed A s in (48). In other words, if A s is also regarded as a decision variable, the synthesisoriented inequality condition in Theorem 4 is theoretically no more conservative than the analysis-oriented inequality condition in Theorem 3 (although appropriate A s is generally difficult to globally search for).

Numerical Example

This section demonstrates effectiveness of the robust stability analysis and synthesis framework developed with the linearization technique using two lemmas, through the comparison with the authors' earlier approach in [START_REF] Hosoe | Robust stability analysis and state feedback synthesis for discrete-time systems characterized by random polytopes[END_REF]. Let us consider the process ξ satisfying Assumption 1 given by

(i.e., Z = 2) and the vertices

(i.e., L = 4) for the open-loop system (3), ( 43)-( 45). This system was dealt with in Hosoe et al. ( 2018) and its A matrix A op (ξ k ; θ) can be identified with the uncertain random matrix (5) under

We first briefly revisit the numerical results obtained in the earlier study [START_REF] Hosoe | Robust stability analysis and state feedback synthesis for discrete-time systems characterized by random polytopes[END_REF] for the above system. In the earlier study, we carried out synthesis and post-synthesis analysis through approximating the distribution of ξ 0 with its samples (i.e., with the empirical distribution), although the details are omitted not to make the arguments verbose. Such a samplebased method then led us in the study to the gain F = [-1.4267, 2.9341, 0.0609] as the optimal solution (i.e., that minimizing λ) for the synthesis under 100 generated samples of ξ 0 (and some restriction on mappings corresponding to W (l) in ( 49) and ( 50), which was needed at the numerical stage in association with the calculation of the confidence level in the post-synthesis analysis). The computation time for this synthesis was 70.06 s under the environment with MATLAB, YALMIP [START_REF] Löfberg | YALMIP: A toolbox for modeling and optimization in MATLAB[END_REF] and SDPT3 [START_REF] Tütüncü | Solving semidefinite-quadratic-linear programs using SDPT3[END_REF] running on a laptop equipped with 8.00 GB RAM and Intel(R) Core(TM) i7-5600U CPU @ 2.60 GHz; we also used this environment for computations in this paper. Since the above synthesis result depended on the generated samples, postsynthesis analysis was additionally performed with 1000 different samples of ξ 0 , and λ est = 0.7492 was obtained as an estimate of the true minimal λ (denoted by λ min ) satisfying the earlier robust stability condition (with a restriction on mappings); this estimate corresponds to that of an upper bound of minimal λ satisfying (6) for all θ ∈ E 4 . With that estimate and some other information,