Antonio Luca

Federico Ricci
email: lricci@imf.org

Trionfetti

Ricci

Federico Trionfetti
email: federico.trionfetti@univmed.fr

Fei Liu

Maria Victoria Fazio

Productivity, networks, and export performance: evidence from a cross-country firm dataset

Keywords: F11, F12, F14 new-new trade theory, export probability

This paper uses a newly assembled multi-country multi-industry …rm-level dataset to test the e¤ect of productivity and networking on the export probability of …rms. Results are in line with the new-new trade theory and with the literature on the information value of networks. Firms are more likely to export if they are more productive, larger, and if they bene…t from foreign networks (ownership and …nancial linkages), domestic networks (chamber of commerce, links to regulation), and communication networks (E-mail, internet). Firms bear a lower probability of exporting if they have state or labor networks. Overall, …rms with better network connections by one standard deviation enjoy a 15% higher probability of exporting.

Introduction

Even a cursory glance at the data reveals that not all …rms are engaged in international trade. This diversi…ed reality about the export performance of …rms has recently found a theoretical collocation in general equilibrium models such are those proposed by [START_REF] Bernard | Plants and Productivity in International Trade[END_REF] and [START_REF] Melitz | The impact of trade on intra-industry reallocations and aggregate industry productivity[END_REF], which rely mainly on …rms'heterogeneity in productivity. The positive e¤ect of productivity on the probability of exporting is well documented in many studies. Evidence in this sense is found inter alia in Aitken, [START_REF] Brian | Spillovers, foreign investment, and export behavior[END_REF], [START_REF] Roberts | The Decision to Export in Colombia: An Empirical Model of Entry with Sunk Costs[END_REF], Bernard andJensen (1999, 2004). Our results con…rm that the probability of exporting is increasing with productivity using …rm level data from many countries and sectors.

A di¤erent strand of literature has studied the role of networks in facilitating international trade; see, e.g., Rodrick (2000), [START_REF] Rauch | Business and social networks in international trade[END_REF], [START_REF] Casella | Anonymous market and group ties in international trade[END_REF], [START_REF] Rauch | Ethnic Chinese networks in international trade[END_REF], [START_REF] Wagner | Immigration and the trade of provinces[END_REF], [START_REF] Rauch | Overcoming Informational Barrieres to International resource Allocation: Prices and Ties[END_REF], [START_REF] Combes | The trade creating e¤ects of business and social networks: evidence from France[END_REF]. Networks and connections reduce the disadvantage presented by informational barriers, and hence reduce the cost of economic transactions. In line with this literature, we …nd that …rms are more likely to export if they bene…t from foreign networks (ownership and …nancial linkages), domestic networks (chamber of commerce, links to regulation), communication networks (E-mail, internet). Firms bear instead a lower probability of exporting if they are a¤ected by state or labor networks.

Our results suggest a possible link between the two literatures cited above. In the …rst literature, productivity and the …xed exporting cost are the key factors generating the partition between exporters and non exporters. For a given …xed exporting cost, the higher the …rm productivity the more likely is that the …rm exports; for a given level of …rm productivity, the lower the …xed exporting cost faced by the …rm the more likely it is that the …rm exports. In the second literature, the barrier to international trade is typically informational (the cost of acquiring knowledge about foreign suppliers, foreign consumers, foreign business and regulation practices, and the cul-tural environment in general). These informational barriers may be seen as a form of …xed exporting cost. Thus, one may postulate that foreign network connection reduce …xed exporting costs thereby increasing the probability of exporting. We are agnostic as to whether this is indeed the causal chain but the empirical evidence gives plausibility to this hypothesis opens the way to interesting research avenues.

The remainder of the paper is as follows. Section 2 derives the estimated relationships, Section 3 presents the data, Section 4 shows the empirical results, and Section 5 concludes.

Estimated relationship

In this section we draw two empirically veri…able conclusions from the reading of the theoretical literature.

In theoretical models in the vein of [START_REF] Melitz | The impact of trade on intra-industry reallocations and aggregate industry productivity[END_REF], a …rm exports if foreign pro…ts are non-negative and does not export otherwise. Foreign pro…ts are given by foreign revenues minus variable cost (including variable trade cost), minus …xed exporting cost. Foreign revenues increase with productivity and variable cost decrease with productivity. Therefore, for given …xed exporting cost:

1. the probability of exporting is increasing with …rm productivity.

In the theoretical literature on networks a reduction in informational barriers increases international trade (possibly because they have lower …xed exporting cost). This possible link between …xed exporting costs and informational barriers is suggestive but not crucial for our purposes. From the literature on networks we simply draw that 2. a …rm with better networking connections has a higher probability of exporting.

These two conclusions may be summarized in the following relationships to be tested empirically:

Probability of Exporting = f(productivity, networking connections)

(1)
The empirical methodology adopts the following probit regression of the export dummy cij for …rm j in industry i of country c.

cij = & ic + X cij + cij (2)
The dependent variable is regressed on various determinants X cij : productivity, employment, and variables related to …rms' connections to networks. Regressions include country-industry …xed e¤ects to control for omitted variables which are speci…c to each sector in each country, and are estimated with robust standard errors.

Data description

We use the Enterprise Survey …rm-level dataset which is based on a World Bank survey of a representative sample of an economy's private sector. The surveys cover a broad range of business environment topics including infrastructure, crime, competition, access to …nance, corruption, and performance measures. 1We derive an export dummy (=1 if the …rm exports, and 0 otherwise) from the information on whether the …rm exports. Tables 1 and2 present the percentage of …rms that export, by country and by industry respectively. Note that about 30 percent of …rms in our sample are exporters but the percentage is 8 percent for smaller …rms, is about 25 percent for medium size …rms, and is 51 percent for larger …rms.

[TABLE 1 AND 2 HERE] Productivity is proxied by the ratio of sales to employment, in log terms, which is commonly used when total factor productivity is not available (sales are converted in US.$ in order to make them comparable across countries). 2We also add employment as a proxy for …rms'size.

There is no standard way to measure all possible network connections which in ‡uence …rms' economic activity. Within the economic indicators available in the dataset we use the following variables: "Firms is in a joint venture", "Firms uses E-mail", "Firm has website", and "Firm belongs to a chamber of commerce", which are dummy variables: 1 if the …rm does, 0 otherwise. "Foreign …nancing", which is measured by the percentage of working capital …nanced by foreign-owned banks. "Foreign participation", which is measured by the percentage of the …rm owned by foreign private sector. "Time spent with gov't regulations", which is measured by the percentage of senior management's time spend dealing with government regulations. "State participation", which is measured by the percentage of the …rm owned by the state. "Unionization" which is measured by the percentage of labor force unionized.

Most of the network variables employed, such as joint venture, chamber of commerce, foreign …nancing and participations directly relate to networks. Some other variables relate to network more indirectly and would also relate to other important economic factors such as technology, institution, culture, and so on. For instance, "E-mail" and "website" may be interpreted as technology variables that reduce the …xed cost of gathering information; but …xed costs are reduced precisely because this technology allows …rms to be part of a network (the web). Another example is "Time spent with gov't regulations" which may be a function of the country institutional and regulatory environment; however, once controlling for country-industry …xed e¤ects, the residual e¤ect at the …rm level is likely to re ‡ect the role …rms' network connections useful to penetrate heavily regulated markets or mar-kets upon which governments have a large in ‡uence. Another example is "Unionization" which is a network of workers re ‡ecting labor market institutions and regulations, and often presenting a negative externality for the …rm. What all these variable have in common is that they can all be related to networking.

In order to gauge the overall e¤ect of network connections (beyond that of any single variable) we construct a synthetic index of all variables that can directly or indirectly characterize a network connection. We name it "network index". Such network index is the sum of all network variables after they are standardized (i.e. demeaned, and divided by their standard deviation), and the sign of each component is adjusted to re ‡ect a positive e¤ect on the probability of exporting. Only the sign of two variables is changed before summing them (the percentage of workforce unionized and state participation), as these two variables tend to appear with a negative sign, and presumably characterize costs or frictions. For convenience, the network variable is normalized (zero mean and standard deviation equal to 1).

For each underlying variable, we carefully check the data for extreme values (dropping few …rms whose employment or sales would be reported as larger than national aggregate indicators at the country level available in standard macroeconomic databases such as World Bank or IFS), consistency checks (dropping …rms where reported categories of employment would not add up to total reported …rm employment, or reported categories of capital would not add up to total reported …rm capital), and for outliers (observations deviating from the country-industry mean by more than three standard deviations, in log terms). After data cleanup, our data set encompasses about 8,000 …rms in 24 industries and 32 countries. 3 Summary statistics for the variables used are in Table 3.

[TABLE 3 HERE] 4 Empirical results.

The e¤ect of productivity, size and network connections on the probability of exporting is explored in Tables 4, 5 and 6. Table 4 uses OLS to check the contribution of each set of variables to the variability of the export dummy. Country e¤ects have a stronger explanatory power than sector e¤ects (with an R 2 of about 13 percent and 7 percent respectively). Both country and industry e¤ects have an explanatory power close to the sum of the two e¤ects separately. With country-industry e¤ects, the R 2 rises to about 27 percent. Productivity variables and network variables, when entered alone, have a similar explanatory power (17 percent versus 15 percent respectively). When all variables enter the regression, the R 2 is about 36 percent. Replacing all network variables with the synthetic network index implies only a minimal reduction in explanatory power (the R 2 is about 35 percent).

