Emmanuel Vinot
email: emmanuel.vinot@ifsttar.fr

TIME REDUCTION OF THE DYNAMIC PROGRAMMING COMPUTATION IN THE CASE OF HYBRID VEHICLE

Keywords: Deterministic Dynamic Programming, Solving Algorithms, Hybrid Vehicle

Deterministic Dynamic Programming is frequently used to solve the management problem of hybrid vehicles (choice of mode and power sharing between thermal and electric sources). However, it is time consuming and thus difficult to use in global sizing optimization or in parametric studies. This paper presents a comparison between three methods to compute the DDP problems. These methods are applied on the well known case of the Toyota PRIUS. It proves that a dense matrix method can reduce the computation time by up to 10 compared to classical solving methods.

Introduction

Hybrid vehicles are an effective solution to reduce CO 2 and pollutant emissions [START_REF] Husain | Electric and hybrid vehicles: design fundamentals[END_REF][2] [START_REF] Miller | Propulsion systems for hybrid vehicles[END_REF]. In such architecture, the electrical driveline provides additional functionality to the drive train. Among these functionalities, regenerative braking and electrical propulsion mode are the most efficient to improve systems efficiency. Combined electrical and thermal propulsion and battery recharge by the engine also allow substantial gain. Finally, as in series parallel hybrid, engine and wheel speed decoupling is another possibility to reduce fuel consumption.

These different capabilities have to be used in the best manner to reduce fuel consumption as much as possible. The choice of the operating mode (i.e. the power energy management) is thus a key point of the efficiency of hybrid architecture [START_REF] Chau | Overview of power management in hybrid electric vehicles[END_REF][5] [START_REF] Koot | Energy management strategies for vehicular electric power systems[END_REF]. It has to fix the electric or hybrid mode, the battery current and the engine operating points.

Moreover, the fuel consumption of an hybrid vehicle can only be announced taking into account the battery charge or discharge over a driving cycle long enough to be representative of real uses conditions. Thus, the power management law cannot be simply developed to minimize fuel consumption on each step of time but has to respect a global amount of charge of the battery.

An optimal method is then very useful to determine the best energy management law over a known in advance representative driving cycle [START_REF] Chau | Overview of power management in hybrid electric vehicles[END_REF] [START_REF] Sciarretta | Control of hybrid electric vehicles[END_REF] [START_REF] Koot | Energy management strategies for vehicular electric power systems[END_REF]. These methods are obviously applied offline as the drive cycle needs to be known. Optimal off-line management remains nevertheless very useful to developed efficient online management laws [START_REF] Koot | Energy management strategies for vehicular electric power systems[END_REF], to perform objective comparisons between hybrid solutions [START_REF] Vinot | HEVs comparison and components sizing using dynamic programming[END_REF], to size components [START_REF] Akli | Integrated optimal design of a hybrid locomotive with multiobjective genetic algorithms[END_REF] [START_REF] Kim | Power management and design optimization of fuel cell/battery hybrid vehicles[END_REF], and in optimal global design process [START_REF] Romaus | Optimal energy management for a hybrid energy storage system for electric vehicles based on stochastic dynamic programming[END_REF] [START_REF] Reinbold | Global optimization of a parallel hybrid vehicle using optimal energy management[END_REF] to avoid the unknown and non-monotonic impact of the management on the fuel economy [START_REF] Kim | Power management and design optimization of fuel cell/battery hybrid vehicles[END_REF] [START_REF] Desai | Optimal design of a parallel Hybrid Electric Vehicle using multi-objective genetic algorithms[END_REF]. Two off-line methods are currently used in order to calculate the optimal management; the Pontryagin's minimum principle [START_REF] Delprat | Control of a parallel hybrid powertrain: optimal control[END_REF] and the Discrete Dynamic Programming (DDP) [START_REF] Bellman | Dynamic Programming[END_REF] [START_REF] Kirk | Optimal control theory[END_REF]. These two methods can take a long time to achieve the optimal management over a drive cycle of several hundreds of seconds. This can be a critical point especially in comparative or parametric studies and even more critical in a design process [START_REF] Romaus | Optimal energy management for a hybrid energy storage system for electric vehicles based on stochastic dynamic programming[END_REF] [START_REF] Reinbold | Global optimization of a parallel hybrid vehicle using optimal energy management[END_REF] when fuel consumption of thousands of solutions can be calculated. It is then a serious issue to reduce the computer effort as much as possible during such calculation.

In the scope of this paper, the authors focus on DDP [START_REF] Bellman | Dynamic Programming[END_REF] [START_REF] Kirk | Optimal control theory[END_REF] [START_REF] Scordia | Global optimisation of energy management laws in hybrid vehicles using dynamic programming[END_REF] [START_REF] Perez | Optimization of power management in an hybrid electric vehicle using dynamic programming[END_REF] [START_REF] Vinot | Improvement of an EVT-based HEV using dynamic programming[END_REF] and investigate the different possibilities to reduce the calculation time. The DDP principle (section 2) allows to find the optimal State of Charge (SOC) trajectory of the battery to minimize the fuel consumption on a drive cycle known in advance. This is performed in a restrictive SOC versus time meshed area limited by battery maximum charge and discharge possibilities. Another specification, when applying DDP to hybrid vehicle is the necessary consideration of the electric mode (discrete state). When sampling the SOC you have to be sure that one of the samples corresponds to the electric mode with no fuel consumption associate (section 3).

Several manners to improve time calculation can then be investigated:

-Reducing the size of the graph as much as possible (section 3) -Building the graph to ensure that the electrical modes correspond to existing edges (section 3). -Using vectorial cost calculation. So that all edge costs of the graph are calculated all at once before solving the graph (section 5). -Improving the manner to solve the DDP problem avoiding loops and big matrix (section 4),

In this paper, the authors present the method they used to reduce the DDP calculation time to a handful of seconds (3 to 30). The four precedent points are treated. Then an example on the well known Toyota Hybrid PRIUS [START_REF] Kamiya | Development of traction drive motors for the Toyota hybrid system[END_REF][20] is presented.

Dynamic Programming applied to SP_HEV

Regarding hybrid vehicle management, an optimization problem can be formulated to minimize the global fuel consumption over a global driving cycle respecting a certain amount of battery discharge [START_REF] Bellman | Dynamic Programming[END_REF][15] [START_REF] Scordia | Global optimisation of energy management laws in hybrid vehicles using dynamic programming[END_REF][17] [START_REF] Vinot | Improvement of an EVT-based HEV using dynamic programming[END_REF]. The state variable is commonly the state of charge (SOC) defined as:

∫ - = T t F batt o o dt t I C t SOC T SOC η ) ( 3600 100 ) ( ) ( ( 1)
Where t 0 and T are respectively the initial and final time of a drive cycle, C batt is the battery capacity in Ah.

The global fuel consumption is:

∫ = T t SOC f o dt t m J ) , ( δ  ( 2)
Where δ SOC is the SOC variation at time t and f m  the fuel rate.

The global constraint on SOC variation has to be respected:

K SOC SOC SOC final initial = - = ∆ ( 3)
Using discrete dynamic programming principle, the problem formulation is the following:

The time is sample in N steps. The SOC is sample in M step between SOC min and SOC max (Fig. 1).

Solving this problem using DDP means to find the best SOC trajectory in a SOC versus time area. The SOC constraint is guaranteed if the limits of the space are built to cross on a specified final state of charge (Fig. 1, section 3).

To solve this problem an iterative method using Bellman principle [START_REF] Bellman | Dynamic Programming[END_REF] is currently used. A cost is associated to each edge between two consecutive points at time k and k+1 (Fig. 10). The fuel consumption on a trajectory between two points is simply the sum of the edge cost (fuel rate) on this trajectory.

Bellman principle stipulates that every sub-trajectory of an optimal trajectory is optimal [START_REF] Bellman | Dynamic Programming[END_REF]. It is therefore easy to find the best trajectory over the cycle using iterative methods.

The edge cost from point (k,i 1 ) to point (k+1,i 2 ) is the instantaneous fuel rate

i2)) 1, + (k i1) ((k, f m 
which depends on the soc variation between the beginning and the end of the edge δsoc. For instance we assume that

i2)) 1, + (k i1) ((k, f m 
depends only on the soc variation and vehicle speed, an example is given in section 5.

Starting at point (1, i 0 ) with a fuel consumption 0 * ) , 1

( 0 = i J
, the general principle is the following: At time k for a state SOC(i 1 ), i.e in point (k,i 1 ), the minimum fuel consumption from the beginning to this point is noted * ) , ( 1 i k J , and is called "cost to go" for the following of this paper. It is known for each point in the column k (i.e for i1 between i lim_min (k) and i lim_max (k) (Fig. 2) Using the Bellman principle the cost to go from beginning to a point (k+1,i 2 ) of time k+1 may easily be determined. The principle is to add the costs of each edge J((k,i 1 ) →(k+1,i 2 )) reaching the point (k+1,i 2 ) to the cost to go of the point (k,i 1 ) and take the minimum. Repeating this operation for each point (k+1,i 2 ) of the column k+1 gives the best trajectories from beginning to time k+1.

The optimal cost to go at time k+1 is thus given by:

) ( ) ( ] [ min ] , [ 2 max 1 2 min 1 1 * ) , ( ) , 1 ( ) , ( 1
* ) , 1 ( 1 2 1 2 i i i i i for J m i J i k i k i k f i i k i + = + = + → +  ( 4)
At each step of time, for each point of the column, the index i 1 which minimizes the fuel consumption is known [START_REF] Chau | Overview of power management in hybrid electric vehicles[END_REF]. Thus it is easy, at the end of the process to determine the best trajectory.

Note that point (k+1,i 2 ) can only be reached by a limited number of edges due to system limitation.

The battery has a limited rate of charge or discharge depending on the battery itself and on the system limitation. Thus the spread of possible edges reaching point i 2 is limited between i 1min (i 2 ) and i 1max (i 2 ). Section 4 of this paper focuses on different manners to solve the DDP problem and on the most efficient one in the cases of hybrid vehicle problems.

Graph Boundaries

An important point to reduce the computation effort is obviously to restrain the area of study. In our case, as the state variable is the battery state of charge the area is limited by a maximum discharge limit and a maximum charge limit Fig. 3.

The maximum discharge limit at each time corresponds to the zero emission mode (electric mode). That means that the vehicle is propelled only by electricity using the energy stocked in the battery. In the case of regenerative braking, this limit corresponds to recharging the battery only by the wheel (no energy is added by the engine).

The maximum charge limit at each time is calculated considering that the engine provides all its possible power, a part is going to the wheel to propel the vehicle, the rest recharges the battery. Note that this remains true in the limits of systems and components capabilities.

The last limits to be applied are those of the battery SOC itself Fig. 3. On a hybrid vehicle the battery remains relatively small, thus it may be rapidly completely charged (SOC=100%) or discharged (k,i 1 )

(k+1,i 2 ) i 1max (i 2 ) i lim_min (k+1) i 1min (i 2 ) i lim_max (k+1) i lim_max (k) i lim_min (k) δsoc (1,i 0 )
Spread of possible edge (SOC=0%). Moreover, for battery ageing considerations, a battery management system usually limits the minimum and maximum value. Common values are (40% and 80 %, [START_REF] Vinot | Model simulation, validation and case study of the 2004 THS of Toyota Prius[END_REF]).

Fig. 3 Battery SOC limitation

Once the limits of the graph are defined, the graph has to be meshed in the two axes. The easiest way is to sample the time in N step of T e starting at zero, and the SOC in M step between SOC min and SOC max starting at initial SOC value SOC 0 , then you obtain a regular mesh (Fig. 5). The drawback of such a grid is that the electrical mode (corresponding to maximum discharge limits) has virtually no chance to correspond to an effective edge. In fact, as the soc variation δsoc elec of this mode is perfectly known at each time by the drive cycle, it has few chances to match the SOC variation of an existing edge. This is of great importance; not considering electrical mode in hybrid vehicles implies great error on fuel consumption and does not allow making the strategy choice. A solution retained in [START_REF] Scordia | Global optimisation of energy management laws in hybrid vehicles using dynamic programming[END_REF] is to then consider that the nearest edge is the one corresponding to the electrical mode and to affect a cost of zero to this edge. Doing this creates errors on the SOC variations (Fig. 4) which are then cumulated along the whole cycle. Thus the final SOC constraint is not respected.

A solution is to mesh the area in such manner that the electrical mode corresponds to the edges of the graph (Fig. 5). For that, in each column, the grid starts from the maximum discharge limit (which is the electrical mode) and goes up. The SOC sample remains the same for each column. Once the graph boundaries have been constructed to minimize the meshed area (taking into account systems limitation) and to include all electrical mode, the cost of each edge of the graph has to be calculated.

Solving algorithms

Once the edge cost calculation is performed, the problem is to solve the DDP for a graph with a weighted oriented edge. Different manners can be imagined to solve this problem, from intuitive manner to complex matrix manner. The fact is that avoiding loops in the algorithms and too big matrix may be a good way to improve the algorithm efficiency especially using Matlab software.

This section deals with three different manners to solve the DDP problem:

(i) A first intuitive method (section 4.A) which consists in applying directly the principle described in section 2, using a loop on each point of a column inside the loop at each time k. (ii)

A matrix method (section 4.B) which treated each point of a column in a vectorial manner but used great sparse cost matrix with a lot of non defined elements. (iii) A third method (section 4.C) based on dense matrix of index of cost and weight.

A. « Point to point » method

This method is one of the more intuitive. The idea, in the iterative principle, is to deal separately with each point of a column at time k+1. A loop is then used for each point of the column (algorithm Fig. 7). On a point i 2 of column k+1, the cost of each edge (fuel rate) which can reach this point is added to the cost to go of these points. The trajectory and the corresponding fuel consumption corresponding to the minimum fuel consumption J * (k+1,i2) is selected.

The advantage of this method is that only the values of "valid" edges are considered. The main drawback is the loop on i 2 included in the time main loop on k (algorithm Fig. 7), and the presences of a call to a minimum function inside the loop. The computer effort used by this algorithm is detailed Table 3 andTable 4. An idea to reduce the computer effort consists in avoiding calling minimum function in the loop on i 2 . For that, a matrix containing all the possible weights of each point of the column is built (section 4.B). 7 Algorithm use in a point to point method B. Matrix method This method Fig. 8 is based on an M*M matrix calculation [START_REF] Mensing | Vehicle trajectory optimization for application in ECO-driving[END_REF] with M the maximum point numbers in a column of the graph. The principle is at each time (k+1) to add the cost to go at time k to the cost of each edge from time k to time k+1. If the edge between two points is not possible (due to the systems operation) an infinite weight is applied. The advantages of this method are its simplicity and you obtain in one operation the best cost and the index of the best trajectory (Fig. 9). For that you only have to add the cost matrix C to the cost to go matrix J * k (composed of the best previous fuel consumption at time k, and take the minimum of each lines of this matrix in one call to a minimum function with a matrix on argument (Fig. 9).

J i i i i i i k i k i i k f i k i k k k i k ] min[ ] , [ index) ing correspond
+ = → = → = ∞ = + → + + + + +  Fig.
The drawbacks are the size of the matrix, typically 1000 2 to 5000 2 elements. This implies great computer efforts to create the matrix and performed operation on it (Fig. 9). Moreover you still need a loop on i 2 inside the main loop to build the cost matrix C which is extracted of the global cost matrix C m . Battery SOC maximum charge limit maximum discharge limit

k k+1 i 2 i 1_sup (i 2 ) i lim_min (k+1) i 1_inf (i 2 ) i lim_max (k+1) i 0 i 1_inf (i 2 ) i 1_sup (i 2 ) i 2
Using such a method leads to a highly sparse matrix (20 defined values on the diagonal on a 1000*1000 matrix for example). Thus another method using as much as possible dense smaller matrix is presented section 4.C.

Fig. 8 Matrix method end C J index J end k i i C i i i C i i i for k k for k M k M k i e i

C. Index matrix method

The idea is instead of using sparse M*M matrix (section 4.B ), using a smaller dense cost matrix. This matrix contains quite only defined value (Fig. 10) and has a size max(S max (k)* M. In fact, the lines are only composed of the defined values of the spread of edge to reach points (k+1,i 2 ). I.e only the edges with effective fuel rate are considered. "Impossible" edges due to system capabilities and limitations are ignored. The key point is then to build the cost plus cost to go matrix (C+J * k+1 ) in an efficient manner (algo Fig. 11) and then take the minimums of each lines. Two index matrix index W and index C are built.

Index W (Fig. 10) contains on each line the indexes of the previous points at column k reaching a point i 2 of column k+1. Each line starts at the first defined value and is completed by non defined values if needed.

Battery SOC SOC max

k k+1 i 2 i 1_sup (i 2 ) i lim_min (k+1) i lim_max (k+1) i 1_inf (i 2 ) i 1_sup (i 2 ) i 2 SOC min 1 M ∞ ∞ J * k +C
The matrix index C contains on each line the indexes in the spread of cost which reach point i 2 . Note that the spread of cost is in fact the column k of a cost matrix C m containing all costs of the graph.

Once these two indexes matrix are built (using a loop) the corresponding cost C m (indexC) and previous best fuel consumption J * k (index W ) are added. The minimum of the lines of the resulting matrix gives the best fuel consumption at instant k J * (k,1:M) and the best trajectories index (Fig. 11). Such a method allows working with small dense matrix which leads to a computer effort gain in building matrix and matrix operations.

Note that an additional gain in computer effort can be performed building the indexes matrix column by column and not line by line. Thus the loop to construct it is smaller (S max iteration instead of M), typically 20 operations instead of 1000. For simplicity reasons the algorithm is presented line by line but the computer effort is indicated using the column by column method. Three methods to solve the DDP problems in the case of HEV have been developed using Matlab Software: a first intuitive one, an intuitive matrix method and a method using dense matrix. To provide Battery SOC SOC max

end index index index index C index J index J end i i i i index i i i i index i i i for C k k for W M k C m W k M k i i i s W i s i s i s C k k ) ( index) ing correspond the (determine )] ( ) (
k k+1 i 2 i 1_sup (i 2 ) i lim_min (k+1) i lim_max (k+1) i 1_inf (i 2 ) i 1_sup (i 2 ) i 2 SOC min 1 M J* k +C ∞ S max 1 i 1_inf (i 2 ) i s_inf (i 2 ) S max 1 i 1_sup (i 2 ) 1 i s_sup (i 2 ) M S max
effective comparison of the computer effort, the example of the well known Toyota PRIUS is presented in the following of the paper.

Application on the Toyota PRIUS-II

The DDP have been applied to different types of hybrid architectures [START_REF] Koot | Energy management strategies for vehicular electric power systems[END_REF][16] [START_REF] Perez | Optimization of power management in an hybrid electric vehicle using dynamic programming[END_REF] including different SP_HEV architectures [START_REF] Vinot | Improvement of an EVT-based HEV using dynamic programming[END_REF]. In the scope of this paper, to provide data on computer effort during the process of optimisation, the well known Toyota PRIUS is chosen as an example.

The Toyota PRIUS II on sale since 2004 has been tested in our labs. Then the car and component characteristic have been inserted in our library of hybrid vehicle simulation [START_REF] Jeanneret | New Hybrid concept simulation tools, evaluation on the Toyota Prius car[END_REF][23] and presented a good set of parameters in the scope of this paper.

In the following, the Toyota Hybrid system architecture used in the PRIUS and the manner to calculate the cost is briefly remind. Thus the computer effort is detailed for each previous method and for different graph sampling.

A. THS architecture

The considered topology is that presented in Fig. 12. It is composed of one engine (ICE) and two electrical machines, one called EM1 and the other called EM2 fed by two electrical inverters (INV1 and INV 2) linked to a battery. A final axle and gear (FG) are linked to the wheels. A planetary gear (PG) enables to split the ICE energy mechanically to the electrical machine EM1 and to the wheels. The ICE is linked to the Planet carrier, the rotor of the electrical machine EM2 to the sun gear. The shaft of the FG ix connected both to the ring of the PG and to the rotor of EM1.

B. Edges costs calculation

In DDP formulation, the cost of the edges, i.e. the fuel rate is performed using a backward approach Fig. 13, [START_REF] Vinot | Improvement of an EVT-based HEV using dynamic programming[END_REF]. This is based on the knowledge of the driving cycle and the vehicle characteristics (inertia and resistant forces). Then one can easily calculate the wheel torque T r with the vehicle speed (wheel speed, ω r ) and acceleration, and go upstream from the wheels to battery and engine.

The vehicle model is an energetic model based on the VEHLIB library [START_REF] Jeanneret | New Hybrid concept simulation tools, evaluation on the Toyota Prius car[END_REF], [START_REF] Trigui | Global Forward-Bacward Approach for a Systematic Analysis and Implementation[END_REF]. The electrical machines are modelled with efficiency maps and the engine with fuel consumption maps. The battery model is a resistive model with parameters dependant on SOC and current. From the wheel torque and speed constraints (T r ,ω r ), one can easily deduce:

f r veh r T dt d J T + = ω (5)
Where T r is the wheel torque, T f is the load torque calculated from the resistant forces, J veh is the overall inertia of the vehicle brought back to the wheels and ω r is the wheel speed. Where T ring is the torque on the planetary gear ring, T em1 is the torque on the EM1 shaft, η is the efficiency of the transmission (axle plus gear), R trans is the transmission ratio, Sign(T r ,ω r ) is the sign of the wheel power. Thus the wheel torque is divided or multiplied by the gear efficiency depending on the direction of the wheel power. Note that the electrical machine and planetary ring revolution is the same.

The Willis relation in the planetary gear imposes: Where ω em2 , ω ring and ω ice are respectively the speeds of the electrical machine EM2 (sun gear), the planetary ring (electrical machine EM1) and the engine (planet carrier). k is the planetary gear ratio (k=-N teeth _ ring /N teeth_sun ).

The torque relations on the planetary gear give:

) 1 ( 2 k T T ice em - - = (8) ) ( 2 k T T plan ring em η - = (9)
Where T ice is the ICE torque and η plan is the efficiency of the planetary gear.

O the electrical node between the electrical machines and the battery: [START_REF] Romaus | Optimal energy management for a hybrid energy storage system for electric vehicles based on stochastic dynamic programming[END_REF] Where I em1 , I em2 and I bat are respectively the DC current in electrical inverters of the EM1 and EM2 electrical machines (considered as receptors) and in battery (considered as a generator).

And finally the relations on the electrical machines impose:

) , ( ) ( 1 1 1 1 1 1 em em em em em em bat T Pt T I RI E ω ω + = - (11) ) , ( ) ( 2 2 2 2 1 2 em em em em em em bat T Pt T I RI E ω ω + = - (12)
Where Pt em1 and Pt em2 are the losses in the electrical machines and inverters, E is the open circuit voltage of the battey and R is the internal battery resistance which depends only on SOC. Note that ω em1 is the speed of the common shaft to EM1 and planetary ring and is known using final gear relation:

r trans ring em R ω ω ω = = 1 ( 13)
Considering ( 4) to [START_REF] Romaus | Optimal energy management for a hybrid energy storage system for electric vehicles based on stochastic dynamic programming[END_REF], and I bat fixed (on an edge) we obtain 7 equations and 8 unknown variables (ω em2 , ω ice , T em2 , T ice , T ring , T em1 , I em2 , I em1 ). That means one degree of freedom to be fixed. This degree of freedom allows choosing, for example, the engine speed to minimize the fuel consumption for a given edge (I bat , T r and ω r imposed), [START_REF] Vinot | Improvement of an EVT-based HEV using dynamic programming[END_REF].

The process of cost calculation is totally independent on the solving of the graph. And is performed once and for all edge of the graph before the DDP process. The computer effort for the calculation of the edges costs and graph solving is then easy to determine separately.

C. Computer effort

This part presents the computer effort to solve DDP for the three previously presented methods. The computer effort is indicated in second and is calculated using the Matlab instruction "cputime". The simulations are performed on cpu with the characteristic presented Table 1. The bogomips calculated by Linux give a good indication of the potential million of instructions per second using this processor. This can allow a comparison with other types of processors.

The computer effort is presented Table 3 and Table 4 for different driving cycles representing urban, road and highway driving condition. The NEDC normalized European driving cycle is also presented.

It is noted that the memory used with the different methods remains small (less than 2Go) and well under the capacity of the cpu.

Two different graph precisions are used. It represents the SOC difference in % between two consecutive points in a column. The time sample is fixed to 1s to each cycle. Table 3 and Table 4 shows that the third method is the most efficient to solve DDP applied to the hybrid vehicle. In fact the time to solve the DDP can be divided by a factor two to six depending on the cycle and the precision. A method using "full" matrix (method B), can improve the time calculation compared to point to point method (method A) if the graph precision is not too high but became definitely too expensive in computer effort as the size of the matrix (graph precision) increases. It is also noted that an increase of the graph precision below 0.01% does not change significantly the fuel consumption. Using such precision, the total time calculation of DDP (edge cost calculation and DDP solving) does not exceed 10 s and can be reduced to 2 s in small cycle. It is then possible to use DDP as a tool to determine fuel consumption in more a complex process [START_REF] Reinbold | Global optimization of a parallel hybrid vehicle using optimal energy management[END_REF].

Conclusion

Different methods to solve dynamic programming applied to hybrid vehicle have been studied; an intuitive method using a "point to point" approach, a matrix based method easy to implement but using full sparse matrix and a method based on dense matrix.

The last method takes advantage of one particularity of DDP applied to the hybrid vehicle: the spread of possible edge from one point of a column to the next column is highly limited by the system capacities. The number of edges from one point is thus small compared to the number of points in a column.

The last method is then the most efficient in terms of computer effort. In fact it allows to use small index matrix (compared to method 2), and avoid to call a minimum function in loops (column loops inserted in time loops).

The tool developed to calculate the fuel consumption by means of DDP has been used and is currently used in more complex studies. With a global time calculation less than 30s (usually less than 5s) it can be used in comparative hybrid architecture studies and even in global optimal design architecture process.

Fig. 1

 1 Fig. 1 DDP applied to Hybrid vehicle

FigFig. 5

 5 Fig. 4 Regular Grid

  Fig. 6 Point to point method

Fig. 9

 9 Fig. 9 Algorithm use in a matrix method

Fig. 10

 10 Fig. 10 index matrix method principle

Fig. 11

 11 Fig. 11 index matrix method algorithm

Fig. 12

 12 Fig. 12 THS architecture

Fig. 13 :

 13 Fig. 13: Backward approach of THS architecture On each edge of the graph, I bat , T r and ω r are known. The equation of each component of the system has then to be used to calculate the engine operating point and thus its fuel consumption.

Fig. 14

 14 Fig. 14 Planetary gear connexions Considering Fig. 13 and Fig. 14 the torque relation on the ring and EM1 shaft is:

Table 1

 1 Cpu characteristic

	Model name	Intel Xeon X5650
	Cpu MHz	2660
	Cache size	12.3 KB
	Memory	24 GB
	bogomips	5320
	Table 2 Graph information
	SOC max	80 %
	SOC min	20 %
	Speed engine sample 50.π/30 rad/s
	Torque engine sample 5 Nm

Table 3

 3 Computer effort for different driving cycle with a graph precision of 0.05%

Table 4

 4 Computer effort for different driving cycle with a graph precision of 0.01%

		URBAN	ROAD	HIGHWAY	NEDC
	graph precision
	time of cycle in s	560	843	1805	1181
	million of edge	3.5	6.3	18	7.9
	time of edge cost calculation in s	1	1.9	5.1	2.5
	fuel consumption in l/100 km	3.772	4.124	5.309	3.720
	A: point to point method	3.8	7.2	18.7	10.4
	B: Matrix method	3.0	4.9	13.5	6.9
	C: Index matrix method	0.6	1.2	4.6	1.9

drive cycle -0.05 time to solve DP in s


João Manoel Gomes Da Silva

Isabelle Queinnec

Alexandre Seuret

Sophie Tarbouriech

Regional stability analysis of discrete-time dynamic output feedback under aperiodic sampling and input saturation*

Keywords: Dynamic output feedback, saturating input, sampled-data systems, aperiodic sampling

This paper addresses the problem of regional stability analysis of a continuous-time plant controlled by a discrete-time dynamic output feedback control law with saturation constraints under aperiodic sampling. To cope with the aperiodic sampling problem, the proposed approach is based on an impulsive system modeling and the use of looped functionals. A generalized sector-based relation is applied to tackle the control saturation effects. From these ingredients, conditions to ensure regional asymptotic stability of the closed-loop system under aperiodic sampling are derived. Based on these conditions, LMI-based optimization problems are proposed to compute estimates of the region of attraction of the closed-loop system or to maximize the bound on the maximal admissible interval between two successive sampling instants, for which the regional stability can be ensured with respect to a given set of admissible states. A numerical example illustrates the application of the results.

I. INTRODUCTION

Issues regarding the implementation of sampled-data systems and the effects of sampling on the stability of closed-loop systems have been intensively studied since the 80's (see, for instance, [START_REF] Aström | Computer Controlled Systems. Theory and Design. Information and System Sciences Series[END_REF], [START_REF] Chen | Optimal sampled-data control systems[END_REF], [START_REF] Chung | Sampled-data observer error linearization[END_REF], [START_REF] Grizzle | Feedback linearization of sampled-data systems[END_REF], [START_REF] Nešić | Formulas relating KL stability estimates of discrete-time and sampled-data nonlinear systems[END_REF], [START_REF] Nešić | A framework for stabilization of nonlinear sampled-data systems based on their approximate discrete-time models[END_REF]). More recently, mainly motivated by networked control applications [START_REF] Zhang | Stability of Networked Control Systems[END_REF], stability of aperiodic sampled-data systems has benefited from a second wave of theoretical developments. The results in this context allow to account for aperiodic samplings in a formal way and basically rely upon four efficient approaches of modeling. The first one is based on an uncertain discrete-time model, which is embedded onto a polytopic model from the use of exponential matrices [START_REF] Cloosterman | Controller synthesis for networked control systems[END_REF], [START_REF] Cloosterman | Stability of networked control systems with uncertain timevarying delays[END_REF]. A second approach refers to a robust analysis of continuous or discrete-time systems. The effects of the aperiodic sampling is embedded into norm bounded uncertainties and the proposed results are based on region-dividing techniques of the admissible inter sampling interval [START_REF] Fujioka | Stability analysis of systems with aperiodic sample-andhold devices[END_REF] or on IQCs and passivity arguments [START_REF] Kao | On stability of systems with aperiodic sampling devices[END_REF], [START_REF] Oishi | Stability and stabilization of aperiodic sampleddata control systems using robust linear matrix inequalities[END_REF]. Another approach regards the modeling of the sampling effects considering a particular time-varying delay on the control signal. This approach, initially proposed in [START_REF] Fridman | Robust sampled-data stabilization of linear systems: An input delay approach[END_REF] and further improved in [START_REF] Fridman | A refined input delay approach to sampled-data control[END_REF] or [START_REF] Naghshtabrizi | Exponential stability of impulsive systems with application to uncertain sampled-data systems[END_REF] is based on the use of Lyapunov-Krasovskii functionals and impulsive systems representation for systems with time-varying delays. Similarly to [START_REF] Sivashankar | Characterization of the l 2 -induced norm for linear systems with jumps with applications to sampled-data systems[END_REF], a mixed continuous-and discrete-time approach has been proposed in [START_REF] Seuret | A novel stability analysis of linear systems under asynchronous samplings[END_REF], [START_REF] Seuret | Taking into account period variations and actuator saturation in sampled-data systems[END_REF] through the concept of looped-functionals and adapted to the case of impulsive systems in [START_REF] Briat | A looped-functional approach for robust stability analysis of linear impulsive systems[END_REF]. This alternative method authorizes a larger class of functionals than in the Lyapunov-Krasovskii approach and has some similarities with the notion of clock-dependent Lyapunov-function, employed for instance in [START_REF] Briat | Convex conditions for robust stability analysis and stabilization of linear aperiodic impulsive and sampled-data systems under dwell-time constraints[END_REF]. The reader may refer to [START_REF] Briat | Theoretical and numerical comparisons of looped functionals and clock-dependent Lyapunov functions -The case of periodic and pseudo-periodic systems with impulses[END_REF] for detailed comparison between these two approaches. It should however be highlighted that the aforementioned methods are restricted to state feedback control laws. *This work was supported by ANR-France projects LimICoS (contract 12 BS03 005 01) and SCIDiS (contract 15-CE23-0014), STICAmSud project ADNEC, CNPq-Brazil (grants Univ. 480638/2012-8 and PQ 305979/2015-9) and CAPES-Brazil (projects PVE-129/2012 and STICAmSud) J.M. Gomes da Silva Jr. is with UFRGS, EE, Department of Automation and Energy (DELAE), Porto Alegre, Brazil jmgomes@ece.ufrgs.br I. Queinnec, A. Seuret and S. Tarbouriech are with CNRS, LAAS, 7 av. du Colonel Roche, F-31400 Toulouse, France and with Univ de Toulouse, LAAS, F-31400 Toulouse, France queinnec,aseuret,tarbour@laas.fr On the other hand, a large attention has been paid to the stability analysis and stabilization of systems with input saturation (see, for example, [START_REF] Hu | Control systems with actuator saturation: analysis and design[END_REF], [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF]). In this context, many methods guaranteeing global or regional (local) asymptotic stability of the closed-loop systems under saturating control laws have been proposed, either in the continuous or discrete-time frameworks. Considering the periodic sampling case, the problem of assessing stability of a continuoustime plant controlled by a discrete-time dynamic output feedback control law has been addressed for instance, in [START_REF] Dai | Output feedback synthesis for sampled-data system with input saturation[END_REF]. In [START_REF] Herrmann | Discrete-time and sampled-data anti-windup synthesis: stability and performance[END_REF], the authors extended the design of discrete-time anti-windup loops to the sampled-data case with a constant period by using the approach proposed in [START_REF] Chen | Optimal sampled-data control systems[END_REF], in order to transpose the problem into a discrete-time framework.

In the present paper, we are interested in analyzing stability of a continuous-time plant controlled by a discrete-time dynamic output feedback control law with saturation constraints under aperiodic sampling. In this context, to cope with the aperiodic sampling problem, the proposed approach is based on the impulsive system modeling as suggested in [START_REF] Goebel | Hybrid Dynamical Systems. Modeling, Stability and Robustness[END_REF] and on the use of looped functionals. Furthermore, a generalized sector-based condition is applied to tackle the control saturation effects. From these ingredients, conditions to ensure regional asymptotic stability of the closed-loop system under aperiodic sampling are derived. Based on these conditions, LMIbased optimization problems are proposed to compute estimates of the region of attraction of the closed-loop system or to maximize the bound on the maximal admissible interval between two successive sampling instants, for which regional stability can be ensured for a given set of admissible states.

Notation. S n denote the set symmetric matrices of R n×n . For given positive scalars 0 < T 1 ≤ T 2 , K is defined as the set of continuous functions from an interval [0, T ] to R n , where T is a positive scalar in [T 1 , T 2 ]. | • | and • stand for the absolute value and the Euclidean norm, respectively. The notation P > 0 for P ∈ S n means that P is a symmetric positive definite matrix. For any positive integer j ≤ n, any vector x ∈ R n and any matrix A ∈ R n×m , A j , x j refer to the jth line of matrix A, the jth component of vector x, respectively. For A ∈ R n×m , He{A} = A + A ′ . I and 0 represent the identity and zero matrices of appropriate dimensions. Co{•} denotes a convex hull. For v ∈ R m , each component of the vector-valued saturation function, sat(v), is defined by sat(v) i = sign(v i ) min{u 0i , |v i |}, i = 1, . . . , m, where ±u 0i denote the symmetric limits on the ith component of v.

II. PROBLEM FORMULATION

A. Plant model

Consider the continuous-time linear plant with saturating inputs described by:

ẋp(t) = Apxp(t) + Bpsat(ū(t)), y(t) = Cpxp(t), (1)
where xp ∈ R np , ū ∈ R m and y ∈ R p represent the state, the input and the output vectors of the plant, respectively. Matrices Ap, Bp, Cp have appropriate dimensions and are supposed to be constant. We assume that the control inputs are sampled, possibly in an aperiodic manner, and the value of ū is kept constant (through a zero-order hold) between two successive sampling instants i.e.:

ū(t) = u(t k ), ∀t ∈ [t k , t k+1 ), (2)
The sequence of increasing positive sampling instants is denoted by

T := {t k } k≥0 . It is assumed that T := {t k } k≥0 is such that k∈N [t k , t k+1 ) = [0, + ∞)
, so that no Zeno phenomena occurs. Moreover, we assume that there exist two positive scalars T1 ≤ T2 such that the difference between two successive sampling instants

T k = t k+1 -t k satisfies 0 < T1 ≤ T k ≤ T2, ∀k ∈ N. (3)
Note that the inter sampling time T k can vary with time, which allows to model an aperiodic sampling policy. The particular case of periodic sampling corresponds to T k = T1 = T2 = δ, ∀k ∈ N.

For output feedback purposes, we consider that the output y(t) is sampled at the same instants of the control input, i.e. at each sampling instant t k , the sample y(t k ) is generated.

B. Discrete-time dynamic output feedback controller

Considering a digital implementation, we assume that system (1) is controlled by a linear discrete-time dynamic output feedback controller, described in state space as follows:

xc(t k+1 ) = Acxc(t k ) + Bcy(t k ) + Ecψ(u(t k )), u(t k ) = Ccxc(t k ) + Dcy(t k ), (4)
where xc ∈ R nc , y ∈ R p and u ∈ R m are the state, the input and the output of the controller, respectively. The matrices Ac, Bc, Cc, Dc, Ec are assumed to be constant and of appropriate dimensions. ψ(u) is a vector-valued decentralized deadzone nonlinearity defined as follows:

ψ(u) = sat(u) -u. (5)
The term Ecψ(u) regards a static anti-windup compensation, which can be appended to the controller to mitigate the saturation effects on performance and stability [START_REF] Gomes Da | Anti-windup design with guaranteed region of stability: an LMI-based approach[END_REF], [START_REF] Gomes Da | Anti-windup design with guaranteed regions of stability for discrete-time linear systems[END_REF].

In this case, at each sampling instant, the plant output is sampled and its value is used to instantaneously update both the controller state and the control input to be applied to the plant between in the interval [t k t k+1 ).

Remark 1: It should be noticed that the method used to design the controller (4) is out of the scope of the present work. Our main goal is to provide a method to assess the stability of the sampled-data closed-loop system defined by the connection between this controller and the continuous plant [START_REF] Aström | Computer Controlled Systems. Theory and Design. Information and System Sciences Series[END_REF]. On the other hand, it is reasonable (and relevant from a practical point of view) to assume that (4) has been designed considering classical techniques based on a periodic sampling paradigm, i.e. with T k = δ, ∀k, and then to analyze the effects on stability of a possibly aperiodic sampling policy with T k ∈ [T1, T2]. In this case, we can assume for instance that matrices Ac, Bc, Cc, Dc, Ec, have been designed considering one of the two approaches below:

1) Discrete-time design. First, an exact discretization of the plant should be performed. Considering t k+1 -t k = δ, ∀k, this leads to the following discrete-time model:

xp(t k+1 ) = A p,d xp(t k ) + B p,d sat(u(t k )), y(t k ) = C p,d xp(t k ) (6)
with

A p,d = e Apδ , C p,d = Cp, B p,d =
δ 0 e Ap(δ-s) dsBp. Thus, based on model [START_REF] Chung | Sampled-data observer error linearization[END_REF], known discrete-time results can be applied to design matrices Ac, Bc, Cc and Ec in order to stabilize (considering periodic sampling) the closed-loop system.

2) Controller discretization. First, a continuous-time stabilizing controller, described by

ẋc(t) = Ācxc(t) + Bcy(t) + Ēcψ(u(t)), u(t) = Ccxc(t) + Dcy(t), (7)
is designed considering the continuous-time plant model [START_REF] Aström | Computer Controlled Systems. Theory and Design. Information and System Sciences Series[END_REF].

Then, assuming a periodic sampling, a discretization of this controller can be performed using classical Euler or Tustin approximations. Alternatively, a discrete-time approximation ensuring that the value of the output of the continuous and the discrete-time controllers coincide at the sampling instants t k = kδ can be considered. This leads to a discrete-time controller as in (4) with:

Ac = e Ācδ , Cc = Cc, Dc = Dc Bc = δ 0 e Āc(δ-s) ds Bc, Ec = δ 0 e Āc(δ-s) ds Ēc (8)

C. Hybrid modeling of the closed-loop system

Considering the controller (4), the state of the plant evolves in continuous-time while the controller state and the plant input have discrete-time updates. These behaviors can be cast in the hybrid dynamical systems framework [START_REF] Goebel | Hybrid Dynamical Systems. Modeling, Stability and Robustness[END_REF]. From the previous assumptions, the closed-loop system dynamics can be represented by an impulsive system described by:

     ẋp = Apxp + BpDcCpxm + BpCcxc + Bpψ, ẋm = 0, ẋc = 0, ∀t ∈ R\T,      x + p = xp, x + m = xp, x + c = Acxc + BcCpxm + Ecψ, ∀t ∈ T, (9)
where x + = x(t k+1 ) and ψ is the deadzone function given by

ψ = sat(Ccxc + DcCpxm) -(Ccxc + DcCpxm)
When t ∈ R\T, the evolution of the plant state follows the first differential equation. The two other variables, representing the held values of the plant state and the controller state, remain constant in the interval [t k , t k+1 ). Since the plant dynamics evolves continuously, if t ∈ T the plant state remains at its current value (i.e. x + p = xp), while the other variables are impulsively reset. We introduced the auxiliary variable xm in [START_REF] Dai | Output feedback synthesis for sampled-data system with input saturation[END_REF] to represent the memory of the current sampled value of the plant state. Actually, it couples the continuous and the discrete-time dynamics through the control signal. Note that

u(t k ) = Ccxc(t k ) + Dcy(t k ) = Ccxc(t k ) + DcCpxm(t k ).

D. Stability Analysis Problems

Considering the nonlinear closed-loop system (9) we are interested in analyzing both the influence of the aperiodic sampling and the control saturation on its asymptotic stability. Since no particular assumption is assumed regarding the stability of the continuous-time plant, i.e. matrix Ap is not necessarily Hurwitz, we mainly focus on the regional stability analysis. Actually, it should be recalled that if the open-loop system is exponentially unstable, only regional stability can be ensured in the presence of control saturation [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF]. In this case, the region of attraction of the origin is defined as the set of all initial states for which the corresponding trajectories converge asymptotically to the origin. Nevertheless, the exact characterization of this region is in general not possible. Hence, it is of major interest the determination of analytical estimates of the region of attraction, which can be seen as regions of asymptotic stability (RAS) [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF]. Hence, we are concerned by the following problems: P1. Given T1 and T2, such that 0 < T1 ≤ T2, determine an estimate of the region of attraction of system (9). P2. Given a set of admissible initial states X0, maximize the length of the interval [T1, T2] of admissible inter sampling times T k , for which asymptotic stability of system ( 9) is ensured for any initial condition in X0.

III. PRELIMINARIES

A. Looped-functional approach

In order to tackle problems P1 and P2, instead of considering the classical hybrid framework to study mixed continuous and discrete dynamics as defined in [START_REF] Goebel | Hybrid Dynamical Systems. Modeling, Stability and Robustness[END_REF], we use an alternative direction by adapting the technique based on a looped functional as developed in [START_REF] Briat | A looped-functional approach for robust stability analysis of linear impulsive systems[END_REF]. This approach relies on the characterization of the trajectories of system (9) in a lifted domain [START_REF] Yamamoto | New approach to sampled-data control systems -a function space method[END_REF], [START_REF] Briat | A looped-functional approach for robust stability analysis of linear impulsive systems[END_REF], [START_REF] Seuret | A novel stability analysis of linear systems under asynchronous samplings[END_REF]. The idea is to view the entire state-trajectory as a sequence of lifted state function

χ k (τ ) ∈ K, χ k : [0, T k ] → ℜ n , with [0, T k ] ⊂ [0,
T2] and the vector χ k (τ ) ∈ R n , with n = 2np + nc, being defined as follows:

χ k (τ ) =   χ p,k (τ ) χ m,k (τ ) χ c,k (τ )   =   xp(t k + τ ) xm(t k + τ ) xc(t k + τ )   with χ k (0) = lim s→t + k   xp(s) xm(s) xc(s)   . (10)
Then, by using [START_REF] Fridman | A refined input delay approach to sampled-data control[END_REF], the hybrid closed-loop system (9) can be alternatively represented as follows:

                                   χp,k (τ ) = Apχ p,k (τ ) + BpDcCpχ m,k (τ ) +BpCcχ c,k (τ ) + Bpψ k (τ ), χm,k (τ ) = 0, χc,k (τ ) = 0, χ p,k+1 (0) = χ p,k (T k ), χ m,k+1 (0) = χ p,k (T k ), χ c,k+1 (0) = Acχ c,k (T k ) + BcCpχ m,k (T k ) +Ecψ k (T k ), (11)
Note that, since χm,k (τ ) = 0 and χc,k (τ ) = 0, it follows that ψk (τ ) = 0, i.e.

ψ k (τ ) = ψ k (0) = ψ(Ccχ c,k (0) + DcCpχ m,k (0)), for all τ ∈ [0, T k ]. Moreover χ m,k (T k ) = χ m,k (0) and χ c,k (T k ) = χ c,k (0).
Based on the lifted function χ k , we recall the following result, which is the basis of the looped-functional approach to assess asymptotic stability of system [START_REF] Dai | Output feedback synthesis for sampled-data system with input saturation[END_REF].

Lemma 1: [START_REF] Briat | A looped-functional approach for robust stability analysis of linear impulsive systems[END_REF] Consider system (9) and let 0 < T1 ≤ T2, be two positive scalars and V : R n → R + be a function for which there exist real scalars 0 < µ1 < µ2 such that

∀χ ∈ R n , µ1||χ|| 2 ≤ V (χ) ≤ µ2||χ|| 2 . ( 12
)
Then, the two following statements are equivalent.

(i) For all k ∈ N, T k ∈ [T1, T2], the forward difference of the function V satisfies ∆V (k) := V (χ k+1 (0)) -V (χ k (0)) < 0;
(ii) There exists a continuous functional V0 : [0, T2] × K → R, which satisfies, for all z ∈ K and

∀T k ∈ [T1, T2], V0(T k , z) = V0(0, z) (13)
and such that, for all k ∈ N,

T k ∈ [T1, T2] and τ ∈ [0, T k ], Ẇ(τ, χ k ) = Λ k T k + d dτ [V (χ k (τ )) + V0(τ, χ k )] < 0, (14)
with

Λ k = V (χ k+1 (0)) -V (χ k (T k )).
Moreover, if one of these two statements holds, then system (9) is asymptotically stable for all T k ∈ [T1, T2] .

Note that from this result it suffices to analyze the behavior of the system in the interval [0, T k ] to assess asymptotic stability of the closed-loop system [START_REF] Dai | Output feedback synthesis for sampled-data system with input saturation[END_REF]. Moreover, differently from other approaches based on Lyapunov functionals, an important feature of the present approach is that V0 is not required to be positive.

B. Generalized sector condition

Based on the variables used to compute the control signal, let us

define χk = χ m,k χ c,k
. From this definition it follows that ψ k (0) =

ψ k (K χk (0)) with K = [DcCp Cc].
Thus, the effects of the saturation nonlinearity can be considered through a generalized sector condition, that applies to deadzone nonlinearities, which is summarized in the next lemma.

Lemma 2: [START_REF] Gomes Da | Anti-windup design with guaranteed region of stability: an LMI-based approach[END_REF], [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF] Consider K, G ∈ R m×(np+nc) and define the set S = { χ ∈ R np +nc ; |(Kj -Gj ) χ| ≤ u0j , j = 1, . . . , m}. If χ ∈ S, then the deadzone nonlinearity ψ(K χ), satisfies the following condition:

ψ(K χ) ′ U (ψ(K χ) + G χ) ≤ 0, (15)
for any diagonal positive definite matrix U ∈ R m×m .

IV. MAIN RESULT

Based on the two previous Lemmas, the following theorem is stated. It provides conditions to determine an ellipsoidal RAS for system (9) under aperiodic sampling.

Theorem 1: For given positive scalars 0 < T1 ≤ T2, assume that there exist symmetric positive definite matrices P, R, symmetric matrices X, S1, a diagonal positive definite matrix U , matrices G, V , S2 and N of appropriate dimensions and a positive scalar α, that satisfy, for i = 1, 2, j = 1, • • • , m

Ψ1(Ti) = Π1 + Ti(Π3 + Π4) TiN TiN ′ -TiR < 0, (16)
Ψ2(Ti) = Π1 + Ti(Π2 -Π3 + Π4) < 0,

P + V M f + M ′ f V ′ [0 Kj -Gj] ′ [0 Kj -Gj ] αu 2 0j > 0, (17)
with

Π1 = M ′ + P M+ -M ′ -P M-+ Π10 -He M ′ 5 U M5 + M ′ 5 U G M3 M4 , Π10 = M ′ 12 S1M12 + He{M ′ 12 S2M2 + N M12}, Π2 = M ′ c RMc + He{M ′ c (S1M12 + S2M2)}, Π3 = M2 M5 ′ X M2 M5 , Π4 = He      Mc 0 0   ′ P   M1 M3 M4      , (19)
where

M1 = I 0 0 0 0 , M2 = 0 I 0 0 0 , M3 = 0 0 I 0 0 , M4 = 0 0 0 I 0 , M5 = 0 0 0 0 I , M12 = M1 -M2, M+ =   0 I 0 0 0 0 I 0 0 0 0 0 BcCp Ac Ec   , M-=   0 I 0 0 0 0 0 I 0 0 0 0 0 I 0   , M f = I -I 0 , Mc = Ap 0 BpDcCp BpCc Bp , (20)
Then, if χ0(0) ∈ E (P, α) = χ ∈ ℜ n ; χ ′ P χ ≤ α -1 , then for any aperiodic sampling satisfying (3), χ k (0) → 0 as k → ∞, i.e. E (P, α) is a RAS for the closed-loop system [START_REF] Dai | Output feedback synthesis for sampled-data system with input saturation[END_REF].

Proof. Let a quadratic Lyapunov function candidate be defined, for any χ in R n , by V (χ) = χ ′ P χ, where P ∈ S n , P > 0. Thus, the function V satisfies [START_REF] Fujioka | Stability analysis of systems with aperiodic sample-andhold devices[END_REF]. Consider now a functional V0 defined for any τ ∈ [0, T k ], as follows:

V0(τ, χ k ) = τ T k (χ p,k (τ )-χ p,k (T k )) ′ S1(χ p,k (τ )-χ p,k (T k )) + 2 τ T k (χ p,k (τ ) -χ p,k (T k )) ′ S2χ p,k (T k ) + (T k -τ ) T k τ χ p,k (T k ) ψ k (0) ′ X χ p,k (T k ) ψ k (0) -τ T k T k τ χ′ p,k (θ)R χp,k (θ)dθ, with S1 ∈ S np , S2 ∈ R np×np , X ∈ S np+m , R ∈ S np such that R > 0. Since V0(τ, χ k
) is equal to zero at τ = 0 and τ = T k , it satisfies condition [START_REF] Goebel | Hybrid Dynamical Systems. Modeling, Stability and Robustness[END_REF]. Moreover, it is continuous at all sampling instants and differentiable over [0, T k ).

Then, from Lemma 1, if we prove that W(τ, χ k ) defined from V (χ k ) and V0(τ, χ k ) above, is such that Ẇ(τ, χ k ) < 0 along the trajectories of system (9), we can conclude that ∆V (k) < 0, which ensures that the trajectories converge asymptotically to the origin.

From Lemma 2, if χk (0) ∈ S, it follows that

Ws(χ k , ψ k ) := 2ψ ′ k (0)U ψ k (0) + G χ m,k (0) χ c,k (0) ≤ 0, ( 21
)
for any diagonal positive definite matrix U . Therefore we can combine ( 21) and ( 14) to relax the condition Ẇ(τ, χ k ) < 0 in a regional context. With this aim, we define

Ẇg(τ, χ k ) := Λ k T k + V (χ k (τ )) + V0(τ, χ k ) - Ws(χ k , ψ k ) T k . (22)
Then, if Ẇg(τ, χ k ) < 0, it follows that Ẇ(τ, χ k ) < 0 provided that χk (0) ∈ S. In order to obtain an upper bound on the functional Ẇg(τ, χ k ), we define the extended vector

ζ k (τ ) = [χ p,k (τ ) ′ χ p,k (T k ) ′ χ m,k (0) ′ χ c,k (0) ′ ψ k (0) ′ ] ′ ,
and we compute each term of Ẇg(τ, χ k ) separately. Let us first focus on

Λ k = V (χ k+1 (0)) -V (χ k (T k )).
From the definition of V (χ k ) and matrices M+ and Min [START_REF] Kao | On stability of systems with aperiodic sampling devices[END_REF], it follows that

Λ k = ζ k (τ ) ′ [M ′ + P M+ -M ′ -P M-] ζ k (τ ). (23)
To compute the second term of Ẇg, V (χ k (τ )), note that χ m,k (τ ) = χ m,k (0), χ c,k (τ ) = χ c,k (0) and ψ k (τ ) = ψ k (0). Hence, considering the auxiliary matrices defined in [START_REF] Kao | On stability of systems with aperiodic sampling devices[END_REF], we get

dV (χ k (τ )) dτ = ζ k (τ ) ′   He      Mc 0 0   ′ P   M1 M3 M4        ζ k (τ ).
Similarly, the sector condition, Ws, can be re-written as follows:

Ws = ζ k (τ ) ′ He M ′ 5 U M5 + G M3 M4 ζ k (τ ). ( 24
)
Let us finally consider the third component of Ẇg, i.e. V0(τ, χ k ). Considering the derivative of V0 with respect to τ and using the matrices Π10, Π2 and Π3 defined in [START_REF] Hu | Control systems with actuator saturation: analysis and design[END_REF], it follows that

V0(τ, χ k ) ≤ 1 T k ζ k (τ ) ′ [Π10 + τ Π2 + (T k -2τ )Π3 +(T k -τ )N R -1 N ′ ζ k (τ ), (25)
where the upper bound has been obtained by noting that, for any matrix N in R (3np +nc+m)×np and any symmetric positive definite matrix R ∈ R np×np , the following inequality holds [START_REF] Seuret | Taking into account period variations and actuator saturation in sampled-data systems[END_REF]:

-

T k τ χp,k (θ) ′ R χp,k (θ)dθ ≤ 2ζ k (τ ) ′ N (χ p,k (τ ) -χ p,k (T k )) +(T k -τ )ζ k (τ ) ′ N R -1 N ′ ζ k (τ ).
Hence, combining the previous expressions, it follows that

Ẇg(τ, χ k ) ≤ 1 T k ζ k (τ ) ′ [Π1 + τ Π2 + T k Π4 +(T k -τ )N R -1 N ′ + (T k -2τ )Π3]ζ k (τ ), (26)
with Π1, Π2, Π3 and Π4 as defined in [START_REF] Hu | Control systems with actuator saturation: analysis and design[END_REF]. A sufficient condition to ensure that Ẇg < 0 consists in guaranteeing that

Π1 + τ Π2 + (T k -2τ )Π3 + T k Π4 +(T k -τ )N R -1 N ′ < 0. (27)
Note now that the inequality ( 27) is affine with respect to the variable τ in [0, T k ]. Then, by convexity, it suffices to ensure that ( 27) is negative for τ = 0 and τ = T k (see [START_REF] Naghshtabrizi | Exponential stability of impulsive systems with application to uncertain sampled-data systems[END_REF] for more details). Thus Ẇg < 0 if both inequalities

Π1 + T k (Π3 + Π4 + N R -1 N ′ ) < 0, Π1 + T k (Π2 -Π3 + Π4) < 0,
are satisfied. Applying the same argument on T k for the interval [T1, T2] and by using the Schur's complement, conditions Ψ1(Ti) < 0 and Ψ2(Ti) < 0, for i = 1, 2, given in ( 16) and ( 17), are obtained. Consider now matrix inequality [START_REF] Herrmann | Discrete-time and sampled-data anti-windup synthesis: stability and performance[END_REF]. From Schur's complement, it follows that:

χ k (0) ′ (P + V M f + M ′ f V - 0 (K j -G j ) ′ α -1 u -2 0j 0 (K j -G j ) ′ ′ )χ k (0) > 0 (28) Since χ m,k (0) = χ p,k (0)
, it follows that M f χ k (0) = 0. Thus, from [START_REF] Sivashankar | Characterization of the l 2 -induced norm for linear systems with jumps with applications to sampled-data systems[END_REF], if χ k (0) ∈ E (P, α), then χk (0) ∈ S. Hence, provided ( 16) and ( 17) are verified, it follows that Ẇ(τ, χ k ) < 0 and, from Lemma 1, it follows that the state trajectory converges asymptotically to the origin. ♦ Remark 2: Considering the particular choice of G = K, it follows that the sector condition (15) is globally valid, i.e. it is valid ∀ χk (0) ∈ R np+nc . In this case, the satisfaction of conditions ( 16) and [START_REF] Grizzle | Feedback linearization of sampled-data systems[END_REF] ensures the global asymptotic stability of the origin of the closedloop system under the aperiodic sampling (see more details in [START_REF] Gomes Da | Stability analysis of dynamic output controllers under aperiodic sampling and input saturation[END_REF]).

V. OPTIMIZATION PROBLEMS

In this section, we propose some optimization schemes based on the conditions of Theorem 1 to provide solutions to problems P1 and P2 stated in Section II.

A. Optimization of the estimate of the domain of attraction (P1)

Since the conditions of Theorem 1 guarantee that E (P, α) is a RAS, an estimate of the actual region of attraction can be obtained by maximizing E (P, α) considering some size criterion [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF]. In this case, given T1 and T2, a solution to P1 can be obtained from the following optimization problem: min trace(P ) + α subject to (16), ( 17), [START_REF] Herrmann | Discrete-time and sampled-data anti-windup synthesis: stability and performance[END_REF]

It should be noticed that relations [START_REF] Gomes Da | Anti-windup design with guaranteed regions of stability for discrete-time linear systems[END_REF] and [START_REF] Grizzle | Feedback linearization of sampled-data systems[END_REF] are nonlinear due to the product between variables G and U . However, since U is a diagonal matrix with dimension equal to the number of inputs, the solution to problem (29) can be computed by iteratively solving convex LMI problems over a grid on the elements of U . For instance, in the mono-input case U reduces to a positive scalar and a line search algorithm can be applied. Another option may be to encapsulate the LMI optimization step in some overall nonlinear optimization procedure (such as the one implemented by fminsearch Matlab function).

B. Optimization of the allowable sampling interval (P2)

The idea is to enlarge the interval [T1, T2] for which regional stability is guaranteed for a given set of admissible initial conditions X0 under an aperiodic sampling policy with T k ∈ [T1, T2]. In this case, assuming that the discrete-time controller (4) has been computed from a nominal periodic sampling period δ, we can consider T1 = δµ and T2 = δ + µ, and the following optimization problem can be formulated: max µ subject to [START_REF] Gomes Da | Anti-windup design with guaranteed regions of stability for discrete-time linear systems[END_REF], [START_REF] Grizzle | Feedback linearization of sampled-data systems[END_REF], [START_REF] Herrmann | Discrete-time and sampled-data anti-windup synthesis: stability and performance[END_REF], X0 ⊆ E (P, α)

The condition X0 ⊆ E (P, α) can be easily expressed as an LMI constraint if X0 is described either by an ellipsoid or by a polytope. Hence, this problem can be solved by iteratively increasing µ and testing LMI conditions.

VI. NUMERICAL EXAMPLE

Let us consider the model of an unstable continuous-time balancing pointer system, expressed as in [START_REF] Aström | Computer Controlled Systems. Theory and Design. Information and System Sciences Series[END_REF] with the following matrices:

Ap = -0 1 1 0 , Bp = 0 -1 , Cp = 1 1
with the saturation level u0 = 5. Note that since Ap is not Hurwitz, in view of the saturation constraint, only local asymptotic stability for the closed-loop system with discrete-time (or continuous-time) control can be guaranteed. We consider a stabilizing continuous-time dynamic controller [START_REF] Cloosterman | Controller synthesis for networked control systems[END_REF], whose matrices are given by ( [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF], p.312): Āc = 0, Bc = 10, Cc = 1, Dc = 2, Ēc = 0.

Next, we illustrate the application of our results under different perspectives of analysis.

A. Influence of δ

Consider a discretized version of the controller obtained from [START_REF] Cloosterman | Stability of networked control systems with uncertain timevarying delays[END_REF]. Let us first consider the influence of the discretization period δ. With this aim, we apply the conditions of Theorem 1 considering a periodic sampling, i.e. T k = δ = T1 = T2, ∀k. It follows that no solution for these conditions can be found for δ larger than 0.08s.

B. Robustness with respect to aperiodic sampling

Consider the discretized version of the controller obtained from [START_REF] Cloosterman | Stability of networked control systems with uncertain timevarying delays[END_REF] with δ = 0.02s, which is given by: Ac = 1; Bc = 0.2; Cc = 1; Dc = 2.

(32)

We evaluate now what happens under aperiodic sampling. For this we consider an allowable intersampling interval T k ∈ [T1, T2] centered at δ, i.e. T1 = δµ and T2 = δ + µ. Considering Problem P1, the influence of µ on the size of the obtained RAS is shown in Figure 1. Note that, as expected, the larger is µ (i.e. the larger is the admissible interval [T1, T2]), the smaller tends to be the region of stability E (P, α). Moreover, the larger region is obtained with µ = 0, which corresponds to the periodic sampling case. Regarding Problem P2, conditions ( 16), ( 17), ( 18) become unfeasible when µ becomes larger than ±24%δ.

C. Robust periodic implementation

Consider now the dynamic discrete-time controller with the matrices given in (32), but that the actual periodic sampling period is different from δ = 0.02s. In other words, we consider T1 = T2 = T = 0.02. It can be observed that conditions ( 16), ( 17) and ( 18) remain feasible when reducing the sampling T much below 0.02, although the size of the stability domain is significantly degraded. Actually this domain shrinks as T decreases from 0.02. On the other hand, when T is increased above the value of 0.02, the problem remains feasible until T = 0.29. This is illustrated in Figure 2 where it can be observed the influence of the sampling period T on the size of the stability domain. For comparison purposes, the RAS obtained (from the application of Proposition 3.1 in [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF]) considering the continuous-time implementation of the controller is also depicted in the figure .  An interesting aspect of this evaluation is that, with T = T1 = T2 larger than δ, initial conditions not belonging to the region of attraction corresponding to the continuous-time closed-loop system (i.e. obtained from the connection between (1) and [START_REF] Cloosterman | Controller synthesis for networked control systems[END_REF] with ū(t) = u(t)) belong to the region of attraction obtained with the discretetime implementation of the controller. This is illustrated in Figure 3, t i m e ( s ) Fig. 3. Time response of the closed-loop system. Left-hand side: with discretized controller given by matrices in (32) and implemented with a periodic sampling T = 0.1. Right-hand side: with the continuous-time controller given by matrices in [START_REF] Zhang | Stability of Networked Control Systems[END_REF].

which compares the time evolution of the closed-loop system with the controller (32) implemented with T = 0.1 and that one with the PI continuous-time controller [START_REF] Zhang | Stability of Networked Control Systems[END_REF], for the initial condition xp(0) = [ -3.3 0 ] and xc(0) = 0. In fact, this initial condition belongs to the RAS obtained for the discrete-time controller, but does not belong to the one associated to the continuous-time implementation.

VII. CONCLUSION

This paper provided theoretical conditions for the stability analysis of a closed-loop system composed of a linear continuous-time plant connected to a sampled-data dynamic output feedback controller through a saturating input. A hybrid model for the closed-loop system has been proposed. Then, considering a looped functional approach, conditions for assessing the regional asymptotic stability of the origin for the closed-loop system under aperiodic sampling have been derived in the form of matrix inequalities. From these conditions, LMI-based optimization problems have been proposed to compute estimates of the region of attraction of the origin for the hybrid nonlinear closed-loop system, as well as to maximize the interval of admissible inter sampling time for which the regional asymptotic stability of the closed-loop system can be ensured under aperiodic sampling.

The extension of the approach for the synthesis of stabilizing controllers taking into account the possibility of aperiodic sampling is an open issue. In this case, from the hybrid nature of the system, it is not possible to apply transformations like the ones proposed for instance in [START_REF] Scherer | Multiobjective output-feedback control via LMI optimization[END_REF] and other works to convexify the problem. Hence, to find convex (or quasi-convex) conditions to address this problem can be seen as a challenge and is the subject of ongoing work.

Fig. 1 .

 1 Fig. 1. Influence of the interval [T 1 , T 2 ] with T 1 = δµ and T 2 = δ + µ on the ellipsoidal regions of stability in the system plane (xc = 0) considering δ = 0.02 and the following cases: µ = 0 (dashed-dotted line), µ = 10%δ (dotted line), µ = 20%δ (solid line) and µ = 24%δ (dashed green line).

Fig. 2 .

 2 Fig. 2. Ellipsoidal regions of stability in the system plane (xc = 0). Influence of a periodic sampling T (T 1 = T 2 = T ) when the controller is discretized with δ = 0.02. Cases with T = δ = 0.02 (dashed-dotted line), T = 0.01 (dotted line), T = 0.1 (solid line) and T = 0.29 (solid line). The continuoustime case is plotted in dashed green line.

  c r e t e -t i m e c o n t r o l l e r , T k = 0 . 1 t i n u o u s -t i m e c o n t r o l l e r