[TABLE 4 HERE.] Tables 5 and6 show probit regressions of the export dummy on productivity, size, and networks variables. 4 In Table 5 robust standard errors are employed, while in Table 6 we cluster errors by country and industry and results are unchanged. The Tables show that both …rms'productivity and …rms'size have a positive and signi…cant e¤ect which con…rms results well established in the literature. The tables also show importance of the explanatory variables related to various forms of …rms'connections to networks (the following e¤ects are based on the …rst column of these tables, when all variables are si- 4 The coe¢ cients in the tables with probit regressions represent marginal e¤ect, i.e. the change in the probability for an in…nitesimal change in each independent, continuous variable. multaneously included). Firms'foreign connections are all positively related to the probability of exporting: an increase in Foreign …nancing or Foreign participation by 10% is associated with a 3% and 2% (respectively) increase in the probability of exporting, while …rms with a Foreign joint venture are 16% more likely to export. Firms are also more likely to export if they have more advanced communication technologies (…rms employing E-mail or websites are associated with a 14% or 7% higher probability of exporting, respectively), or stronger domestic connections (…rms which are member of a chamber of commerce are associated with a 13% higher probability of exporting) or by the time spent with government regulations (a 10% larger share of senior management's time spent dealing with government regulations is associated with a 2% higher probability of exporting). Younger …rms tend to be associated with a higher export probability (every additional 10 years of existence are associated with a 2% decrease in the probability of exporting). 5 Lastly, …rms a¤ected by state or labor networks tend to export less: a 10% increase in state participation is associated with a 3% lower probability of exporting; the e¤ect of unionization is very small when controlling for all variables (in column 1 of both table 5 and6), but in column 12 of both tables a 10% increase in state participation is associated with a 1% lower probability of exporting). 6 Overall, an improvement the network index by one standard deviation is associated with a higher probability of exporting of about 15%.

[TABLE 5 and 6 HERE] 5 As …xed e¤ects control for the (slight) di¤erence in years of the survey across countries, the variable related to "Year …rm began operations", is automatically equivalent to the opposite of the age of the …rm. Note that this result is in contrast with the literature on advanced economies (see for example Barba Navaretti et al., 2010), where older …rms are more likely to export, and may re ‡ect the more vibrant nature of the countries in our sample. 6 The variables "Firms is in a joint venture", "Firms uses E-mail", "Firm has website", and "Firm belongs to a chamber of commerce", are dummy variables: 1 if the …rm does, 0 otherwise. "Foreign …nancing" is measured by the percentage of working capital …nanced by foreign-owned banks. "Foreign participation" is measured by the percentage of the …rm owned by foreign private sector. "Time spent with gov't regulations" is measured by the percentage of senior management's time spend dealing with government regulations. "State participation" is measured by measured by the percentage of the …rm owned by the state. "Unionization" is measured by the percentage of labor force unionized.

Conclusions.

What determines the export performance of …rms? Recent theoretical developments in the trade literature point at productivity as a key determinants, while a separate strand of literature has highlighted over the past decade the importance of networking connections. The empirical literature has extensively documented the importance of productivity for the export performance of …rms, while evidence of the role of networks is very limited. This paper makes use of a cross-country and multi-sector …rm level data to jointly con-…rm these theories. We …nd that …rms are more likely to export if they are more productive and larger, but also if they bene…t from foreign networks (…nancial, ownership, or joint-venture links), domestic networks (chamber of commerce, links to regulation), and communication networks (E-mail, Website), while their export performance su¤ers from state and labor networks (government participation or unionization).

dia 2000 ;

 2000 Kosovo 2003; Kyrgyzstan 2003; Lithuania 2004; Madagascar 2005; Malawi 2005; Mauritius 2005; Morocco 2004; Nicaragua 2003; Pakistan 2002; Philippines 2003; Poland 2003; Serbia 2003; South Africa 2003; Sri Lanka 2004; Tajikistan 2003; Tanzania 2003; Thailand 2004; Turkey 2005; Uzbekistan 2003; Vietnam 2005; Zambia 2002.

SIZE

TABLE 4 OLS regression of the export dummy on fundamentals controlling for country*industry effects

 4

		1	2	3	4	5	6	7	8	9	1 0	1 1	1 2	1 3
	VARIABLES
	Log of productivity	0.06	0.051	0.055	0.046	0.05	0.04	0.046	0.052	0.051	0.054	0.049	0.052	0.064
		(0.000)*** (0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***
	Log of employment	0.119	0.128	0.137	0.123	0.135	0.11	0.121	0.111	0.129	0.134	0.147	0.142	0.135
		(0.000)*** (0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***
	Foreign Financing	0.003		0.004
		(0.023)**		(0.000)***
	Foreign Joint Venture	0.159			0.13
		(0.001)***			(0.000)***
	Foreign Participation	0.002				0.002
		(0.000)***				(0.000)***
	Firm uses Email	0.14					0.201
		(0.000)***					(0.000)***
	Firm has Website	0.067						0.131
		(0.016)**						(0.000)***
	Firm part of Chamber of Commerce	0.128							0.093
		(0.000)***							(0.000)***
	Time spent with gov't regulations	0.002								0.001
		(0.014)**								(0.021)**
	State Participation	-0.003									-0.002
		(0.029)**									(0.000)***
	Percentage of workforce unionized	0										-0.001
		-0.359										(0.000)***
	Year firm began operations	0.002											0.003
		(0.001)***											(0.000)***
	Network index													0.155
														(0.000)***
	Observations	2311	7526	6591	7422	7422	7183	7067	5733	7198	5859	5643	7064	2311
	Pseudo R2:	0.314	0.315	0.308	0.323	0.316	0.335	0.323	0.339	0.322	0.315	0.291	0.312	0.308
	Robust p-values in parentheses
	*** p<0.01, ** p<0.05, * p<0.1

TABLE 5 Probit regression of the export dummy on fundamentals controlling for country*industry effects 15

 5 )*** (0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)*** )*** (0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***(0.000)***

		1	2	3	4	5	6	7	8	9	1 0	1 1	1 2	1 3
	VARIABLES
	Log of productivity	0.06	0.051	0.055	0.046	0.05	0.04	0.046	0.052	0.051	0.054	0.049	0.052	0.064
	(0.000Log of employment 0.119	0.128	0.137	0.123	0.135	0.11	0.121	0.111	0.129	0.134	0.147	0.142	0.135
	(0.000Foreign financing 0.003		0.004
		(0.038)**		(0.000)***
	Foreign joint venture	0.159			0.13
		(0.001)***			(0.000)***
	Foreign participation	0.002				0.002
		(0.000)***				(0.000)***
	Firm uses Email	0.14					0.201
		(0.004)***					(0.000)***
	Firm has website	0.067						0.131
		(0.019)**						(0.000)***
	Firm part of chamber of commerce	0.128							0.093
		(0.000)***							(0.000)***
	Time spent with gov't regulations	0.002								0.001
		(0.033)**								(0.028)**
	State participation	-0.003									-0.002
		(0.012)**									(0.002)***
	Unionization	0										-0.001
		-0.394										(0.000)***
	Year firm began operations	0.002											0.003
		(0.001)***											(0.000)***
	Network index													0.155
														(0.000)***
	Observations	2311	7526	6591	7422	7422	7183	7067	5733	7198	5859	5643	7064	2311
	Pseudo R2:	0.314	0.315	0.308	0.323	0.316	0.335	0.323	0.339	0.322	0.315	0.291	0.312	0.308
	Robust p-values in parentheses
	*** p<0.01, ** p<0.05, * p<0.1

TABLE 6 Probit regression of the export dummy on fundamentals controlling for country*industry effects: errors clustered at country*industry level

 6

The survey is conducted by private contractors and is answered by business owners and top managers. "The manufacturing and services sectors are the primary business sectors of interest. This corresponds to …rms classi…ed with ISIC codes 15-37, 45, 50-52, 55, 60-64, and 72 (ISIC Rev.3.1). . . . Formal (registered) companies with 5 or more employees are targeted for interview. Services …rms include construction, retail, wholesale, hotels, restaurants, transport, storage, communications, and IT. Firms with 100% government/state ownership are not eligible to participate in an Enterprise Survey. . . . In each country, businesses in the cities/regions of major economic activity are interviewed." (http://www.enterprisesurveys.org/Methodology/).

See, for example, Helpman et al. (2004), Eaton et al. (2008), and the survey o¤ered by Wagner (2007).

The industries in our sample are: Textiles; Leather; Garments; Agroindustry; Food; Beverages; Metals and machinery; Electronics; Chemicals and pharmaceutical; Construction; Wood and furniture; Non-metallic and plastic materials; Paper; Sport goods; IT services; Other manufacturing; Accounting and …nance; Advertising and marketing; Other services; Retail and wholesale trade; Transport; Mining and quarrying; Auto and auto components; Other transport equipment. The countries in the sample (with the respective year of survey) are: Algeria 2002; Bangladesh 2002; Chile 2004; China 2003; Ecuador 2003; Egypt 2004; El Salvador 2003; Ethiopia 2002; Guatemala 2003; Honduras 2003; In-

(1) Dummy va ri a bl e: 1 if the fi rm does , 0 Otherwi s e.

(2) Percenta ge fi nance for working ca pi ta l : forei gn-owned banks ;

(3) Percenta ge of fi rm owned by forei gn private sector;

(4) Percenta ge of senior management's ti me spent dealing with government regul a ti ons ;

(5) Percenta ge of fi rm owned by the state (6) Percenta ge of labor force unionized.

(7) See text for definition.

Table 3 Summary Statistics


Ibtissem Boulanouar

Abderrezak Rachedi

Stephane Lohier

Gilles Roussel

Energy-Aware Object Tracking Algorithm using Heterogeneous Wireless Sensor Networks

Keywords: Wireless Sensor Networks, Wireless Multimedia Sensor Network, tracking, mobile target. I

Mobile target tracking consists in detecting and locating targets in a surveillance region. Wireless Multimedia Sensor Networks (WMSN) can provide more information about the moving target than classical Wireless Sensor Networks (WSN), but their process and the transmission of multimedia data are costly which reduces the network's lifetime. In this paper, we propose a low-cost new solution for tracking a mobile target called EAOT (Energy Aware Object Tracking). It consists of a distributed cooperative algorithm that runs in heterogeneous Wireless Sensor Networks composed of both scalar and multimedia sensors. The scalar sensors (MS) are equipped with a motion detector; their role is to detect the target and then activate the camera sensors (CS) through message exchanges. We conducted simulations to evaluate the performance of the proposed mechanism in terms of tracking accuracy, number of exchanged messages, energy consumption, and spatial average deviation. The results are compared to two existing solutions: one based on simple camera sensors which always run in active mode and the other one based on the optimal camera node selection (OCNS) to track the target. The obtained results show that for the same number of nodes, EAOT consumes between 2.6% and 44% less energy than OCNS depending on the distribution of MS/CS.

INTRODUCTION

Wireless Sensor Networks (WSN) are small devices networks with resource constraints which can communicate with each other on short distances to harvest and process data.

To make them multifunctional, these sensors have been equipped with small cameras and microphones in order to capture and retrieve multimedia contents. These Wireless Multimedia Sensor Networks (WMSN) allow a wide range of applications such as: intrusion detection, target tracking, habitat and health care monitoring. In addition to the common characteristics shared with WSN, the WMSN have special features: high bandwidth demand, specific QoS requirements, sector sensing range, etc. Moreover, wireless Camera Sensors (CS) unlike classical video systems can be used in both indoor and outdoor environments, where energy and network infrastructure are not available and where no human intervention is possible. They offer a wider panel of applications whether for environmental, industrial or military monitoring [START_REF] Akyldiz | A survey on Wireless Multimedia Sensor Networks[END_REF]. Most recent studies in WMSN focus on increasing the network lifetime by using an energy aware routing protocol, providing Quality of Service (QoS) and security in communications and by placing nodes in order to increase the coverage area.

In this paper, we focus on object tracking in WMSN and more particularly on non-communicating object tracking like human beings, animals, etc. Mobile target tracking consists in retrieving the target's successive coordinates during its moving through the monitoring area. This topic has been well studied in classical WSN, particularly in the case of communicating object tracking. With the emergence of small CMOS cameras, tracking mobile object using video sensor wireless networks has received more attention. The main difficulty is to handle the tradeoff between ensuring the tracking accuracy and maximizing the network lifetime. Indeed, the solution for the best tracking accuracy is to keep all camera sensors active (naive technique), but this is impossible because of energy consumption due to processing and transmission of the multimedia data retrieved. The best solution is to only activate the camera sensors able to localize the target.

We distinguish three main classes of tracking techniques [START_REF] Jin | Dynamic clustering for object tracking in wireless sensor networks[END_REF]: naive [START_REF] Heinzelman | Energy-Efficient Communication Protocol for Wireless Microsensor Networks[END_REF], predictive based [START_REF] Pahalawatta | Optimal sensor selection for video-based target tracking in a wireless sensor network[END_REF] [5] and dynamic clustering [START_REF] Zaho | Collaborative Signal and Information Processing:An Information Directed Approach[END_REF]. The naive technique means that all nodes are always in active mode, what maximizes the energy consumption. The second one consists in using a predictive model to anticipate the future position of the target. This technique is not appropriate for some applications like intrusion detection. Indeed, the predictive model used cannot ensure 100% accuracy of the calculated trajectory. Finally, in the dynamic clustering, a cluster of nodes is selected at each step of the target progression in the area of interest. While the selected cluster of nodes performs the tracking, the other nodes remain in sleeping mode. So this technique saves the network's energy. In order to find a better tradeoff between accuracy and energy saving, we introduce a Heterogeneous Wireless Sensor Network (HWSN) composed of both scalar and camera sensors. The scalar sensor consumes less energy than the multimedia sensor which allows us to keep them always in active mode. We propose an Energy Aware Object Tracking (EAOT) protocol based on message exchange, where the scalar sensors handle the object detection and only activate the camera sensors that can localize the target object. The main characteristics of the proposed EAOT are: 1) immediate detection of the target when it enters the monitoring area 2) performing the tracking in a low-cost manner 3) object localization based on three possibilities: the real coordinates, a visual observation of the target and as worst case approximate coordinates.

The rest of this paper is organized as follows: in section II, we review the related work in WMSN research topics. Section III describes the proposed protocol called EAOT. We present and discuss the simulation results obtained from our solution EAOT and other existing solutions in section IV. Finally, section V concludes the paper and presents our future works.

II. RELATED WORK

The mobile object tracking is one of the Wireless Multimedia Sensor Networks (WMSNs) applications. This application consists in locating the mobile target at every step of its progression in the surveillance area. The tracking solutions can be classified in three main categories [START_REF] Jin | Dynamic clustering for object tracking in wireless sensor networks[END_REF]: naive, predictive based and dynamic clustering. Because of the high cost of the naive technique, it is unreasonable to use such a method in WMSNs where energy is a precious resource.

In the predictive based technique, a predictive model is used to predict the future position of the mobile target. An adapted Kalman Filter is used in [START_REF] Pahalawatta | Optimal sensor selection for video-based target tracking in a wireless sensor network[END_REF] to calculate the future sensor utility depending on past data collected. In [START_REF] Zaidi | Mobility Tracking based on Autoregressive Models[END_REF] the authors retrieve the mobility parameters of the target and use an autoregressive model to integrate them and predict the future trajectory.

Dynamic clustering [START_REF] Zaho | Collaborative Signal and Information Processing:An Information Directed Approach[END_REF] is the most used technique in literature. A cluster of nodes is selected at every step of the evolution of the mobile target in the surveillance region. A cluster head is selected depending on predefined criteria. The advantage of dynamic clustering is that the sensors which are not selected to be cluster members remain in sleeping mode. Thus, this technique allows saving energy.

Papers described here correspond to the third technique .The authors in [START_REF] Liu | Optimal Node Selection for Target Localisation in Wireless Camera Sensor Netwoks[END_REF] propose a distributed solution based on node collaboration to select the optimal subset of camera sensors that participate in the target location process. This solution consists of a two phase algorithm that maintains the adequate number of nodes in active mode. When the target is detected by a camera sensor, it broadcasts its location and orientation to the neighbor. Each of them calculates the probability of detecting the target. If the detecting probability reaches a predefined threshold, the concerned node activates its camera to take part in the locating process. This collaborative algorithm is also merging with a wake up algorithm that maintains the adequate number of nodes in active mode. This solution presents a good tradeoff between energy consumption and tracking accuracy but the target can enter the surveillance region without being immediately detected because of the wake up algorithm. SensEye [START_REF] Kulkarni | SensEye: A Multitier Camera Sensor Network[END_REF] is the first solution that introduces the concept of heterogeneous network. The authors propose three-tier camera sensors; every tier supports a specific task. The first tier assumes target detection and localization while the second one performs target recognition and the last tier assumes target tracking.

SensEye can only be used in indoor environments where energy is available which considerably limits the types of applications where we can use it.

The authors of [START_REF] Newell | Distributed Collaborative camera actuation for redundant data elimination in Wireless Multimedia Sensors[END_REF] also used the concept of heterogeneous networks but for different objectives. Indeed, the activation goal is the event coverage, they used the scalar sensors to determine the event boundary and actuate the necessary camera sensors. The main objective is the elimination of data redundancy.

III. ENERGY AWARE OBJECT TRACKING

In this section, we present the proposed protocol called Energy Aware Object Tracking (EAOT). Its main objectives are as follows: 1) having a visual observation of the target, which permits its identification 2) having target coordinates to use them for specific applications like intrusion detection 3) saving energy to increase the network's lifetime.

A. Preliminaries

In this subsection we give some definitions and assumptions. The heterogeneous concept of sensor nodes is described with the characteristics and role of each type of node. 2) Definition 2: The Camera Sensor (CS) is a wireless multimedia sensor equipped with both motion detector and video camera. Each object located in CS's Field of View (FoV) can be visually detected. FoV is a CS's directional of view and it is assumed to be a cone with angle and radius V (as illustrated in figure 1.B).

A B

The main difference between MS and CS is not only limited to the services but also to the energy consumption. CS requires more energy to run its cameras compared to MS.

3) Definition 3: We introduce the concept of approximate coordinates (x, y). Equations 1 and 2 show how these coordinates can be calculated.

Where PMS represents the probability that the target is detected by MS and N is the total number of MSs that detect the target in the given area at the same time.

The target is detected by MS when it is in its FoD. In other words: the probability PMS that a target at distance RMS is detected is the probability that RMS ≤ D, where D is the circle radius of FoD.

In [START_REF] Stolkin | Probability of detection and optimal sensor placement for threshold based detection systems[END_REF], the authors provide a probabilistic analysis of a simple detection system. Using their concept, we can calculate PMS as follows:

Where and are respectively the mean and the standard deviation. is the distribution function of standard normal random variable. PMS decreases while the RMS increases.

We assume the random deployment of scalar and multimedia sensors in the monitoring region. All deployed motion sensors (MS) have the same field of detection radius (FoD). The same assumption is valid for all deployed camera sensors (CS) which have the same field of view (FoV) parameters. Moreover, we assume that the nodes are static, each MS knows its location and each CS knows its location and orientation using one of the existing methods described in literature [START_REF] He | Range-free localization schemes for large scale sensor networks[END_REF]. A single mobile target at a time is expected to cross the monitored region by taking a random path. Every sensor can communicate with another one as long as it is in its transmission range, independently of the type of sensor. Initially we supposed that all cameras are inactive and the motion detectors are always active.

B. Solution overview

The proposed solution EAOT is a cooperative protocol based on two types of packets called: DETECTION and LOCALIZATION. The DETECTION packet is sent only by the MS node to inform the neighboring nodes, particularly CS nodes, about the presence of any object in its FoD area. The LOCALIZATION packet is sent only by the CS nodes to the first MS that informs it once they receive a certain threshold M of DETECTION packets from the MS nodes. Both packet types contain the coordinates of sender nodes. In addition, the packet size of DETECTION and LOCALIZATION is similar and equal to 64Kb. Figure 2 illustrates an example of object tracking with the proposed solution EAOT. MS1 is the first node that detects the target (the target trail is illustrated by the dashed line). MS1 broadcasts a DETECTION message containing its own coordinates to the nodes which are in its transmission range.

The life span of DETECTION message T is fixed to 20 seconds, chosen according to the target speed. Once MS2 detects the target, it also broadcasts a DETECTION message. When CS1 receives the second message from MS2 and if threshold M is set to 2, then it actuates its camera to localize the target. In this case, the target is in its FoV. CS2 retrieves the real target's coordinates and sends them only to MS1 via LOCALIZATION message. In the following step, CS2 activates its camera after the reception of the DETECTION messages from MS2 and MS3 but the target is not in its FoV. In this case, MS2 is responsible for assigning approximate target coordinates using previous equations.

Successively, MS1 and MS2 send the results to the sink (we suppose that all nodes have a route toward the sink).

The added value of the protocol EAOT consists in giving the target coordinates closest to its real position and in the worst case the approximate coordinates can be generated. In addition, the visual tracking is exploited with minimal energy consumption. When the CS receives at least M DETECTION messages: M plays an important role in tracking accuracy, when M increases the tracking accuracy also increases. This CS activates its camera to capture images of the target and checks if the target is in its FoV using background subtraction [START_REF] Piccardi | Background substraction[END_REF]. If so, the concerned CS performs the localization via any existing solution [START_REF] Oztrarak | Lightweight object localization with single camera in Wireless Multimedia Sensor Networks[END_REF] and sends it to the first MS. We suppose that this calculated coordinates are exact and error-free. Finally, it disables the camera.

A. MS flowchart

IV. PERFORMANCE EVALUATION

In this section, we present the simulation settings, evaluation metrics and we discuss the obtained results.

A. Simulation setting

We have used the NS-2 simulator [14] in order to evaluate EAOT. In all scenarios, the nodes have been randomly placed in a rectangular two dimensional area. The parameter values of the simulation are summarized in Table 1.

Detection range (D) of MS 15m

Depth of view (V) of CS 20m

Angle of view of CS /3

Size of message 64kb

M (number of activation messages)

messages

Life span of "Detection " message T

20s

EAOT with different distributions of nodes is compared to two other solutions (Table 2). In both of these only CSs are used. The first one is the naive solution which consists in deploying CSs. These CSs are always in active mode: capturing/processing multimedia streaming. The second one is the OCNS (Optimal Camera Node Selection), the operating principle of OCNS [START_REF] Liu | Optimal Node Selection for Target Localisation in Wireless Camera Sensor Netwoks[END_REF] is explained in section II.

B. Simulation results

Three metrics are used to compare our mechanism to other solutions: tracking precision, number of messages exchanged and energy consumption. Another metric is applied only for EAOT distributions: spatial average deviation.

1) Tracking precision

We measure it by considering the number of points collected for a defined trajectory. One point is represented by a two dimensional coordinate. For Naive and OCNS the retrieved points are the only points where we have an image of the target; we call them real target coordinates. In our solution the retrieved points can be a real target coordinates (extracted from the captured image) or an approximate target coordinates obtained from MS (equation 1 and 2). Figure 4 shows the average number percentage of retrieved points for every scenario depending on the number of nodes. 100 % represents the ideal tracking precision: for a pedestrian, the best number of retrieved points is 1 every meter.

Table 3 shows the real target coordinates percentage depending on the MS/CS distribution. We note that the best tracking accuracy is obtained by the naive approach because all CS nodes are always running and the visual target detection has a positive impact on the precision. Among the three EAOT schemes, EAOT2 represents the best compromise between real target coordinates and approximate ones and has a better target precision than OCNS. EAOT3 also has more retrieved points than OCNS but more approximate target coordinates. The main reason is the distribution of MS/CS. In EAOT3, we have more MS so it results in more approximate coordinates. EAOT1 has fewer points because there are not enough MSs to activate the CS.

2) Spatial average deviation

This metric represents the average deviation of the obtained trajectory vs. the real one. It is obtained by calculating the average of the distances between approximate target coordinates and real ones.

We suppose that the coordinates calculated by the CSs are exact and error-free. So, they are not considered in average deviation results. Distance d is calculated as follows:

(Xc, Yc) are the calculated coordinates and (Xr, Yr) are the real ones. Average deviation cannot be considered as tracking precision metric because it is only applied to EAOT schemes. It is the only solution where target coordinates can be real or approximate. Figure 5 shows the average deviation depending on the number of nodes and on the distribution. The results are obtained by calculating the distance between every retrieved point of EAOT schemes and the corresponding one in a real trajectory.

We observe that EAOT1 has the best results compared to other EAOT distributions due to the greater number of CSs than in other schemes. The average deviation of EAOT3 increases with the number of nodes. The main reason is the increase of MS: the more MSs we have, the more approximate coordinates are obtained.

3) Number of exchanged messages

We suppose that the message size is the same in all solutions. This metric shows the number of cooperative messages exchanged. For EAOT distributions it consists of both DETECTION/LOCALIZATION messages exchanged during the simulation. In OCNS, in addition to the cooperative messages, we add the messages of the wake up algorithm which does not exist in EAOT. Indeed, we assume that the MSs are always in active mode. Moreover, the messages transmitted to the sink are neglected in this metric for all solutions.

Figure 6 shows the average number of messages depending on the number of nodes. We remark that the number of exchanged messages increases with the number of nodes whatever the scenario. Only EAOT3 has more exchanged messages than OCNS, because of the high number of MSs which broadcast "DETECTION" messages when they detect the target.

4) Energy consumption

The percentage of energy consumption includes the energy consumed by cooperative message exchanges and the activation of MSs and CSs during one simulation. Figure 7 shows the average energy consumption percentage of the entire network depending on the number of nodes. 100% represents the total network energy available before the simulation starts. It is obvious that Naive consumes much more energy than the other solutions. Due to the use of scalar sensors for the detection, whatever the distribution of CS/MS EAOT always consumes less energy than OCNS. These results confirm the efficiency of EAOT for energy aware object tracking.

From the presented experiments and results, we can deduce that:

 Naive is an unrealizable solution. Despite the fact that it performs the tracking with the best precision, this solution is too costly and considerably reduces the network's lifetime.

 EAOT2 is the distribution which presents the best trade-off between tracking accuracy and energy consumption. It performs tracking with better precision, less energy and a smaller number of exchanged messages than OCNS.

 The simulation results can be improved by optimizing the node's positioning.

V. CONCLUSION AND FUTURE WORKS

Saving energy is a crucial issue in WMSNs whatever the research topic and the application field. Multimedia tracking of mobile targets being a costly application in terms of energy, the need for a low cost solution increases. In this paper, we proposed a new approach using Heterogeneous Wireless Sensor Networks (HWSNs). The main objective of our solution called EAOT for Energy Aware Object Tracking is to handle the tradeoff between tracking precision and energy consumption.

We conducted some simulations with different scenarios and distributions of MS/CS to evaluate our solution.

The results show that EAOT increases the tracking precision and reduces the energy consumption compared to OCNS where only CSs are used. The simulations performed also show that the distribution of MS/CS has an un-negligible impact on the obtained results. Our future works will consist in improving EAOT by adding the targets' identification and classification as well as the multi-target tracking.

Figure 1 .

 1 Figure 1. FoD and FoV 1) Definition 1: The Motion Sensor (MS) is a wireless scalar sensor equipped with motion detector. Each object located in the Field of Detection (FoD) of MS can be detected. FoD is represented by a circle with radius D (as illustrated in figure 1.A).

Figure 2 .

 2 Figure 2. An example of object tracking with EAOTEAOT is illustrated in figure3. We present two flowcharts: one for MS nodes and the other one for CS nodes. As explained before, the object tracking process starts when an MS detects the target. Each MS that has the target in its FoD broadcasts a DETECTION message. If it receives a LOCALIZATION message, it extracts the target coordinates and transmits them directly to the sink. Otherwise, after a predefined period of time T corresponding to the total time of transmitting/processing a DETECTION message, localization process, transmitting/processing a LOCALIZATION message and the travel time of the target between two nodes. The MS acts like a fusion center and handles the localization. It calculates the approximate coordinates. The obtained results are then transmitted to the sink by the MSs. The sink handles the reconstitution of the trajectory by piecing together the received coordinates.

Figure 3 (

 3 Figure3(B) illustrates the CS flowchart. When the CS receives at least M DETECTION messages: M plays an important role in tracking accuracy, when M increases the tracking accuracy also increases. This CS activates its camera to capture images of the target and checks if the target is in its FoV using background subtraction[START_REF] Piccardi | Background substraction[END_REF]. If so, the concerned CS performs the localization via any existing solution[START_REF] Oztrarak | Lightweight object localization with single camera in Wireless Multimedia Sensor Networks[END_REF] and sends it to the first MS. We suppose that this calculated coordinates are exact and error-free. Finally, it disables the camera.

B. CS flowchart

Figure 3 .

 3 Figure 3. EAOT flowchart

Figure 4 .

 4 Figure 4. Tracking precision/number of nodesTABLE 3. REAL TARGETCOORDINATES % Real target coordinates % EAOT1 74% EAOT2 66,66% EAOT3 41%

Figure 5 .

 5 Figure 5. Average deviation

Figure 6 .

 6 Figure 6. Number of exchanged messages according to the number of nodes

Figure 7 .

 7 Figure 7. Average energy consumption according to the number of nodes

TABLE 1 .

 1 SIMULATION PARAMETER

	Mac Layer protocol	IEEE 802.15.4
	Access Mode	CSMA/CA non-beaconed
	Number of nodes	20, 30, 40, 50
	Area size	100×100m²
	Target speed	1,38m/s (pedestrian)
	Simulation Time	150s
	Transmission range	3m

TABLE 2 .

 2

		SCENARIOS
	EAOT1	30% MS/70% CS
	EAOT2	50% MS/50% CS
	EAOT3	70% MS/30% CS
	naïve	100% CS, always in active mode
	OCNS	100% CS, solution proposed in
		[7]


Mathieu Claeys

Denis Arzelier

Didier Henrion

Jean-Bernard Lasserre

Measures and LMI for impulsive optimal control with applications to space rendezvous problems

This paper shows how to find lower bounds on, and sometimes solve globally, a large class of nonlinear optimal control problems with impulsive controls using semidefinite programming (SDP). This is done by relaxing an optimal control problem into a measure differential problem. The manipulation of the measures by their moments reduces the problem to a convergent series of standard linear matrix inequality (LMI) relaxations. After providing numerous academic examples, we apply the method to the impulsive rendezvous of two orbiting spacecrafts. As the method provides lower bounds on the global infimum, global optimality of the solutions can be guaranteed numerically by a posteriori simulations, and we can recover simultaneously the optimal impulse time and amplitudes by simple linear algebra.

Introduction

Optimal control problems are still an active area of current research despite the availability of powerful theoretical tools such as Pontryagin's maximum principle or the Hamilton-Jacobi-Bellman approach, that both provide conditions for optimality. However, numerical methods based on such optimality conditions rely on a certain number of assumptions that are often not met in practice. In addition, state constraints are particularly hard to handle for most of the methods.

On the other side, many numerical methods have been developed that deliver locally optimal solutions satisfying sufficient optimality conditions. However, the users of these methods are often left to wonder if a better solution exists. For example, in the particular case of impulsive controls, it is often not known if more regular solutions could provide a better cost. For a recent survey on impulsive control see e.g. [START_REF] Kurzhanski | Dynamic programming for impulse controls[END_REF] and the references therein. See also [START_REF] Gajardo | Minimal time sequential batch reactors with bounded and impulse controls for one or more species[END_REF] for a recent application and for more references 1 . For historical works see e.g. [START_REF] Neustadt | Optimization, a moment problem and nonlinear programming[END_REF][START_REF] Rischel | An extended Pontryagin principle for control systems whose control laws contain measures[END_REF][START_REF] Schmaedeke | Optimal control theory for nonlinear vector differential equations containing measures[END_REF] and also [START_REF] Bensoussan | Contrôle impulsionnel et inéquations variationnelles[END_REF]. This paper presents a method based on [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF][START_REF] Henrion | Nonlinear optimal control synthesis via occupation measures[END_REF] but covering a larger class of solutions, including impulsive controls. This algorithm provides a sequence of non-decreasing lower bounds on the global minimizers of affine-in-the-control polynomial optimal control problems. In particular, it may assert the global optimality of local solutions found by other methods. As importantly, the algorithm is also able to provide numerical certificates of infeasibility or unboundedness for ill-posed problems. Finally, in some cases, it is also possible to generate the globally optimal control law.

At the end of the paper, this method is successfully applied to the problem of coplanar space fuel-optimal linearised rendezvous. We show with two different examples from the literature that the proposed algorithm is able to retrieve the impulsive optimal solution conjectured by running a direct approach based on the solution of a Linear Programming (LP) problem. Without assuming the nature of the propulsion (continuous or impulsive), the obtained impulsive solution is certified to be a global fuel-optimal solution.

Contributions

The paper improves the method presented in [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF][START_REF] Henrion | Nonlinear optimal control synthesis via occupation measures[END_REF] in the following ways:

• Impulsive control can now be taken into account.

• Because controls are represented by measures and not by variables, the size of SDP blocks is significantly reduced. This allows to handle larger problems in terms of number of state variables as well as to reach higher LMI relaxations.

• Total variation constraints can be handled very easily.

These three improvements make it altogether possible to tackle problems such as consumption minimization for space rendezvous, the other significant contribution of this paper.

Notations

Integration of a function f : R n → R with respect to a measure µ on a set X ⊂ R n is written X f (x) dµ(x). The Lebesgue or uniform measure on X is denoted by λ whereas the Dirac measure concentrated at point x is denoted by δ x . A measure µ is a probability measure whenever dµ = 1. The support of measure µ is denoted by supp µ. The indicator function of set X (equal to one in X and zero outside) is denoted by I X .

F (X) is the space of Borel measurable functions on X, whereas BV (X) is the space of functions of bounded variation on X. R[z] is the ring of polynomials in the variable z. B(X) denotes the Borel σ-algebra associated with X.

If k ∈ N n denotes a vector of indices then x k with x ∈ R n is the multi-index notation for

x k i i . The degree of the index k is deg k = k i . Finally, N n d is the set of all indices for which deg k ≤ d, k ∈ N n .

The optimal control problem

This paper deals with the following nonlinear optimal control problem

V (x 0 ) = inf u(t)∈F ([0,T ]) m I(x 0 , u) = T 0 h(t, x(t))dt + T 0 H(t)u(t)dt + h T (x(T )) s.t. ẋ(t) = f (t, x(t)) + G(t)u(t), t ∈ [0, T ] x(0) = x 0 ∈ X 0 , x(T ) ∈ X T , x(t) ∈ X ⊂ R n .
(

) 1
where the dot denotes differentiation w.r.t. time and the prime denotes transposition. Criterion I(x 0 , u) is affine in the control u, and V is called the value function. It is assumed that all problem data are polynomials, meaning that all functions are in R[t, x], and that all sets are compact basic semialgebraic. Recall that such sets are those which may be written as {x : a i (x) ≥ 0, i = 1, . . . , m} for some family

{a i } m j=1 ⊂ R[x].
A mild technical condition (implying compactness of X) must be satisfied [13, Assumption 2.1], but it is often met in practice (for instance, an additional standard ball constraint

x 2 i ≤ r 2 enforces the condition). The reason for making these assumptions will be apparent in the later sections.

Without additional assumptions and constraints, the infimum in problem (1) is generally not attained in the space of measurable functions [START_REF] Young | Lectures on the calculus of variations and optimal control theory[END_REF]. For this reason, in this paper we consider problems for which controls are allowed to be generalized functions, i.e. measures, thereby extending the original formulation as follows:

V R (x 0 ) = inf w(t)∈BV ([0,T ]) m I(x 0 , w) = T 0 h(t, x(t))dt + T 0 H(t)dw(t) + h T (x(T )) s.t. dx(t) = f (t, x(t)) dt + G(t)dw(t), t ∈ [0, T ] x(0) = x 0 ∈ X 0 , x(T ) ∈ X T , x(t) ∈ X ⊂ R n (2
) where V R stands for the relaxed value function. In particular, in problem (2) controls may be impulsive: the (vector) control can been seen as a (vector) distribution of the first order and it is therefore the distributional derivative dw(t) of some (vector) function of bounded variation w(t) ∈ BV ([0, T ]) m , see e.g. [START_REF] Schmaedeke | Optimal control theory for nonlinear vector differential equations containing measures[END_REF] and [18, §4] or also [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF]Prop. 8.3].

The measure problem

In this section, we formulate problem (2) into a measure differential problem, a necessary step towards obtaining a tractable SDP problem. Optimal control problems involving measures have been introduced to accept solutions that are ruled out or ill-defined in classical optimal control, see e.g. [START_REF] Young | Lectures on the calculus of variations and optimal control theory[END_REF]. Multiple solutions, impulsive or chattering controls can be handled naturally by the associated measure problem. This section, rather than providing rigorous proofs, outlines the main ideas behind this transformation.

A few remarks are worth pointing out. First of all, it is crucial that G(t) be a matrix of smooth functions, an hypothesis automatically fulfilled by polynomials. As a matter of fact, multiplying distributions with such functions is a well-defined operation (unlike e.g. the product of two distributions). Therefore, except for some very particular cases [START_REF] Miller | The generalized solutions of ordinary differential equations in the impulse control problems[END_REF], G cannot be a function of states x j that could potentially present jump discontinuities. To simplify notations, we have simply assumed that G depends on t only2 . Secondly, in the absence of impulses, the distributional differential is the traditional differential, and the dynamics are classical differential equations with controls dw(t) = u(t)dt which are absolutely continuous with respect to the Lebesgue measure. Finally, it must be noted that state trajectories x(t) are themselves functions of bounded variations, being the sum of two such functions, and that this is their broadest class in the sense that there is no more general distribution solutions for the states [START_REF] Schmaedeke | Optimal control theory for nonlinear vector differential equations containing measures[END_REF].

Because distributional derivatives of functions of bounded variation on compact supports can be identified with measures [18, §50], the dynamics in problem (2) may be interpreted as a measure differential equation. As X ⊂ R n is assumed to be compact, by one of the Riesz representation theorems [10, §36.6], these measures can be put in duality correspondence with all continuous functions v(t, x(t)) supported on [0, T ] × X. We will use these test functions to define linear relations between the measures. Note that because continuous functions on compact sets can be uniformly approached by polynomials by virtue of the Stone-Weierstrass theorem, it is enough to consider polynomial test functions

v(t, x(t)) ∈ R[[0, T ] × X].
By Lebesgue's decomposition theorem [10, §33.3], we can split the control measures w(dt) into two parts: their absolutely continuous parts with density u : [0, T ] → R m (with respect to the Lebesgue measure) and their purely singular parts with jump amplitude vectors u t j ∈ R m supported at impulsive jump instants t j , j ∈ J, with J a subset of Lebesgue measure zero of [0, T ], not necessarily countable3 . We write

w(dt) = u(t)dt + j∈J G(t j )u t j δ t j (dt) to model jumps in state-space x + (t j ) = x -(t j ) + G(t j )u t j , ∀ j ∈ J.
Now, given an initial state x 0 ∈ X 0 and given a control w(t) ∈ BV ([0, T ]) m , denote by x(t) ∈ BV ([0, T ]) n the corresponding feasible trajectory. Then for smooth test functions

v : [0, T ] × R n → R, it holds T 0 dv(t, x(t)) = v(T, x(T )) -v(0, x(0)) = T 0 ∂v ∂t + ∂v ∂x ′ f dt + T 0 ∂v ∂x ′ Gudt + j∈J v(t j , x + (t j )) -v(t j , x -(t j )). (3)
We are going to express the above temporal integration (3) along the trajectory in terms of spatial integration with respect to appropriate and so-called occupation measures. For this purpose, define:

• The time-state occupation measure µ[x 0 , w(t)](A × B) = A I B (x(t)) dt, ∀A ∈ B([0, T ]), ∀B ∈ B(X)
which measures the occupation of A × B by the pair (t, x(t)) all along the trajectory. Note that we write µ[x 0 , w(t)] to emphasize the dependence of µ on initial condition x 0 and control w(t). However, for notational simplicity, we may use the notation µ. By a standard result on Borel measures on a cartesian product, the occupation measure µ can be disintegrated into

µ(A × B) = A ξ(B | t) dt, where ξ(dx | t) is the distribution of x ∈ R n , conditional on t ∈ [0, T ]. It is a stochastic kernel, i.e., -for every t ∈ [0, T ], ξ(• | t) is a probability distribution on X, and
-for every B ∈ B(X), ξ(B | •) is a Borel measurable function on [0, T ].
In our case, since the initial state x 0 and the control w(t) are given, the stochastic kernel ξ(dx | t) is well defined along continuous arcs of the trajectory as

ξ(B|t) = I B (x(t)) = δ x(t) (B), ∀t ∈ [0, T ] \ J, ∀B ∈ B(X). (4)
On the other hand, at every jump instant t j ∈ J, we let

ξ(B | t j ) = λ(B ∩ [x -(t j ), x + (t j )]) λ([x -(t j ), x + (t j )]) , ∀t j ∈ J, ∀B ∈ B(X).
This means that the state is uniformly distributed along the segment linking the state before and after the jump, the above denominator ensuring that ξ(• | t) has unit mass for all t.

• The control-state occupation measure

ν[x 0 , w(t)](A × B) = A ξ(B | t) dw(t), ∀A ∈ B([0, T ]), ∀B ∈ B(X).
• The final state occupation measure

µ T [x 0 , w(t)](B) = I B (x(T )), ∀B ∈ B(X T ).
With these definitions, Eq. ( 3) may be written in terms of measures as:

X T v(T, x) dµ T (x) -v(0, x 0 ) = (5) [0,T ]×X ∂v ∂t + ∂v ∂x ′ f dµ(t, x) + [0,T ]×X ∂v ∂x ′ Gdν(t, x) = [0,T ] X ∂v ∂t + ∂v ∂x ′ (f + Gu) ξ(dx | t) dt + j∈J v(t j , x + (t j )) -v(t j , x -(t j )).
Similarly, the criterion in (2) to evaluate the trajectory and the control reads:

I(µ, ν, µ T ) = [0,T ]×X h dµ + [0,T ]×X H dν + X T h T dµ T .
In view of the above formulation with occupation measures, one may now define a relaxed version (or weak formulation) of the initial (measure) control problem [START_REF] Bensoussan | Contrôle impulsionnel et inéquations variationnelles[END_REF]. First note that

V R (x 0 ) = inf w(t) I(µ[x 0 , w(t)], ν[x 0 , w(t)], µ[x 0 , w(t)])
where the infimum is taken over all the occupation measures defined above, corresponding to a given initial condition x 0 and control w(t). Second, instead of searching for a control w(t), we search for a triplet of measures that solves the infinite dimensional problem:

V M (x 0 ) = inf µ,ν,µ T I(µ[x 0 ], ν[x 0 ], µ T [x 0 ])
under the trajectory constraints ( 5) for all v ∈ R[t, x] and the support constraints supp

µ = supp ν = [0, T ] × X, supp µ T = X T .
The measures now depend only on initial condition x 0 , since they just have to satisfy linear constraints [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF]. This motivates the notation µ

[x 0 ], ν[x 0 ], µ T [x 0 ]
in the above problem. This problem is an obvious relaxation of problem (2) which is itself a relaxation of (1), hence

V M (x 0 ) ≤ V R (x 0 ) ≤ V (x 0 ).
In the remainder of the paper, we will deal with this relaxed version of the occupation measures problem. However, for a well-defined control problem (2) one expects that in fact V M (x 0 ) = V R (x 0 ) and that an optimal solution of the relaxed problem will be the triplet of occupation measures corresponding to an optimal trajectory of problem (2) with given initial state x 0 and control w(t). Note that for the standard polynomial optimal control problem (1), without impulsive controls, and under additional convexity assumptions, it has been proved in [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF] that indeed

V M (x 0 ) = V R (x 0 ) = V (x 0 ).

Initial state with a given distribution

Recall that the occupation measures defined in the previous section all depend on x 0 . Observe that if µ 0 is a given probability measure on X 0 ⊂ R n and if one now defines:

µ(A × B) = X 0 µ[x 0 ](A × B) dµ 0 (x 0 ), ν(A × B) = X 0 ν[x 0 ](A × B) dµ 0 (x 0 ), µ T (B) = X 0 µ T [x 0 ](B) dµ 0 (x 0 )
for all A ∈ B([0, T ]) and B ∈ B(X), then

I(µ[µ 0 ], ν[µ 0 ], µ T [µ 0 ]) = X 0 I(µ[x 0 ], ν[x 0 ], µ T [x 0 ])dµ 0 (x 0 )
becomes the expected average cost associated with the trajectories and with respect to the probability measure µ 0 on X 0 .

Therefore, the relaxed problem with measures now reads as follows:

V M (µ 0 ) = inf µ,ν,µ T I(µ[µ 0 ], ν[µ 0 ], µ T [µ 0 ]) = hdµ + Hdν + h T dµ T s.t. v dµ T -v dµ 0 = ∂v ∂t + ∂v ∂x ′ f dµ + ∂v ∂x ′ Gdν supp µ = supp ν = [0, T ] × X, supp µ T = X T . (6)
Note that in this case, the stochastic kernel ξ(dx|t) along continuous arcs of the trajectory is generally not a Dirac measure as in (4), unless µ 0 is a Dirac measure supported at x 0 and the optimal control w is unique.

By solving this relaxed problem we expect that its optimal value satisfies

V M (µ 0 ) = X 0 V M (x 0 )dµ 0 (x 0 ),
i.e. that V M (µ 0 ) is the expected average cost associated with optimal trajectories and with respect to the probability measure µ 0 on X 0 .

Free initial state

In this case, in addition to the control we also have the freedom of choosing the best possible initial state. For this purpose introduce an unknown probability measure µ 0 on X 0 . Then the relaxed problem with measures, analogue of ( 6), reads almost the same except that:

• we now optimize over µ, ν, µ 0 , µ T ;

• in the support constraints we introduce the additional constraint supp µ 0 = X 0 .

By solving this relaxed problem we now expect that its optimal value denoted V M (X 0 ) satisfies:

V M (X 0 ) = inf µ 0 V M (µ 0 ) = inf x 0 ∈X 0 V M (x 0 ).

Decomposition of control measures

All measures in (6) are positive measures, except for the signed measures ν which deserve special treatment for our purposes. Using the Jordan decomposition theorem [10, §34], these measures may be split into a positive part ν + and negative part ν -, that is ν = ν +ν -, both being positive measures.

This decomposition has the added benefit of providing an easy expression for the L 1 norm of the control, which is sometimes to be constrained or optimized in some problems. Indeed, define the total variation control measure by

|ν| = ν + + ν -.
The total variation norm of the measure ν is just the mass of |ν|, i.e., ν T V = d|ν|.

Handling discrete control sets

It is often desirable to restrict the set of admissible controls to be a subset of R. Here we will limit ourselves to the very important case of handling discrete control sets. Let us assume that controls u are only allowed to take their values in U = {u 1 , ..., u m }. Define ν i as the probability measures of choosing controls u i . Then clearly, the total probability of choosing one of the controls in U must be 1 at each time along the trajectory. Then the control measures ν are simply the linear combination of the probability measures weighted by their respective control values:

ν = i u i ν i .
Using the same method as in §3, we have

∀ v(t) ∈ R[t]: i v(t) u i dν i (t, x) = v(t) dν(t, x).
Note that with this substitution, all measures involved in the measure problem are now positive; there is no need to implement the trick of §3.3. Using these extra constraints, it is now possible to solve bang-bang control problems.

Summary

To summarize, the advantages for introducing the relaxed control problem (6) with measures are the following:

• controls are allowed to be measures with absolutely continuous components and singular components including impulses;

• state constraints are easily handled via support constraints;

• the initial state has a fixed given distribution on some pre-specified domain;

• a free initial state in some pre-specified domain is also allowed.

The moment problem

So far, the hypothesis of polynomial data has not been used, but its crucial importance will appear in this section, where measures will be manipulated through their moments. This will lead to a semi-definite programming (SDP) problem featuring countably many equations.

Define the moments of measure µ as

y µ k = X z k dµ(z). (7)
Then, with a sequence y = (y k ), k ∈ N n , let L y : R[z] → R be the linear functional

f = k f k z k → L y (f ) = k f k y k , f ∈ R[z].
Define the moment matrix of order d ∈ N associated with y as the real symmetric matrix M d (y) whose (i, j)th entry reads

M d (y)[i, j] = L y z i+j = y i+j , ∀i, j ∈ N n d .
Similarly, define the localizing matrix of order d associated with y and h ∈ R[z] as the real symmetric matrix M d (h y) whose (i, j)th entry reads

M d (h y)[i, j] = L y h(z) z i+j = k h k y i+j+k , ∀i, j ∈ N n d .
As a last definition, a sequence y µ = (y µ k ) is said to have a representing measure if there exists a finite Borel measure µ on X, such that relation [START_REF] Gajardo | Minimal time sequential batch reactors with bounded and impulse controls for one or more species[END_REF] holds for every k ∈ N n . Now comes the crucial result of the section: a sequence of moments y µ has a representing measure defined on a semi-algebraic set X µ = {x : Theorem 3.8]. This has the very practical implication that the measure problem defined in (6) has an equivalent formulation in terms of moments. Indeed, because all problem data were assumed to be polynomial, the criterion in (6) can be transformed into a linear combination of moments to be minimized:

a µ i (x) ≥ 0, i = 1, 2, . . .} if and only if M d (y µ ) 0, ∀ d ∈ N and M d (a µ i y µ ) 0, ∀ d ∈ N and ∀a µ i defining set X µ [13,
V m = inf y (b µ ) ′ y µ + (b ν ) ′ y ν + (b µ T ) ′ y µ T = b ′ y (8)
where the infimum is now over the aggregated sequence y of moments of all the measures. Because the test functions were also restricted to be polynomials, the constraints in ( 6) can be turned into countably many linear constraints on the moments:

A µ y µ + A ν y ν + A µ 0 y µ 0 + A µ T y µ T = Ay = 0. ( 9
)
The only non-linear part are the SDP constraints for measure representativeness, to be satisfied ∀d ∈ N:

M d (y µ ) 0, M d (a µ i y µ ) 0, M d (y ν ) 0, M d (a ν i y ν ) 0, M d (y µ 0 ) 0, M d (a µ 0 i y µ 0 ) 0, ( 10
) M d (y µ T ) 0, M d (a µ T i y µ T ) 0.

LMI relaxations

The final step to reach a tractable problem is relatively obvious: we simply truncate the problem to its first few moments. Let d 1 ∈ N be the smallest integer such that all criterion monomials belong to N n+1 2d 1 . This is the degree of the so called first relaxation. For each relaxation, we reach a standard LMI problem that can be solved numerically by off-the-shelf software by simply truncating Eq. ( 8), ( 9) and ( 10) to involve only moments in N n+1 2d , with d ≥ d 1 the relaxation order. Observe that

d j > d i ⇒ V d j M ≥ V d i M .
Therefore, by solving the truncated problem for ever greater relaxation orders, we will obtain a monotonically non-decreasing sequence of lower bounds to the true cost. In the examples below, we will see that in practice, the optimal cost is usually reached after a few relaxations.

Academic examples

In this section, many examples are presented to showcase the different features of the method. Ex. 1 to 5 are variations of the same basic problem to give a thorough tour of the method's capabilities. Ex. 6 is taken from the literature and shows how the method compares to, or rather nicely complements, existing optimal control algorithms. All examples use GloptiPoly [START_REF] Henrion | GloptiPoly 3: Moments, Optimization and Semidefinite Programming[END_REF] for building the truncated LMI moment problems and SeDuMi [START_REF] Sturm | Using SeDuMi 1.02, a Matlab toolbox for optimization over symmetric cones[END_REF] for their numerical solution.

Before proceeding to the examples, define the marginal M d (y, z) of a moment matrix with respect to variable z as the moment matrix of the subsequence of moments concerning polynomials of z only.

Example 1 (Basic impulsive problem).

V = inf u(t) 2 0 x 2 (t) dt
(t) = u(t) x(0) = 1, x(2) = 1 2 x 2 (t) ≤ 1.
In this introductory example, it is straightforward to notice that the optimal solution consists of reaching the turnpike x(t) = 0 by an impulse at initial time t = 0, and likewise, departing from it by an impulse at final time t = T = 2, see Fig. 1.

The associated measure problem reads:

V M = inf µ,ν [0,T ]×X x 2 dµ such that X T v dµ T - X 0 v dµ 0 = [0,T ]×X ∂v ∂t dµ + [0,T ]×X ∂v ∂x dν ∀ v ∈ R[t, x] µ 0 = δ 0 X 0 = {1} µ T = δ 1 2 X T = 1 2 X = x ∈ R : 1 -x 2 ≥ 0 .
Using the procedure outlined above, one obtains a series of truncated moment problems that can be solved by semi-definite programming. Letting y µ ij = t i x j dµ, the first LMI relaxation is Example 3 (Discrete control set with chattering). We take again Ex. 1 with the additional constraint that the control u(t) takes its value in the set U = {±1}, using the method explained in Section 3.4. The solution to this problem is easy to infer: reach the turnpike x(t) = 0 as quickly as possible by applying the negative control until t = 1, then chatter with equal probability to remain on the turnpike until t = 3 2 , after which the positive control must be applied until t = 2 (see Fig. 3). This solution has an optimal cost of 3 8 ≈ 0.375. Compare this value with those of Table 1, which presents the evolution of the criterion with respect to the relaxation order of the truncated problem. After the fourth relaxation, the marginal w.r.t. x of the control measure corresponding to the control u(t) = +1 closely approaches the positive measure 1 2 I [1, 3 2 ] (dx) + I [ 3 2 ,2] (dx) while the marginal w.r.t. x of the control measure of u(t) = -1 converges to I [0,1] (dx) + 1 2 I [1, 3 2 ] (dx), as expected.

V 1 M = inf y y µ 02
Example 4 (Infeasible problem). If the problem is infeasible, it may be detected by the infeasibility of one of the LMI relaxations. Take Ex. 1 with the additional total variation that puts the end point out of reach from the starting point. Indeed, at the first relaxation, the LMI problem is flagged as infeasible with a Farkas dual vector, providing a certificate of infeasibility of the original problem.

Example 5 (Unbounded problem). If the problem is unbounded, it will be detected at the first LMI relaxation. Consider the problem of maximizing the total variation of a linear control problem: sup

u(t) 1 0 |u(t)| dt (11)
such that

ẋ(t) = u(t) x(0) = 0, x(1) = 0 x 2 (t) ≤ 1.
As expected, the LMI problem from the first relaxation on is flagged as unbounded because its dual is flagged as infeasible.

Example 6 (Bang-bang control of the Vanderpol equation). Consider the following timeoptimal problem of the Vanderpol equation: inf

u(t)∈U T such that ẋ1 (t) = x 2 (t) ẋ2 (t) = -x 1 (t) -(x 2 1 (t) -1) x 2 (t) + u(t) x(0) = -0.4 -0.6 ′ , x(T ) = 0.6 0.4 ′ U = {±1} .
In [START_REF] Simakov | Computations for time-optimal bang-bang control using a Lagrangian formulation[END_REF], this problem is solved by applying a gradient-based optimization technique on a parametrization of admissible trajectories, with a minimum time of 2.14. However, this method can only prove the local optimality of solutions. Applying our method, we obtain a cost of 2.15 at the fifth relaxation, providing a (numerical) certificate of global optimality for that local solution.

7 The fuel-optimal linear impulsive guidance rendezvous problem

In this section, the proposed approach is applied to the far-range rendezvous in a linearised gravitational field. This problem is defined as a fixed-time minimum-fuel impulsive orbital transfer between two known circular orbits. Under Keplerian assumptions and for a circular rendezvous, the complete rendezvous problem may be decoupled between the out-of-plane rendezvous problem for which an analytical solution may be found [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] and the coplanar problem. Therefore, only coplanar circular rendezvous problems based on the Hill-Clohessy-Wiltshire equations and associated transition matrix [START_REF] Clohessy | Terminal guidance system for satellite rendezvous[END_REF] are considered for numerical illustration of the proposed results. The general framework of the minimumfuel fixed-time coplanar rendezvous problem in a linear setting is recalled in [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] and [START_REF] Arzelier | Using polynomial optimization to solve the fuel-optimal linear impulsive rendezvous problem[END_REF] where an indirect method based on primer vector theory is proposed. Considering the necessity of easy-to-implement numerical solution for on-board guidance algorithms, direct methods based on linear programming (LP) problem may be used as in [START_REF] Mueller | Collision avoidance manoeuver planning with robust optimization[END_REF]. For an a priori fixed number of impulsive manoeuvres and using a classical transcription method [START_REF] Betts | Practical methods for optimal control using nonlinear programming[END_REF] [14], the genuine infinite-dimensional problem may be converted into a finite-dimensional approximation given by the following LP problem:

V LP = min u N i=1 u θ i 1 s.t. x(θ f ) = Φ(θ f , θ 1 )x(θ 1 ) + N i=1 Φ(θ f , θ i )Bu θ i x(θ 1 ) = x 0 , x(θ f ) = x f ( 12
)
where Φ is the Hill-Clohessy-Wiltshire transition matrix, B = 0 2×2 1 2 ′ and u θ i is the vector of velocity increments at θ i in the local vertical local horizontal (LVLH) frame [START_REF] Arzelier | Using polynomial optimization to solve the fuel-optimal linear impulsive rendezvous problem[END_REF]. Time has been changed to the true anomaly θ for the independent variable as is usual in the literature [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF], and it ranges in the interval [θ 1 , θ f ]. Note that this formulation implies that only the impulsive solution of the general linear rendezvous problem may be obtained for a fixed number of velocity increments.

To be consistent with our previous notations we let t = θ, θ 0 = 0 and θ f = T . Our impulsive optimal control problem (2) writes

V M = inf w(t) T 0 |dw 1 |(t) + |dw 2 |(t) s.t. dx =     0 0 1 0 0 0 0 1 0 0 0 2 0 3 -2 0     x(t)dt +     0 0 0 0 1 0 0 1     dw(t) x(0) = x 0 , x(T ) = x f
where state components model positions (X, Z) = (x 1 , x 2 ) in the orbital plane, and their respective velocities ( Ẋ, Ż) = (x 3 , x 4 ). It clearly encompasses formulation [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF] since it allows to consider continuous or impulsive thrusters as well. In both cases, the fuel consumption is measured by the one-norm of vector function [θ 1 , θ f ] → dw 1 [START_REF] Ross | Space trajectory optimization and L 1 optimal control problems[END_REF] whereas the two-norm of this vector is used in general in the literature, see [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF], [START_REF] Arzelier | Using polynomial optimization to solve the fuel-optimal linear impulsive rendezvous problem[END_REF] and references therein.

For the sake of comparison between these two approaches, two academic examples taken from [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] are presented.

Example 7 (In-plane rendezvous 1). Consider the first case presented in [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF]. It consists of a coplanar circle-to-circle rendezvous with zero eccentricity. The rendezvous manoeuvre must be completed in one orbital period with boundary conditions x 0 = [ 1 0 0 0 ] ′ and x f = [ 0 0 0 0 ] ′ . This type of rendezvous is usually difficult to handle by numerical methods because of its singularity due to the high number of symmetries involved.

With a grid of N = 50 points, the LP algorithm gives a two-impulse solution at the initial and final times of the rendezvous without interior impulse nor initial or final coasting period. The optimal impulses are both horizontal and opposite u 0 = -u 2π = 0.05305 0 ′ . The fuel cost is given by V LP = 0.1061. The LMI method has no difficulty to recover the optimal solution given by the LP algorithm. A cost of 0.1061 is obtained for each relaxation. It is then easy to extract from the matrices that the optimal solution for the first control consists of two symmetric impulses of magnitude 0.0531 at the initial and final times, while the second control is identically 0. The optimal trajectory in the orbital plane is depicted in Figure 4 where + indicates the 50 points of discretization. Example 8 (In-plane rendez-vous 2). As a second example, the third case of [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] is revisited. The rendezvous is nearly identical to the previous one except for the final condition that imposes to reach the target with relative velocity of 0.427 in the Z direction, namely x 0 = [ 1 0 0 0 ] ′ and x f = [ 0 0 0 0.427 ] ′ . Again, a grid of N = 50 points is used when running the LP algorithm. It converges to a four-impulse trajectory depicted in Fig. 6. The numerical results are summarized in Table 3. Using our algorithm, we reached the same criterion (within numerical tolerance) after the fourth relaxation (see Tab. 3). As usual, the controls can be inferred from the moment matrix of the ν measures. Indeed, ν 1 converges to the measure (u θ i ) 1 δ θ i with impulse amplitudes (u θ i ) 1 and anomaly θ i taken from Table 2, while ν 2 converges to an all zero measure. Not only does this result prove the global optimality of the conjectured solution within the class of all impulsive solutions no matter the number of impulses, but it also shows that it is optimal over all measure thrust solutions.

Conclusion

The focus of this work is on actual computation of optimal impulsive controls for systems described by ordinary differential equations with polynomial dynamics and polynomial (semialgebraic) constraints on the state. State trajectory and controls are measures which are linearly constrained, resulting in an infinite-dimensional linear programming (LP) problem consistent with the formalism of our GloptiPoly software [START_REF] Henrion | GloptiPoly 3: Moments, Optimization and Semidefinite Programming[END_REF]. This LP problem on measures can then be solved numerically via a hierarchy of linear matrix inequality (LMI) relaxations, for which off-the-shelf semi-definite programming (SDP) solvers can be used. The optimal impulse sequence can then be retrieved by simple linear algebra, and global optimality can be verified by a posteriori simulation or comparison with suboptimal control sequences computed by alternative techniques.

For space rendezvous, our technique can be readily adapted to cope with state (e.g. obstacle avoidance) constraints, as soon as they are basic semialgebraic. Other criteria than the total variation can also be handled. Smoother solutions can be expected, maybe consisting of a mix of absolutely continuous and singular controls, including impulsive controls.

Figure 1 :

 1 Figure 1: Trajectory for Ex. 1

Figure 2 : 2 Figure 3 :

 223 Figure 2: Trajectory for Ex. 2

Figure 4 :

 4 Figure4: Trajectory in the orbital plane (X, Z) in LVLH: Case 1 of[START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF]

Figure 5

 5 Figure 5 shows position, velocity and impulses history versus true anomaly.

Figure 5 :

 5 Figure 5: Positions (X solid, Z dashed), velocities ( Ẋ solid, Ż dashed) and impulses (on X axis): Case 1 of [5]

Figure 6 :

 6 Figure 6: Trajectory in the orbital plane (X, Z) in LVLH: Case 3 of [5]

Table 3 :

 3 Criterion as a function of LMI relaxation order for Ex. , velocity and impulses history are illustrated in Figure7. Note the symmetry of the optimal four-impulse solution.

Figure 7 :

 7 Figure 7: Positions (X solid, Z dashed), velocities ( Ẋ solid, Ż dashed) and impulses (on X axis): Case 3 of [5]

Table 1 :

 1 Criterion as a function of LMI relaxation order for Ex. 3

		d	1	2	3	4
		V d M	0.000 0.288 0.368 0.372
	constraint	2 0 |u(t)| dt ≤ 1 4

Table 2 :

 2 Impulse times and amplitudes for Ex. 8

		LMI method		LP method
	θ i	(u θ i ) 1	(u θ i ) 2	θ i	(u θ i ) 1	(u θ i ) 2
	0	-0.0386	0	0	-0.0392	0
	1.791 +0.109	0	1.795 +0.109	0
	4.495 -0.109	0	4.488 -0.109	0
	6.283 +0.0389	0	6.283 +0.0392	0

CNRS; LAAS; 7 avenue du colonel Roche, F-31077 Toulouse; France.

Université de Toulouse; UPS, INSA, INP, ISAE; UT1, UTM, LAAS; F-31077 Toulouse; France

Faculty of Electrical Engineering, Czech Technical University in Prague, Technická 2, CZ-16626 Prague, Czech Republic

Institut de Mathématiques de Toulouse, Université de Toulouse; UPS; F-31062 Toulouse, France.

We are grateful to Térence Bayen for pointing out this reference to us.

In all rigour, it could be possible to include state jumps in G, but this requires a careful definition of what is meant by integration, as done e.g. for studying stochastic differential equations. This goes well beyond the scope of this paper.

We suspect however that for the control problems studied in this paper, subset J can be assumed countable without loss of generality.

subject to the linear constraints associated to the dynamics: y µ T 00y µ 0 00 = 0 y µ T 10y µ 0 10 = y µ 00 y µ T 01y µ 0 01 = y ν + 00y ν -

It turns out that the optimal value V M = 0 is estimated correctly (to numerical tolerance) from the first relaxation on and that the optimal trajectory x(t) = 0 can easily be recovered. Indeed, the marginal M d (y µ , x) is the length of the time interval multiplying a truncated moment matrix of a Dirac measure concentrated at x = 0, while its marginal with respect to t equals a truncated Lebesgue moment matrix on the [0, 2] interval. More importantly, one can recover the optimal controls as the marginal M d (y ν , t) is the weighted sum of Dirac measures located at the impulse times, the weights being the impulse amplitudes. In summary, we can recover numerically the optimal measures µ(dt, dx) = I [0,2] (dt)δ 0 (dx), ν(dt, dx) = -δ 0 (dt)I [0,1] (dx) + δ 2 (dt)I [0, 1 2 ] (dx). Example 2 (Total variation constraints). We take back Ex. 1 with an additional constraint on the total variation of the control:

Clearly, the solution of Ex. 1, with a total variation of 3 2 , violates this extra constraint, so the algorithm should converge to another solution. Again, from the first relaxation on, the cost of the associated truncated moment problem is 1 8 . It is also plain to see that M d (y µ , x) is the truncated moment matrix of a Dirac located at x = 1 4 , hinting a trajectory x(t) = 1 4 . On the control side, starting from the second relaxation, it also becomes evident that M d (y ν , t) is the truncated moment matrix of the signed measure -3 4 δ 0 + 1 4 δ 2 , revealing the times and amplitudes of impulses compatible with admissible controls. This leads to the trajectory of Fig. 2, which therefore is an optimal solution of the problem.


Thomas Jacob

John Wahr

Tad Pfeffer

Sean Swenson

Recent contributions of glaciers and ice caps to sea level rise, from GRACE

Glaciers and ice caps (GIC) are important contributors to present-day global mean sea-level rise (SLR). Most previous global mass balance estimates for GIC rely on interpolation of sparse mass balance measurements, representing a small fraction of the GIC area. Instead, we here perform a global, simultaneous inversion of mass change over all ice-covered regions larger than 100 km 2 using monthly GRACE-derived satellite gravity fields spanning January 2003 to December 2010. This is the first GRACE-based study where every such GIC region is considered and the results quantified. We conclude that GIC, excluding the Greenland and Antarctic peripheral glaciers and ice caps (PGIC), lost mass at a rate of 148 ± 30 giga tonnes per year (Gt/yr) during this period, contributing 0.41 ± 0.08 mm/yr to SLR. This rate is significantly smaller than previous estimates that rely on extrapolation of mass balance measurements. Our high mountain Asia (HMA) result, in particular, shows only a modest mass loss; far smaller than previously published estimates. For completeness, we also consider Greenland and Antarctica for this same time period, and find that those ice sheets, including their PGIC lost 384 ± 71 Gt/yr, contributing 1.06 ± 0.19 mm/yr to SLR. The total SLR contribution from all ice-covered regions was then 1.48 ± 0.26 mm/yr, agreeing with independent estimates of SLR from new water to within the respective error bars.


Yvon Maday

AND Julien Salomon

Kamel Riahi
email: riahi@ann.jussieu.fr

AN INTERMEDIATE TARGETS METHOD FOR TIME PARALLELIZATION IN OPTIMAL CONTROL

In this paper, we present a method that enables to solve in parallel the Euler-Lagrange system associated with the optimal control of a parabolic equation. Our approach is based on an iterative update of a sequence of intermediate targets and gives rise independent sub-problems that can be solved in parallel. Numerical experiments show the efficiency of our method.

Dans cet article, on présente une méthode permettant une parallélisation en temps de la résolution des équations d'Euler-Lagrange associées à un problème de contrôle optimal dans le cas parabolique. Notre approche est basée sur une mise à jour itérative de cibles intermédiaires et donne lieu à des sous-problèmes de contrôle indépendants. Les résultats numériques prouvent l'efficacité de la méthode.

Introduction:

In the last decade, time domain decomposition has been exploited to accelerate the simulation of systems ruled by time dependent partial differential equations [START_REF] Lions | Virtual and effective control for distributed systems and decomposition of everything[END_REF]. Among others, the parareal algorithm [START_REF] Lions | Résolution d'EDP par un shéma pararréel[END_REF] or multi-shooting schemes [START_REF] Bellen | Parallel algorithms for initial value problems for nonlinear vector difference and differential equations[END_REF] have shown excellent results. In the framework of optimal control, this approach has been used to control parabolic systems [START_REF] Maday | A parareal in time procedure for the control of partial differential equations[END_REF], [START_REF] Mathew | Analysis of block parareal preconditioners for parabolic optimal control problems[END_REF]. In this paper, we introduce a new approach to tackle such problems. The strategy we follow is based on the concept of target trajectory that has been introduced in the case of hyperbolic systems in [START_REF] Maday | Parareal in time control for quantum systems[END_REF]. Because of the irreversibility of parabolic equations, a new definition of this trajectory is considered. It enables us to define at each bound of the time sub-domains relevant initial conditions and intermediate targets, so that the initial problem is split into independent optimization problems.

We now introduce some notations. Given d ∈ N, we consider the optimal control problem associated with a heat equation defined on a compact set Ω ⊂ R d and a time interval interval I = [0, T ]. The control applies on a subset Ω c ⊂ Ω and belongs to the Hilbert space H := L 2 (I; L 2 (Ω c )) whose scalar product and norm are denoted by ., . H and . H respectively. We also denote by . and . c the norms associated with the spaces L 2 (Ω) and L 2 (Ω c ) respectively. Given a function ϕ defined on I × Ω, we denote its restriction to a sub-domain

I ′ × Ω ⊂ I × Ω by ϕ | I ′ .
The paper is organized as follows. The optimal control problem is presented in Sec. 2 ; Sec. 3 is devoted to the description of the time parallelization setting. The algorithm is given in Sec. 4. Finally, numerical results are presented in Sec. 5.

Optimal control problem

Given α > 0, T > 0, a target state y target ∈ L 2 (Ω), we consider the optimal control problem :

Find v ⋆ ∈ H v ⋆ := argmin v∈H J(v),
where J is the quadratic cost functional defined by ( 1)

J(v) := 1 2 y(T ) -y target 2 + α 2 v 2
H . The state variable y depends linearly on the control v through the evolution equation [START_REF] Maday | A parareal in time procedure for the control of partial differential equations[END_REF] ∂ t y -ν∆y = Bv, with an initial data y 0 ∈ L 2 (Ω). In this equation ∆ denotes the Laplace operator, ν is the diffusion coefficient and B is the injection from

L 2 (Ω c ) to L 2 (Ω), so that Bv ∈ L 2 (I; L 2 (Ω)).
The corresponding optimality system reads

∂ t y -ν∆y = Bv on I × Ω y(0) = y 0 , (3)
∂ t p + ν∆p = 0 on I × Ω p(T ) = y(T ) -y target , (4)
αv + B * p = 0, (5)
where B * is the adjoint operator of B. Note that for any α > 0, the functional J is strictly convex, so that the system (3-5) has a unique solution v ⋆ . We denote by y ⋆ , p ⋆ the associated state and adjoint state.

Time parallelization setting

We now aim at solving in parallel the coupled system corresponding to Equations [START_REF] Bal | A parareal time discretization for non-linear PDEs with application to the pricing of an american put[END_REF][START_REF] Bellen | Parallel algorithms for initial value problems for nonlinear vector difference and differential equations[END_REF][START_REF] Lions | Résolution d'EDP par un shéma pararréel[END_REF]. To do this, we decompose the interval I into subintervals and introduce a set of intermediate target states so that the initial problem is replaced by a set of smaller independent optimal control problems. The resolution of these problems is achieved using an inner loop, whereas the intermediate states are updated by an outer loop. We start with the definition of the target states. Given a control v ∈ H and its corresponding state y and adjoint state p, we define the target trajectory by: ( 6)

χ(v) = y(v) -p(v) on I × Ω
In what follows and for the sake of simplicity, we omit the dependence of y, p and χ on the control v in the notations. The introduction of this trajectory is motivated by the following result.

Lemma 1. We keep the previous notations. Let τ ∈]0, T [, and the optimal control problem: Find w ⋆ τ ∈ H such that w ⋆ τ := argmin w∈H J τ (w), where

J τ (w) := 1 2 y(τ ) -χ ⋆ (τ ) 2 + α 2 w 2 L 2 ([0,τ ];L 2 (Ωc))
with y(τ ) the solution of the Equation (2). We have

w ⋆ τ = v ⋆ | [0,τ ]
Proof: Thanks to the uniqueness of the solution of the optimization problem associated to J τ , it is sufficient to show that v ⋆ | [0,τ ] is a solution of its optimality system. This one is given by Equation (3) (with v = w τ ) restricted to [0, τ ] × Ω and ( 7)

∂ t p + ν∆p = 0 on [0, τ ] × Ω p(τ ) = y(τ ) -χ(τ ), (8)
αw τ + B * p = 0. First, note that y ⋆ | [0,τ ] obviously satisfies Equation (3) restricted to [0, τ ] × Ω with v = v ⋆ | [0,τ ] .
It directly follows from the definition of χ ⋆ (see ( 6)), that: [START_REF] Mathew | Analysis of block parareal preconditioners for parabolic optimal control problems[END_REF]. Finally, Equation ( 8) is a consequence of ( 5). The result follows.

p ⋆ (τ ) = y ⋆ (τ ) -χ ⋆ (τ ), so that p ⋆ | [0,τ ] satisfies
Given N ≥ 1, we decompose the interval I = [0, T ] into subintervals I = ∪ N -1 n=0 I n , where I n = [t n , t n+1 ], t 0 = 0 < t 1 < ... < t N -1 < t N = T . We also introduce the spaces H n := L 2 (I n ; L 2 (Ω c )) and the corresponding scalar product ., . Hn and norm . Hn . In this framework, given v n ∈ H n we define v ⋆ n as follows (9)

v ⋆ n := argmin vn∈Hn J v n (v n ), with (10) J v n (v n ) := 1 2 y n (t n+1 ) -χ(t n+1 ) 2 + α 2 v n 2 Hn ,
where χ is associated to v through the definition [START_REF] Maday | Parareal in time control for quantum systems[END_REF]. In this functional, the state y n is defined by

(11) ∂ t y n -ν∆y n = Bv n on I n × Ω y n (t n ) = y(t n ).
These subproblems have the same structure as the original one and are also strictly convex. Note also that their definitions depend on the control v through the target trajectory, hence the notation J v n . The optimality system associated with these minimization problems are given by Equation ( 11) and ( 12)

∂ t p n + ν∆p n = 0 on I n × Ω p n (t n+1 ) = y(t n+1 ) -χ(t n+1 ), ( 13
) αv n + B * p n = 0.
Lemma 2. We keep the previous notations. Denote by χ ⋆ the target trajectory defined by Equation (6) with y = y ⋆ and p = p ⋆ and by y ⋆ n , p ⋆ n , v ⋆ n the solutions of Equations (11-13) associated with v ⋆ . One has:

v ⋆ n = v ⋆ |I n .
The proof of this result follows the lines of Lemma 1 and is left as an exercise to reader.

Algorithm

We are now in the position to propose a time parallelized procedure to solve Equations ( 3)- [START_REF] Lions | Résolution d'EDP par un shéma pararréel[END_REF]. Consider an initial control v 0 and assume that, at step k, a control v k is known. The computation of v k+1 is achieved as follows:

(1) Compute y k , p k and the associated target trajectory χ k according to Equations (3), ( 4) and ( 6) respectively. (2) For n = 0, ..., N -1, solve the sub-problems (9) in parallel and denote by ṽk+1 n the corresponding solutions. (3) Define ṽk+1 as the concatenation of the sequence (ṽ k+1 n ) n=0,...,N -1 . (4) Update the control variable by ( 14)

v k+1 = v k + θ k (ṽ k+1 -v k ),
where the value θ k is chosen to minimizes J(v k + θ k (ṽ k+1v k )). We have not detailed Step 2 as we rather aim at presenting a structure for a general approach. However, because of the strict convexity of the problems we consider, a small number of conjugate gradient method steps can be used to achieve the resolution of these steps.

Numerical experiments

In this section, we test the efficiency of our method. We consider a 2D example, where Ω = [0, 1] × [0, 1] and

Ω c = [ 1 3 , 2 3 ] × [ 1 3 , 2 3
] . The parameters related to our control problem are T = 6.4, α = 10 -2 and ν = 10 -2 . The time interval is discretized using a uniform step δt = 10 -3 , and an Implict-Euler solver is used to approximate the solution of Equations [START_REF] Bal | A parareal time discretization for non-linear PDEs with application to the pricing of an american put[END_REF][START_REF] Bellen | Parallel algorithms for initial value problems for nonlinear vector difference and differential equations[END_REF]. For the space discretization, we use P 1 finite elements. Our implementation make use of the freeware FreeFem [START_REF] Pironneau | Free soft : FreeFem++-mpi[END_REF] and the parallelization is achieved thanks to the Message Passing Interface library. The independent optimization procedures required in Step 2 are simply carried out using one step of an optimal step gradient method. The results are presented in Figure 1.

In the first plot, we consider the evolution of the cost functional values with respect to the iterations and do not take into account the parallelization. The result reveals that our algorithm significantly accelerates the optimization process. This outcome may indicates that the splitting introduced in our approach acts as a preconditionner during the numerical optimization. This will be the purpose of some further investigation [9], in the same spirit as in [START_REF] Mathew | Analysis of block parareal preconditioners for parabolic optimal control problems[END_REF].

In a second plot, we represent the evolution of the cost functional values with respect to the number of matrix-vector product. Parallel computations that are done in Step 2 are only counted once. When comparing with a standard optimal gradient step method, we observe speed-up approximatively equal to 3.